<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">%~MoulinÃ© par MaN_auto v.0.32.4 2024-08-22 11:39:02
\documentclass[CRPHYS, Unicode, XML, screen, biblatex]{CEDRAM}

\addbibresource{crphys20240229.bib}


\usepackage{xcolor}
\usepackage{soul}
\usepackage{bm}
%\usepackage{ulem}
\usepackage[squaren, Gray, cdot]{SIunits}

\newcommand{\deriv}[2]{\frac{\dd #1}{\dd #2}} 
\newcommand{\dderiv}[2]{\frac{\dd^2 #1}{{\dd #2}^2}} 
\newcommand{\dpart}[2]{\frac{\partial #1}{\partial #2}} 
\newcommand{\ddpart}[2]{\frac{\partial^2 #1}{{\partial #2}^2}}
\newcommand{\ud}{\,\mathrm{d}}
\newcommand{\dd}{\mathrm{d}} 

\newcommand{\QBO}{\mathrm{QBO}}
\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bbC}{\mathbb{C}}
%\newcommand{\bbZ}{\mathbb{Z}}
\newcommand{\clL}{\mathcal{L}}
\newcommand{\clH}{\mathcal{H}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

\let\oldtilde\tilde
\renewcommand*{\tilde}[1]{\mathchoice{\widetilde{#1}}{\widetilde{#1}}{\oldtilde{#1}}{\oldtilde{#1}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\CDRsetmeta{articletype}{review}

%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{Geophysical and astrophysical fluid dynamics in the Laboratory}
%\def\cdr@specialtitle@french{Dynamique des fluides gÃ©ophysiques et astrophysiques au laboratoire}
%\makeatother{}


\title{Quasi-biennial oscillation: laboratory experiments}

\alttitle{Oscillation quasi-biennale : expÃ©riences de laboratoire}

\author{\firstname{BenoÃ®t} \lastname{Semin}\IsCorresp}
\address{PMMH, CNRS, ESPCI Paris, PSL University, Sorbonne UniversitÃ©, UniversitÃ© Paris-CitÃ©, F-75005, Paris, France}
\email[B. Semin]{benoit.semin@espci.fr}

\author{\firstname{Fran\c cois}
\lastname{P\'{e}trelis}}
\address{LPENS, CNRS, ENS Paris, PSL University, Sorbonne UniversitÃ©, UniversitÃ© Paris-CitÃ©, F-75005, Paris, France}
\email[F. Petrelis]{petrelis@lps.ens.fr}



%\newcommand{\deriv}[2]{\frac{\dd #1}{\dd #2}}
%\newcommand{\dderiv}[2]{\frac{\dd^2 #1}{{\dd #2}^2}}
%\newcommand{\dpart}[2]{\frac{\partial #1}{\partial #2}}
%\newcommand{\ddpart}[2]{\frac{\partial^2 #1}{{\partial #2}^2}}
%\newcommand{\ud}{\,\dd}
%\newcommand{\dd}{\dd}

\keywords{
\kwd{Fluid mechanics}
\kwd{internal gravity waves}
\kwd{wave-mean flow interaction}
\kwd{instability}
\kwd{non-linear physics}}

\altkeywords{
\kwd{MÃ©canique des fluides} 
\kwd{ondes internes de gravitÃ©}
\kwd{interaction ondes-Ã©coulement moyen} \kwd{instabilitÃ©}
\kwd{physique non-linÃ©aire}}

\begin{abstract}
The quasi-biennial oscillation (QBO) is an oscillation of the wind in the equatorial stratosphere. This wind is a mean flow induced by atmospheric waves, including internal gravity waves, which explain that the period (28 month) is not linked to any astrophysical forcing. This oscillation has only been reproduced in 3 laboratory experiments, which share a similar geometry. We present the details of our experimental set-up, and we explain which improvements allowed us to obtain quantitative measurements during long times. We show experimentally the feedback of the mean flow on the waves, which is one of the key ingredient of the oscillation. The details of the analytical resolution of the 1D model of Plumb and McEwan are given. We compare experimental results to analytical and numerical results, and found a qualitative agreement. The period decreases when the forcing increases, and the amplitude of the mean flow is not monotonic with respect to height and displays two local maxima as a function of height close to the threshold. The bifurcation is always a Hopf one, but can be subcritical or supercritical depending on the dominant dissipation mechanism of the mean flow which can be tuned experimentally by changing the Brunt--V{\"a}is{\"a}l{\"a} frequency. We argue that an investigation of the bifurcation in general circulation models (GCM) is of interest to better understand the evolution of the QBO due to climate change.
\end{abstract}

\begin{altabstract}
L'oscillation quasi-biennale (QBO) est une oscillation du vent dans la stratosphÃ¨re Ã©quatoriale. Ce vent est un Ã©coulement moyen engendrÃ© par des ondes atmosphÃ©riques, notamment des ondes internes de gravitÃ©, ce qui explique que la pÃ©riode (28 mois) n'est liÃ©e Ã&nbsp; aucun forÃ§age astrophysique. Cette oscillation n'a Ã©tÃ© reproduite qu'Ã&nbsp; 3 reprises par des expÃ©riences analogues de laboratoire, et ces 3 expÃ©riences ont Ã©tÃ© rÃ©alisÃ©es dans la mÃªme gÃ©omÃ©trie. Nous prÃ©sentons en dÃ©tail notre dispositif expÃ©rimental, et en particulier les amÃ©liorations qui nous ont permis d'obtenir des mesures quantitatives sur de longues durÃ©es. Nous montrons expÃ©rimentalement la rÃ©troaction du l'Ã©coulement moyen sur les ondes, qui est l'un des ingrÃ©dients clÃ©s de la QBO. Nous donnons les dÃ©tails de la rÃ©solution analytique du modÃ¨le 1D de Plumb et McEwan. Nous comparons les rÃ©sultats expÃ©rimentaux Ã&nbsp; ceux obtenus analytiquement et numÃ©riquement, et nous montrons qu'ils sont en accord qualitatif. La pÃ©riode de l'Ã©coulement moyen diminue lorsque le forÃ§age augmente, et l'amplitude d'Ã©coulement moyen en fonction de la hauteur n'est pas monotone et prÃ©sente au voisinage du seuil deux maxima locaux. La bifurcation est toujours une bifurcation de Hopf, mais elle peut Ãªtre sous-critique ou supercritique en fonction du mÃ©canisme dominant de dissipation de l'Ã©coulement moyen, qui peut Ãªtre changÃ© expÃ©rimentalement en modifiant la frÃ©quence de Brunt-V{\"a}is{\"a}l{\"a}. Nous concluons en expliquant qu'une Ã©tude de la bifurcation dans les modÃ¨les de circulation gÃ©nÃ©rale (GCM) serait intÃ©ressante pour mieux comprendre l'Ã©volution de la QBO liÃ©e au changement climatique.
\end{altabstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\dateposted{2024-09-30}
\begin{document}
\maketitle



\section{Introduction}

The quasi-biennial oscillation (QBO) is the periodic reversal of the wind in the equatorial stratosphere~\cite{Baldwin2001} (the stratosphere is the second layer of the atmosphere, between $\unit{16}{\kilo\meter}$ to $\unit{50}{\kilo\meter}$ at the equator). The period of the oscillation is $28$ months on average, and is not directly linked to any astrophysical forcing. This wind is known to be generated by a broad spectrum of atmospheric waves, including small scale internal gravity waves and planetary waves like Kelvin waves, Rossby-gravity waves and inertia-gravity waves~\cite{Baldwin2001}. The interaction between the wind and the waves leads to the periodic reversal of the wind.

A key feature of the QBO is its downward propagation without loss of amplitude~\cite{Baldwin2001}. The position where the mean flow changes sign drifts towards the location from where the waves are emitted, i.e. the troposphere (the troposphere is the bottom layer of the atmosphere, between $\unit{0}{\kilo\meter}$ to $\unit{16}{\kilo\meter}$ at the equator).

The QBO has been mainly studied through observations on Earth, but similar wind oscillations exist on other planets like Saturn~\cite{Fouchet2008} or Jupiter~\cite{LiX2000,Showman_2019}. Internal waves can also be sustained in stably stratified layers in stars~\cite{EJK}, and models of the Sun tachocline also observe periodic reversals of a mean flow similarly to the QBO mechanism~\cite{PINCON}. Overall this mechanism is robust and generic.

The QBO used to be one of the most predictable fluctuations of the large-scale circulation of the Earth's atmosphere~\cite{Anstey2022}. However, two disruptions of the QBO have been observed recently in 2015/16 and 2019/20, after more than 50 years without observed disruption~\cite{Wang2023}. These disruptions are thought to be due to strong extra-equatorial waves forcing.

A key interest of the QBO is the correlation between phenomena outside the tropical stratosphere, the so-called teleconnections~\cite{Butchart2022,Anstey2022}. The QBO influences the convection in the troposphere close to the equator~\cite{Collimore2023} and modulates the mixing and the transport of chemical constituents in the atmosphere, including ozone among others~\cite{Bowman1989}. Even if this wind oscillation is located close to the equator, it has some influence in the extratropics. This wind oscillation is correlated with hurricanes in North America~\cite{Baldwin2001} and affects winter conditions in Europe~\cite{Marshall2009}. It also influences the Madden-Julian oscillation~\cite{Martin2021,Hood2023}, an eastward propagation of large regions of both enhanced and suppressed tropical rainfall, observed over the Indian and Pacific oceans.


The QBO has been simulated using general circulation models (GCM) for more than 20 years now~\cite{Takahashi1999}. It used to be difficult for the models used for climate change forecast to produce a QBO. For instance, most of the models included in CMIP5 (Coupled Model Intercomparison Project Phase 5) do not reproduce the QBO~\cite{Butchart2018}. Now, many climate models simulate QBO-like oscillations, but often with an amplitude weaker than the observed one in the lower stratosphere~\cite{Anstey2022}. Moreover, these models rely on the parametrization of small-scale tropical waves which cannot be resolved~\cite{lott, Anstey2022}, and the uncertainty on the response of the QBO to climate change is at least partly due to uncertainties on internal gravity waves~\cite{Achatz2024}.

Simple 1D models have been developed when the basics mechanisms of the QBO have been understood~\cite{Lindzen1968, Plumb1977}. These models have inspired the analog experiments of the QBO (see~\cite{Plumb1978} and discussion below), and reciprocally the experiments have shown the relevance of the physical ingredients used in these simple models. We will discuss this aspect later in the present article. These simple models are still useful to extract some basics mechanisms, as for example in~\cite{Renaud2019} where the successive bifurcations of the Plumb and McEwan model are investigated and where the existence of a second bifurcation is shown, leading to quasi-periodic trajectories. In~\cite{Renaud2020} the effect of the boundary condition is investigated. The excitation with a spectrum of waves has been studied in~\cite{Leard2020b}, where the authors show that the spreading of the wave energy across a wide frequency range leads to more regular oscillations.

Couston et al. 2018~\cite{Couston2018} have shown that an oscillating flow can emerge from internal gravity waves generated by turbulent convective motion, using direct numerical simulations (DNS), i.e. without parametrization. A 2D model coupling the troposphere and the stratosphere and taking into account the phase change of water has been used to investigate the effect of the QBO on precipitations~\cite{Bui2017}. A comparison between 2D and 3D direct numerical simulations reproducing the experiments has been performed for different forcing parameters~\cite{Wedi2006}.



Three analog experiments of the QBO have been performed: the first has been performed by Plumb and McEwan in 1978~\cite{Plumb1978}, the second by Otobe et al. in 1998~\cite{Otobe1998} and the third by Semin et al. in 2018~\cite{Semin2018}. These $3$ experiments share a similar geometry: a linearly stratified aqueous solution is placed between two concentric vertical cylinders, of size about $0.5$~m. The forcing is achieved using $16$ membranes placed at the bottom~\cite{Plumb1978} or top~\cite{Otobe1998,Semin2018} of the fluid. Each membrane oscillates in opposite phase with its two neighbors, so that the forcing is a pattern of azimuthal standing waves, which is equivalent to two counter-propagating waves. The choice of the geometry, size and number of membranes has thus been done by Plumb and McEwan in 1978~\cite{Plumb1978}.

In the first experiment by Plumb and McEwan~\cite{Plumb1978}, the flow was measured using tracers. Since this study dates back before particle image velocimetry (PIV) became widely available, the particles were tracked manually from ``Cine film''. The average velocity on a cycle was measured. Tracers generated by electrolysis on a vertical wire were also used. This first experiment has validated the wave/mean flow theory, the instability of mean flow, the reversal of the mean flow and the origin of the downward propagation. This was an essential validation of a mechanism which is far from being obvious. A stability diagram is displayed in this article, which shows a very qualitative agreement with the theory. Only measurements over one period are presented.

The second experiment by Otobe et al.~\cite{Otobe1998} has confirmed the results of Plumb and McEwan~\cite{Plumb1978}. The mean flow is measured with particles, and the waves with a MoirÃ© technique. The distortion of the waves by the mean flow is presented. The difficulty to obtain reproducible results is emphasized: the authors state for example ``The period of the oscillation is $1$ to $2$ hours, but there is no apparent correlation with the experimental parameters''.

The third experiment has been done by Semin et al.~\cite{Semin2018}. We will present this experiment in more detail in the present article. This is a quantitative experiment, where both the waves and the mean flow are measured using PIV. Continuous slow injection of new (and heavier) fluid at the bottom of the tank and slow withdrawal of fluid at the top prevented the weakening of the stratification at the top where the membranes move, and thus allows very long measurements. We have determined the nature of the bifurcation. The bifurcation is always a Hopf bifurcation, but can be supercritical or subcritical depending on the main dissipation of the mean flow.

An experiment close to the analog QBO experiments has been performed by LÃ©ard et al.~\cite{Leard2020} (see also~\cite{LeBars2020}). The authors take advantage of the non-monotonic change of the density of water around $\unit{4}{\celsius}$ to excite directly the waves in a stratified layer by the convection in the lower layer. A mean flow appears in the upper stratified layer, but with a dynamics different from the QBO, explained by the authors by the difference in Prandtl numbers between water and air. Another experiment has been reported by LÃ©ard et al. 2021~\cite{Leard2021}. The waves in a stratified fluid are forced by $12$ jets in a upper unstratified layer. A mean flow of very low intensity is obtained in this set-up where the forcing is turbulent and statistically symmetric, contrary to the previously reported experiments of mean flow generation by internal waves~\cite{King2009, Bordes2012}. However, no reversal has been obtained, probably due at least partly to the mixing at the interface between the two layers. These experiments justify a posteriori the use of membranes to force the fluid in the analog QBO experiments.


In this article, we present in Section~\ref{Experimental set-up} the details of our experimental set-up, including many tricks not discussed in our previous articles in~\cite{Semin2016,Semin2018}. The details of the analytical resolution of the Plumb and McEwan model is given in Section~\ref{Theoretical approach} and in Appendix~\ref{Linear and weakly nonlinear solution}. Numerical simulations are also discussed. The experimental results are presented in Section~\ref{Experimental results}.



\section{Experimental set-up}\label{Experimental set-up}



\begin{figure}[!htbp]
\centering
\includegraphics[width=0.75\textwidth]{fig1.pdf}
\caption{Schematic view of the experimental set-up.}
\label{fig:schema_Plumb}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=0.4\textwidth]{fig2L.pdf}
\includegraphics[width=0.35\textwidth]{fig2R.pdf}
\caption{Left: picture of the experimental set-up. The green arrows indicate the motion of the membranes, and the red and blue arrow the direction of the mean flow. Right: picture of a membrane.}
\label{fig:photos}
\end{figure}


A schematic view of the experimental set-up is shown in Figure~\ref{fig:schema_Plumb} and a picture in Figure~\ref{fig:photos} left. It has been already described in~\cite{Semin2016}. It is made of two transparent concentric vertical cylinders. The outer diameter of the inner cylinder is $\unit{365}{\milli\meter}$, the inner diameter of the outer cylinder is $\unit{600}{\milli\meter}$, and the height is $H=\unit{410}{\milli\meter}$.

The gap $h$ between these two cylinders is filled with a linearly density-stratified solution of NaCl or MgCl$_2$ in water, obtained using the double-bucket method~\cite{Otobe1998}. We use NaCl for water softener (Laugil) and magnesium chloride ice melt (France DÃ©neigement). For NaCl, we use $\unit{62}{\liter}$ of tap water in bucket $1$, and $\unit{48.5}{\liter}$ of tap water and $\unit{15}{\kilo\gram}$ of NaCl in bucket $2$. For MgCl$_2$, we use $\unit{55.9}{\liter}$ of tap water in bucket $1$, and $\unit{19.5}{\liter}$ of tap water and $\unit{36.4}{\kilo\gram}$ of MgCl$_2$ in bucket $2$ (MgCl$_2$ is already partially hydrated when we buy it). Moreover $\unit{40}{\milli\liter}$ of chlorine bleach is added in each bucket to avoid algea proliferation. During the injection, bucket $1$ is constantly mixed using a VWR VOS overhead stirrer, as required by the double-bucket method. The density profile is obtained from conductivity, measured using a conductivity probe (Mettler-Toledo SevenGo Pro) positioned at different heights using a linear manual translation stage. We have used a pycnometer to measure the density of solution that we used to calibrate the conductivity. Note that the conductivity of MgCl$_2$ is not monotonic with respect to density: a given value of conductivity may correspond to two values of density. We use the stably stratified property of the fluid to select the right value. The stratification is quantified by the Brunt--V\"{a}is\"{a}l\"{a} frequency, $N=\sqrt {-(g/\rho_0) \, \dd \rho_0 /\dd z}$, where $g$ is the acceleration of gravity, $\rho_0$ is the background density at the top and $z$ is the vertical axis. The typical value of $N$ is $\unit{1.5-2.2}{\radian \cdot \reciprocal \second}$, the highest values are obtained using MgCl$_2$.

The fluid motion is forced using $16$ silicone membranes (Profiles Market, thickness $\unit{0.3}{\milli\meter}$). The membranes are in contact with the top of the fluid and are fixed using nylon and stainless steel screws on a laser-cut PMMA annulus (see Figure~\ref{fig:photos} right). We used Boehm hollow punches to make the fixation holes in the membranes. However, it happened several times that a membrane broke due to tearing, which then induced mixing at the top and unreliable experimental results, which were then discarded. The silicone membranes do not degrade in salty water. A gap of a few millimeters is located between the annulus and the cylinders, so that air can be removed from the membrane by pushing gently with a plastic stick.



Each membrane can be individually moved up and down by a linear stepper motor (Nanotec L2818-OV15601 with SMCI36 controller), which is computer controlled using Labview. The membranes are fixed to the motor stick using two rubber plates (typical dimensions $ \unit{45}{\milli\meter}$, in black in Figure~\ref{fig:photos} right) which are sandwiched between two transparent and rigid PMMA plates. The motion of each motor is a sine function, and two neighboring membranes are driven by motors in opposition of phase so that the forcing is a standing wave whose azimuthal wavelength $\lambda_x \simeq \unit{200}{\milli\meter}$ is twice the curvilinear distance between $2$ motors. We note $M$ the amplitude of the motion of the membranes and we report here on measurements performed with a forcing period $T_f=\unit{15}{\second}$. The length of the stick of the motor number $m$ is:
\begin{equation}
h(t)=M \sin(\omega t+m \pi)+h_{0,m}
\end{equation}
where $\omega=2\pi/T_f$ and $h_{0,m}$ is the initial value of the stick length.

The fluid velocity is measured using particle image velocimetry (PIV). Since the velocities are slow, we use a low-cost laser sheet (Apinex, $\unit{5}{\milli\watt}$) and an industrial camera (Allied Vision Stingray) at $10$ or $12$ frames per second. The particles are imaged in the whole height on a width of about $\unit{30}{\milli\meter}$. Movies of $900$ images are acquired every
$\unit{100}{\second}$ or $1500$ images every $\unit{250}{\second}$. We made the seeding particles using wax and alumina powder, to obtain particles with a range of average densities similar to the range of densities in the fluid. The protocol to make the beads is detailed in Appendix~\ref{Fabrication of the seeding beads}. The beads are mixed with hand washing liquid and water ($50\%$) before being gently poured at the top of the fluid, to avoid air bubble accumulation at their surface.
We use OpenPiV~\cite{OpenPIV} on Matlab to determine the velocity field from two successive images. To take into account the small lens effect due to the cylinder, a different calibration length is used on the vertical and horizontal axis. The relaxation time of the particle is $t_b=\rho_b (2r_b)^2/(18 \eta)$, where $\rho_b$ is the density of the particle (slightly above that of water), $r_b \sim \unit{0.5}{\milli\meter}$ is the radius of the particle and $\eta$ the dynamic viscosity of water. An estimate of this time is thus $t_b \sim \unit{0.06}{\second} \ll T_f$: despite their large size, the particles are expected to follow the flow.
The typical interrogation windows has a width of $ 64 $ pixels and a height of $32 $ pixels (corresponding to $\unit{11}{\milli\meter}$), and the overlap between the interrogation windows is $50\%$. In each interrogation window, the
velocity as a function of time is fitted using a sine function to obtain the amplitude and the phase of the wave (oscillation at the forcing angular frequency $\omega$) and the mean flow. We have checked that the mean flow is azimuthal and only depends on the height $z$, except very close to the boundaries.


One of the main improvements between our set-up and the previous ones is the continuous injection of fluid at the bottom and removal at the top using a peristaltic pump (Ismatec Reglo ICC 4 channels). Even if the membrane motion at the top tends to slowly mix the fluid, this ensures that the stratification remains linear (see Figure~\ref{fig:stratification}). No measurement of mean flow lasting more than $2$ periods had been reported before our experiments. In the first experiments that we performed, without pump, we also observed only $1$ or $2$ oscillations of the mean flow, and the measured stratification at the top was very weak at the end of the experiment. Since the stratification affects the emission and propagation of waves in the system, the weakening of the stratification is consistent with the impossibility to generate a mean flow any more. To avoid such a weakening of the stratification due to mixing, salty fluid is continuously injected at the bottom (around $\unit{9 }{\milli\liter\cdot\reciprocal\minute}$) and removed at the top of the experiment in $3$ points: between the annulus of membranes and the cylinder (around $\unit{21 }{\milli\liter\cdot\reciprocal\minute}$) and on top of two membranes (around $\unit{14 }{\milli\liter\cdot\reciprocal\minute}$). This configuration ensures that no liquid accumulates at the top. The typical vertical velocity associated with this flow is $\unit{1}{\micro\meter\cdot\reciprocal\second}$, which is significantly smaller than the other velocities involved in the experiment. This modification allows us to perform quantitative measurements in a nearly steady state (variation of $N$ less than $\unit{0.1}{\radian \cdot \reciprocal \second}$ during a run of $\unit{90\, 000}{\second}$, see Figure~\ref{fig:stratification}). However, since the saturated concentration of the salt is finite, this system imposes a limit on the duration of the experiment. When the saturated fluid is injected, the experiment has to be stopped.


\begin{figure}[!htbp]
\centering
\includegraphics[width=0.5\textwidth]{fig3.pdf}
\caption{Stratification before (red dots) and after (blue diamonds) a $\unit{90 \, 000}{\second}$ long experiment. The orange dashed line is a linear fit of the red dots, and the cyan dotted line is a linear fit of the blue diamonds. Stratification using NaCl, $T= \unit{15}{\second}$, $M= \unit{13.5}{\milli\meter}$. The Brunt--V\"{a}is\"{a}l\"{a} frequency computed from the linear fits is $N=\unit{1.48}{\radian \cdot \reciprocal \second}$ before the experiment, and $N=\unit{1.40}{\radian \cdot \reciprocal \second}$ after.}
\label{fig:stratification}
\end{figure}



\section{Theoretical approach}\label{Theoretical approach}

\subsection{Description of the model}

The aim of this section is to present the model of~\cite{Plumb1978,Semin2016}, which is a minimal model of the QBO. In the next section, we present several analytical results obtained on this model. Some of them were mentioned in the supplementary material of~\cite{Semin2018}


In this model, a mean flow is generated by two counter-propagating internal gravity waves. These two waves propagate through a linearly stratified fluid, characterized by a Brunt--V\"{a}is\"{a}l\"{a} frequency, $N=\sqrt {-(g/\rho_0) \, \dd \rho_0 /\dd z}$, where $g$ is the acceleration of gravity, $\rho_0$ is the background density and $z$ is the upwards vertical axis. The angular frequency of the waves is $\omega$, one of the wave is propagating in the azimuthal direction $+x$, and the other in the direction $-x$ ; the horizontal wavevector is written $k_x$. These waves verify the linear dispersion relation for internal gravity waves:
\begin{equation}\label{eq:dispersionRelation}
\omega^2=N^2 (\sin(\theta))^2,
\end{equation}
where $\theta$ is the angle between the wavevector $\bm{k}$ and the vertical. Due to the structure of the internal gravity waves, the horizontal wave component $u'$ and the vertical component $v'$ verify:
\begin{equation}\label{eq:linkVelocity}
v'=u' \tan(\theta),
\end{equation}

The equation of the model are obtained using the Navier--Stokes equation, the Boussinesq approximation and the Wentzel--Kramers--Brillouin approximation. As in the original work of Plumb and McEwan~\cite{Plumb1978}, we consider upward propagating waves in the calculation. The following dimensionless equation is obtained for the mean flow $u$:
\begin{equation}\label{eq:dudtPlumb}
\dpart{u}{t}=- \dpart{D}{z}+\Lambda_1 \ddpart{u}{z}-\Lambda_2 u,
\end{equation}
where:
\begin{equation}\label{eq:fluxPlumb}
D=F \left(c_1 \exp \left[- \int_0^z \left(\frac{\alpha_1}{(u-c_1)^2} + \frac{\alpha_2}{(u-c_1)^4} \right) \ud z' \right] +c_2 \exp \left[- \int_0^z \left(\frac{\alpha_1}{(u-c_2)^2} + \frac{\alpha_2}{(u-c_2)^4} \right) \ud z' \right]
\right),
\end{equation}
traces back to the forcing of the mean flow by the waves. $\nu$ is the kinematic viscosity, $c_1=1$ and $c_2=-1$, $\Lambda_1$, $\Lambda_2$, $\alpha_1$, $\alpha_2$ are positive, $\alpha_1+\alpha_2=1$ and $F$ is the control parameter. An experiment is characterised entirely by $3$ dimensionless parameters, for instance $\Lambda_1$, $\Lambda_2$, and $\alpha_1$. However, we keep $F$ as a free parameter here because we will expand $F$ close to the threshold to apply the perturbation method. The boundary conditions are:
\begin{align}
u(0) &amp;= 0, \label{eq:boundarycondition1}
\\
u(+\infty) &amp; = 0,\label{eq:boundarycondition2}
\end{align}
and corresponds to a no slip boundary condition at $z=0$ like in~\cite{Plumb1978} and a vanishing velocity at infinity. We chose the latter condition to simplify the analysis and because the mean flow is actually localised close to the forcing boundary.

To make the link between the dimensionless variables and the physical variables, we note in the following set of equation and only in this set the dimensional variables with a $\sim$:
\begin{equation}\label{eq:defparam}
\begin{split}
\tilde{F} &amp;= \left&lt;\tilde{u'} \, \tilde{v'} \right&gt;, \\
\tilde{c} &amp;= \frac{\tilde{\omega}}{\tilde{k}_x}, \\
\tilde{d} &amp;= \left(\frac{\tilde{N} \tilde{\gamma}}{\tilde{k}_x \tilde{c}^2} + \frac{\tilde{N}^3 \tilde{\nu}}{\tilde{k}_x \tilde{c}^4} \right)^{-1}, \nonumber \\
z &amp;= \frac{\tilde{z}}{\tilde{d}}, \\
u &amp;= \frac{\tilde{u}}{\tilde{c}}, \\
\alpha_2 &amp;= \frac{\tilde{N}^3 \tilde{\nu} \tilde{d} }{\tilde{k}_x \tilde{c}^4},\\
\Lambda_1 &amp;= \frac{\tilde{\nu}\tilde{c}}{\tilde{F}(0) \tilde{d}},\\
\Lambda_2 &amp;= \frac{\tilde{\gamma}\tilde{c} \tilde{d}}{\tilde{F}(0) },\\
\end{split}
\end{equation}
where $\tilde{F}$ the momentum flux, $\tilde{c}$ the horizontal phase velocity, $\tilde{d}$ the vertical dissipation length, $\tilde{\nu}$ is the kinematic viscosity, $\tilde{\gamma}$ the dissipation due to the walls, $\tilde{F}(0)$ the momentum flux at $\tilde{z}=0$.


From the dimensionless solution of equation~\eqref{eq:dudtPlumb}, it is possible to obtain the dimensional variables: the dimensionless lengths should be multiplied by $\tilde{d}$ to obtain the dimensional lengths, the dimensionless times should be multiplied by $\tilde{c} \tilde{d}/\tilde{F}(0)$ and the dimensionless mean flow by $\tilde{c}$.



\subsection{Linear and weakly non-linear solutions of the equations}

For the dimensionless system of equations~\eqref{eq:dudtPlumb} and~\eqref{eq:fluxPlumb}, with the boundary conditions~\eqref{eq:boundarycondition1} and~\eqref{eq:boundarycondition2}, we have obtained an explicit solution of the linear problem from which an equation for the threshold is obtained. We have also derived the amplitude equation in the vicinity of the threshold. Details of the calculation are given in Appendix~\ref{Linear and weakly nonlinear solution} and we sum up the results here.


We introduce the stream function $\psi$
\begin{equation}\label{eq:defpsi}
\psi=\int_0^z u \ud z
\end{equation}
so that
\begin{equation}\label{eq:ufctpsi}
u=\dpart{\psi}{z}
\end{equation}

We are considering eigenmodes such that:
\begin{equation}
\dpart{\psi}{t}= \mu \psi
\end{equation}
with $\mu \in \bbC$.

We define the parameters
\begin{align}
a &amp;= \frac{F \times 2 (2\alpha_1+4\alpha_2) }{\Lambda_1}, \label{eq:defatext}
\\
b &amp;= \frac{\mu + \Lambda_2}{\Lambda_1}, \label{eq:defbtext}
\end{align}
that appear in the expression of $\psi$
\begin{equation}\label{eq:solgenetext}
\psi(z)=\psi_0 (z)\left(1+K f(z)\right).
\end{equation}
$\psi_0$, $K$ and $f$ are given by
\begin{align}
\psi_0 &amp;=J\left(2\sqrt{b}, 2 \sqrt{a} \exp(-z/2)\right),\label{eq:sollinhomogtext}
\\
K &amp;=\frac{- 1}{f(0)}= \frac{- 1}{ \int_{z_0}^0 \frac{\int_{+ \infty}^{z} \psi_0 (z') \ud z'}{(\psi_0(z))^2} \ud z}.\label{eq:valueofKtext}
\\
f(z) &amp;= \int_{z_0}^z \frac{\ud z' \int_{+ \infty}^{z'} \psi_0 (z'') \ud z''}{(\psi_0(z'))^2},\label{eq:valueftext}
\end{align}
where $J$ is the Bessel function of the first kind, of order $ 2\sqrt{b}$, and of argument $2 \sqrt{a} \exp(-z/2)$. $z_0$ is an arbitrary positive constant ; the final solution does not depend on the value chosen for $z_0$.



Using the boundary conditions and properties of the Bessel functions, we obtain an integral equation that relates the parameters of the problem and the growth rate
\begin{equation}\label{eq:linthresh4text}
\int_0^{1} \left[J\left(2\sqrt{b}-1, 2 \sqrt{a} v\right) + J\left(2\sqrt{b}+1, 2 \sqrt{a} v\right)\right] \ud v=0.
\end{equation}

We have solved equation~\eqref{eq:linthresh4text} numerically using a Matlab minimization routine and Bessel functions of complex order and argument~\cite{cbesselJ}. We fix the value of $\alpha_1$ (and thus $\alpha_2=1-\alpha_1$). We set $F=1$ which amounts to changing $\Lambda_i$ into $\tilde{\Lambda}_i=\Lambda_i/F$ and dropping the $\tilde{.}$ to ease notation. We then compute $a$ from equation~\eqref{eq:defatext}. Using equation~\eqref{eq:linthresh4text}, we obtain the corresponding value for $b$ at the threshold. Since $\mu$ is imaginary at the threshold and $\Lambda_2$ is real, we obtain from equation~\eqref{eq:defbtext} the value of $\Lambda_2=\Lambda_1 \Re(b) $ and $\mu= \Lambda_1 \Im(b) $. The eigenmodes are obtained using equation~\eqref{eq:solgenetext} with equations~\eqref{eq:valueftext}, \eqref{eq:valueofKtext} and~\eqref{eq:sollinhomogtext}, and their time evolution are given by multiplying the eigenmodes by $\exp(\mu t)$.



We have also derived the amplitude equation for the unstable mode in the vicinity of its onset of instability. We note
\begin{equation}
F= F_0 + \varepsilon F_1
\end{equation}
where $F_0$ is the value at the instability threshold ($\Lambda_1$ and $\Lambda_2$ being fixed) and $\varepsilon$ a small parameter (dimensionless distance to the threshold).


We introduce a slow timescale $T=\varepsilon t$ and write $u$ as:
\begin{equation}\label{eq:expansionofutext}
u=\sqrt{\varepsilon} \left(u_0 + \varepsilon u_1+ \varepsilon^2 u_2+\dots\right)
\end{equation}


At order $\varepsilon^{1/2}$, we recover the linear problem so that $u_0$ is the product of the eigenmode of the linear problem $u$ with a slowly varying amplitude:
\begin{equation}\label{eq:defmodeneutre}
u_0(t,T,z)=\frac{1}{2} \left(A(T) u(t,z) + A^*(T) u^*(t,z) \right),
\end{equation}
where $u$ is given by equation~\eqref{eq:ufctpsi} and~\eqref{eq:solgenetext}.


At order $\varepsilon^{3/2}$, using the solvability condition, we obtain the amplitude equation for $A$
\begin{equation}\label{eq:Lu1Hu0proj2text}
\dpart{A}{T}= \alpha A F_1 + \beta |A|^2 A,
\end{equation}
where explicit expressions for $\alpha$ and $\beta$ are given in the annex.



This is the amplitude equation of a pitchfork Hopf bifurcation. Indeed, this is constrained by the oscillatory nature of the unstable mode that breaks invariance under shift in the origin of time.


The coefficients can be calculated numerically using the expressions given in Appendix~\ref{Linear and weakly nonlinear solution}. In particular, we have calculated the following expression using a Matlab routine
\begin{equation}\label{eq:defStext}
S=\frac{\Re (\beta)}{\Re(\alpha)},
\end{equation}
where $\Re$ denotes the real part. The nature of the bifurcation is determined by the sign of $S$: if $S &gt;0$, the bifurcation is subcritical, and it is supercritical if $S &lt;0$.


The value of $S$ is determined by $a$ and $b$ through equation~\eqref{eq:sollinhomogtext} and the equations mentioned after. Since $\mu \in i \bbR$ at the threshold, the value of $b$ at the threshold only depends on the mean flow dissipation through $\Lambda_1$ and $\Lambda_2$ (see equation~\eqref{eq:defbtext}). The wave dissipation, i.e. the values of $\alpha_1$ and $\alpha_2$, only changes the value at the threshold of $F$. The nature of the bifurcation thus depends on the dissipation of the mean flow, whereas the nature of the dissipation of the waves only changes a prefactor which has always the same sign.



\subsection{Numerical simulations}

We have solved numerically the system of equations~\eqref{eq:dudtPlumb} and~\eqref{eq:fluxPlumb} using a centered second-order finite difference method with grid size $\delta z=10^{-2}$ (or $\delta z=4 \times 10^{-3}$ for the stability diagram). The total length is $4$ in dimensionless unit, i.e. $4d$ in physical units: this is sufficient to mimic an infinite domain, since the waves are almost entirely damped for $z=4$. Since the optimization of the numerical simulations is not the core of the present article, we have used an explicit Euler method in time with $\delta t=10^{-4}$ (or $\delta t=10^{-5}$ for the stability diagram), which is simpler than the schemes used by~\cite{Yoden1988} and~\cite{Couston2018}. The boundary conditions are no-slip at both the bottom and upper boundary, according to equations~\eqref{eq:boundarycondition1} and~\eqref{eq:boundarycondition2}. This numerical scheme is implemented in Matlab. The equations and the boundary conditions are not exactly the same as thee ones in the works of~\cite{Yoden1988} and~\cite{Couston2018}, who only consider bulk dissipation of the mean flow an a stress-free boundary condition at the upper boundary.

The numerical system can be solved using explicit Euler method, is numerically rather stable, and rapid to solve using a laptop (a few minutes with a $2$~GHz processor laptop, using a Matlab script). We point out that this system can be used as a non-trivial example of numerical resolution at the undergraduate level.



\section{Experimental results}\label{Experimental results}

\subsection{Interaction between wave and mean flow}



\begin{figure}
\centerline{\includegraphics[width=1\textwidth]{fig4T.pdf}}
\centerline{\includegraphics[width=1\textwidth]{fig4B.pdf}}
\caption{Measurements of the wave and mean flow for an experiment with $M=\unit{14.5}{\milli\meter}$ with MgCl$_2$ ($N=\unit{2.16}{\radian \cdot \reciprocal \second}$). The first row of the figure corresponds to the beginning of the forcing, when the mean flow is still small (sub-figures a,b,c). The second row corresponds to a time at which the mean flow is established. We chose a time at which the mean flow is negative everywhere for the simplicity of the interpretation (sub-figures d,e,f). The first column corresponds to the root-mean square value of the horizontal wave velocity $\left&lt;\sqrt{(u')^2}\right&gt;$ (sub-figures a,d). The second column corresponds to the wave phase, obtained from the horizontal wave velocity $u'$ (sub-figures b,e). The third column corresponds to the mean flow (sub-figures c,f).}
\label{fig:2017_07_25_Time_16_34_18_17}
\end{figure}


\begin{figure}
\centering
\includegraphics[width=0.4\textwidth]{fig5.pdf}
\vspace*{-10pt} 
\caption{Red dots: wave root-mean-square value averaged in the $x$ direction, from Figure~\ref{fig:2017_07_25_Time_16_34_18_17}$\MK$(a), i.e. $M=\unit{14.5}{\milli\meter}$ with MgCl$_2$, no mean flow. Dashed line: exponential fit, with $\mu_0=\unit{5.45}{\milli\meter \cdot \second^{-1}}$, $d_e=\unit{81}{\milli\meter}$. }\label{fig:2017_07_25_Time_16_34_rms1Donde}
\vspace*{-10pt} 
\end{figure}


The mechanism of the QBO is based on the interaction between the waves and the mean flow. We investigate this interaction in this section.

The characteristic of the waves and the mean flow $\unit{2}{\minute}$ after the beginning of the forcing at $M=\unit{14.5}{\milli\meter}$ is shown in Figures~\ref{fig:2017_07_25_Time_16_34_18_17}$\MK$(a,b,c). The waves are established, but the mean flow is still zero. This is in agreement with the study of a single wave in the same set-up~\cite{Semin2016}, and confirms that the characteristic time-scale of the mean flow is much larger than the one of the waves.

The initial phase of the wave is shown in Figures~\ref{fig:2017_07_25_Time_16_34_18_17}$\MK$(b). This phase corresponds to the expected stationary wave in the $x$ direction. The measured vertical wavelength $\lambda_z=\unit{46}{\milli\meter}$ is in reasonable agreement with the expected value $\lambda_z=\lambda_x \tan (\theta)\simeq \unit{40}{\milli\meter}$. The local amplitude of the wave is given by Figure~\ref{fig:2017_07_25_Time_16_34_18_17}$\MK$(a). The amplitude of the wave is not constant in the $x$ direction, as expected for the sum of two counter-propagative waves. The wave amplitude does not vanish, which is linked to the choice of the position of the field of view, and to the fact that it is shorter than half the wavelength in the $x$ direction $\lambda_x/2=\unit{100}{\milli\meter}$. By averaging the amplitude in the $x$ direction, we obtain Figure~\ref{fig:2017_07_25_Time_16_34_rms1Donde}. The amplitude of the wave decreases with height. This is due to the dissipation of the wave, mainly due to the dissipation by viscosity in the bulk: $\alpha_2=0.965$ where $\alpha_2$ is the dimensionless quantification of the bulk dissipation whereas $\alpha_1=0.035$ is the dissipation due to the wall. We can fit the curve with the function:
\begin{equation}\label{eq:fitampwaves}
u=\mu_0 \exp \left(-z/(2d_e) \right)
\end{equation}
where $d_e$ is the experimental decay length and $\mu_0$ is related to the amplitude of the wave at the top through $u_0=\mu_0/q$, where $q=1.6$ in our experiment takes into account that the measured amplitude is linked to the sum of two counter-propagating waves, measured in a finite field of view. The fit is good, indicating that the model is satisfactory. The value of $d_e$ from the fit is $d_e=\unit{81}{\milli\meter}$, is in rather good agreement with the theoretical value $d=\unit{95}{\milli\meter}$, obtained using $\gamma=\unit{10^{-3}}{\second}$ (see~\cite{Semin2016}), $\nu=\unit{10^{-6}}{\meter^2 \cdot \second^{-1}}$ and the measured value of $N$. The small discrepancy is probably due to the difference in viscosity, the solution of MgCl$_2$ being more viscous than water or NaCl solution. We did not perform a quantitative analysis because the concentration and thus the viscosity is not homogeneous in the solution.


We now consider the same forcing parameters but after $\unit{75}{\minute}$ in Figures~\ref{fig:2017_07_25_Time_16_34_18_17}$\MK$(d,e,f). The mean flow is well established at that time. We have chosen a snapshot at which the mean flow has the same sign everywhere so that the qualitative explanations are easier. The phase displayed in Figure~\ref{fig:2017_07_25_Time_16_34_18_17}$\MK$(d) is different from the one when there is no mean flow (Figure~\ref{fig:2017_07_25_Time_16_34_18_17}$\MK$(b): there is a kink in the phase at the top. At the bottom, the phase is tilted with an angle $\unit{10}{\degree}$, where the theoretical angle for a single wave is $\theta=\arcsin(\omega/N)=\unit{11}{\degree}$. This figure shows that there is a feedback of the mean flow on the wave, and that one of the wave is dominant at the bottom, the other being of much smaller amplitude. The wave at the bottom has a positive horizontal wave velocity (see for instance~\cite{Dauxois2017} for the basic geometry of internal gravity waves). This feedback is consistent with the denominator in equation~\eqref{eq:fluxPlumb}, and consistent with the qualitative explanation of the QBO~\cite{Baldwin2001} for which the wave which has a horizontal phase velocity of the same sign as the mean flow accelerates the mean flow more efficiently than the wave which has a horizontal phase velocity of the opposite sign, and is thus more attenuated. The modification of the waves by the mean flow has been shown in~\cite{Otobe1998} using the MoirÃ© technique and numerically in~\cite{Wedi2006}, but our method is more direct because we extract the wave phase. The feedback of the mean flow on the wave is also seen by comparing the amplitude of the wave without and with a mean flow in Figures~\ref{fig:2017_07_25_Time_16_34_18_17}$\MK$(a) and (d). In the latter case, the amplitude is lower everywhere but very close to the membrane. This is consistent with an increased attenuation of one of the waves due to the mean flow. The feedback of the mean flow on the waves for a single wave has been already shown experimentally in the geometry of a dipole and a wave beam in~\cite{Godoy-Diana_2006} and in a similar geometry in our previous work~\cite{Semin2016}. It has also been shown with two waves in~\cite{Otobe1998}, but with a more qualitative measurement of the waves.



\subsection{Characteristics of the mean flow}


\begin{figure}
\centerline{
\includegraphics[width=0.5\textwidth]{fig6.pdf} }
\caption{Example of mean flow as a function of time for a forcing amplitude $M=\unit{13.5}{\milli\meter}$, a period $T_f=\unit{15}{\second}$ and $ N=\unit{1.5}{\radian \cdot \reciprocal \second}$, stratification using NaCl.}\label{fig:spatios}
\end{figure}

\begin{figure}
\centerline{
\includegraphics[width=0.65\textwidth]{fig7.pdf} }
\caption{Mean flow as a function of time at $z=\unit{320}{\milli\meter}$ for the same experiment as the one shown in Figure~\ref{fig:spatios}: $M=\unit{13.5}{\milli\meter}$,$T_f=\unit{15}{\second}$, $ N=\unit{1.5}{\radian \cdot \reciprocal \second}$, NaCl. Blue points: experimental data, blue dashed line: curve linking the experimental points, black line: sine fit.}
\label{fig:2014_09_04_ecmoy_H320_b}
\end{figure}


An example of a space-time diagram of the mean flow, known as Hovm\"oller diagram in the atmospheric sciences, is shown in Figure~\ref{fig:spatios}. In this example, the mean flow is zero during a very long time and then suddenly appears without any apparent cause. This is consistent with a subcritical transition to the QBO for this value of $N$, as discussed in~\cite{Semin2018} and in Section~\ref{Stability diagram}. More than $9$ periods of the mean flow are observed in this example, which is much higher than in any previous experiments~\cite{Plumb1978,Otobe1998}. As mentioned in Section~\ref{Experimental set-up}, this has been achieved by the continuous injection of fluid at a very small velocity using an peristaltic pump.

The value of the mean flow at a given height for the same experiment is shown in Figure~\ref{fig:2014_09_04_ecmoy_H320_b}. This figure shows that the mean flow is actually quite regular, and can be fitted at a given height by a sine function during more than $9$ periods. This sinusoidal behaviour of the mean flow is consistent with the weakly nonlinear methods used in the theoretical approach of Section~\ref{Theoretical approach}.

\begin{figure}
\centerline{
\includegraphics[width=0.4\textwidth]{fig8L.pdf}
\includegraphics[width=0.45\textwidth]{fig8R.pdf} }
\caption{Left: relative phase of the mean flow, $M=\unit{12}{\milli\meter}$, MgCl$_2$. The horizontal dashed line corresponds to $\pi$. Right: amplitude of the mean flow as a function of height $z$. Blue dots: experiment, $M=\unit{12}{\milli\meter}$ (threshold at $M=\unit{11.3}{\milli\meter}$), MgCl$_2$; magenta line and diamond symbols: numerical resolution of the 1D model, $M=\unit{10.0}{\milli\meter}$ (threshold at $M=\unit{9.1}{\milli\meter}$), orange line: analytical linear profile, height coefficient chosen to match the data.}
\label{fig:2017_07_25_M13_MeanFlow}
\end{figure}



\begin{figure}
\centerline{
\includegraphics[width=0.45\textwidth]{fig9L.pdf}
\includegraphics[width=0.45\textwidth]{fig9R.pdf} }
\caption{Left: period of the mean flow as a function of the forcing amplitude $M$, for MgCl$_2$. Blue points: experiment, the forcing is changed every $\unit{8\, 000}{\second}$ to reach a quasi-steady state, and the period is measured in the last $\unit{4\, 000}{\second}$. See~\cite{Semin2018} for the space-time diagram. Magenta line and diamond symbols: numerical simulation. Right: amplitude as a function of $M$. Same color code. Height $z=\unit{320}{\milli\meter}$ for the experimental data, i.e. about $d_e$ from the forcing, and distance $d$ for the simulations. The blue dashed line and the violet dashed dotted line are of the form $\propto \sqrt{M-M_t}$, where $M_t$ is the threshold, which is different for the experiments and the simulations.}
\label{fig:2017_07_25_ecmoy_H_T_subexp}
\end{figure}


Since the mean flow is sinusoidal in time, it is possible to plot the relative phase of the mean flow as a function of height, as done in Figure~\ref{fig:2017_07_25_M13_MeanFlow}$\MK$left. The phase difference between the top and the bottom is close to $\pi$ in the field of view. This is consistent with the usual view of the QBO and the theoretical calculations where the mean flow at the top and at the bottom are in opposition of phase.

The amplitude of the mean flow as a function of height is displayed in Figure~\ref{fig:2017_07_25_M13_MeanFlow}$\MK$right. Experimental results, numerical and analytical results are shown in this figure. The chosen parameters are slightly different, due to the difference in the thresholds between the experiment and the numerical simulation. The value of $d$ for scaling the simulation is taken from Figure~\ref{fig:2017_07_25_Time_16_34_rms1Donde}, and the value of $F$ is extrapolated from the value of $\mu_0$ obtained from this figure, and by using the fact that $F$ is quadratic in $M$.
Since the wave velocity varies by about $20\%$ across the width~\cite{Semin2016}, and that the previous hypotheses may not be perfectly valid, we cannot expect a perfect quantitative agreement. However, Figure~\ref{fig:2017_07_25_M13_MeanFlow}$\MK$right shows that the values of the velocities and the shape of the mean flow profile are in qualitative agreement. In particular, the variation of the mean flow with height is not monotonic, leading to a local minimum in the amplitude. Experimentally, this effect is robust and has been observed irrespective of the stratification. Numerically, it is observed provided that the distance to the threshold is not too large. This shape is also observed in the linear analytical shape.

The period as a function of the forcing $M$ is shown in Figure~\ref{fig:2017_07_25_ecmoy_H_T_subexp}$\MK$left. The set of experiments used includes the one used in Figure~\ref{fig:2017_07_25_M13_MeanFlow}, and the numerical simulations are performed using the same parameters. The period is finite at the threshold and decreases when $M$ increases. The experiments and the simulations are in qualitative agreement. The period of the mean flow, about $T_{\QBO}=\unit{2\, 500}{\second}$ is much larger than the wave period $T_f=\unit{15}{\second}$ as expected. Note that obtaining a reasonable value of the mean flow period is not easy experimentally. It is achieved successfully here thanks to the continuous injection of new fluid, as shown by the failure of~\cite{Otobe1998} to obtain meaningful values for this period.

The mean flow amplitude as a function of the forcing $M$ is shown in Figure~\ref{fig:2017_07_25_ecmoy_H_T_subexp}$\MK$right. The value of the thresholds between experiment and theory are different but close. In both cases, the transition is continuous and the first point can be fitted by a form $\propto \sqrt{M-M_t}$. This behavior and the finite value of the period at the threshold shows that the bifurcation is a supercritical Hopf bifurcation in that experiment. The nature of the bifurcation will be investigated in more details in the next section.

\subsection{Stability diagram}\label{Stability diagram}

\begin{figure}
\centering
\includegraphics[width=0.55\textwidth]{fig10.pdf}
\vspace*{-10pt}
\caption{Stability diagram, from analytical solutions ($\alpha_2=1$), simulations, with experimental values. The solution is stable~(i.e. no mean flow) for large values of $\Lambda_1$ and $\Lambda_2$, and unstable (oscillating mean flow) for small values of $\Lambda_1$ and $\Lambda_2$. Red line and red circles: analytical values for the threshold, supercritical bifurcation; black circle: analytical value of the tricritical point; blue line and blue circles: analytical values for the linear threshold, subcritical bifurcation. Red squares: threshold from numerical simulations, case where the bifurcation is supercritical; blue squares: threshold from numerical simulations, case where the bifurcation is subcritical. Dashed orange line: line corresponding approximately to MgCl$_2$ during an experiment where the wave flux is changed; dashed violet line: same for NaCl. Plain orange diamond: experimental values where a mean flow is observed for MgCl$_2$; empty orange diamonds: same but without mean flow; orange asterisk: experimental values for which a mean flow is observed, with value of forcing interpolated. Violet diamonds: same as orange diamonds, but for NaCl.}
\label{fig:stabilitydiagram}
\end{figure}


The stability diagram is shown in Figure~\ref{fig:stabilitydiagram}. The bifurcation is subcritical at small values of $\Lambda_1$, and supercritical for large values of $\Lambda_1$. The tricritical point is at about $\Lambda_1=0.12$ and $\Lambda_2=1.88$.
The values of the thresholds are in agreement to within $3\%$ between the numerical and analytical solutions. The nature of the bifurcation is roughly in agreement.
The behaviour of the bifurcation for the two different salts depends on the ratio $\Lambda_2/\Lambda_1=\gamma d^2/\nu$. This ratio depends on $N$, and in our experiments of the salt because we use the more soluble and more expensive MgCl$_2$ only to achieve high values of $N$. The two lines displayed in Figure~\ref{fig:stabilitydiagram} correspond to two experiments, one with NaCl and a low value of $N$ and one with MgCl$_2$ and a larger value of $N$. During an experiment with MgCl$_2$ where the amplitude of the forcing is changed, the ratio $\Lambda_2/\Lambda_1$ changed slightly due to the change in the stratification (see Section~\ref{Experimental set-up}) and to the fact that the kinematic viscosity of a MgCl$_2$ solution depends on its density. The ratio $\gamma /\nu$ is independent of the viscosity, but $d$ decreases when the viscosity is increased. More precisely, it is possible to evaluate the value of $\Lambda_i$ for the point where there is no mean flow, because the wave amplitude can be fitted by an exponential (see an example in Figure~\ref{fig:2017_07_25_Time_16_34_rms1Donde}). This is the case for the largest forcing amplitude (plain diamonds), for which the mean flow is initially zero. It is also possible for the points at which there is no more mean flow (empty diamonds). In between, it is possible to interpolate the values of $\Lambda_i$ for the other values of the forcing for which there is a mean flow, assuming that $F(0)$ is quadratic in the forcing amplitude and that the viscosity increases linearly with time (the latter hypothesis is reasonable because the change is small). The dashed lines in Figure~\ref{fig:stabilitydiagram} correspond to a guide to the eye of where the experimental points should be without viscosity effects.

For MgCl$_2$, we have a rather good agreement between the experimental threshold and the theoretical one, without any adjustable parameter. This is not the case with NaCl. The reason is not clear, but we can note that the decay length $d$ for NaCl is different from the calculated one.


\section{Conclusion}

The quasi-biennial oscillation (QBO) is a the periodic reversal of the wind in the equatorial stratosphere, due to the forcing by atmospheric waves. The period, of about $28$ months, is not linked to any external astrophysical forcing, which is a remarkable feature~\cite{Hamilton1998}.

Three analog experiments have been performed in 1978 by Plumb and McEwan~\cite{Plumb1978}, in 1998 by Otobe et al.~\cite{Otobe1998} and in 2018 by Semin et al.~\cite{Semin2018}. The geometry is similar in all three experiments: a linearly stratified aqueous solution between two cylinders is forced using $16$ membranes. Simple $1$D model~\cite{Plumb1977} and a careful consideration of the dimensionless numbers of the problem have lead to the choice of this particular set-up.

The first experiment~\cite{Plumb1978} has provided a proof of the concept of the wave induced mean flow oscillation. The last experiment~\cite{Semin2018} has shown that the bifurcation between a state without mean flow and a state with a mean flow depends experimentally on the value of the Brunt--V\"{a}is\"{a}l\"{a} frequency. The model shows that the nature of the bifurcation depends on the main dissipation mechanism of the mean flow, i.e. in the bulk or in the boundary layer close to the wall of the cylinders. This main dissipation mechanism changes in our experiments when the Brunt--V\"{a}is\"{a}l\"{a} frequency changes.

In this article, we have presented the details of the experimental set-up discussed in~\cite{Semin2016,Semin2018}, including many experimental tricks.

We also give the details of the analytical resolution of the $1$D model of Plumb and McEwan~\cite{Plumb1978}. This model describe the generation of a mean flow by two counter-propagating internal gravity waves. We have presented an analytical solution of this model up to the third order, in order to determine the threshold, the eigenvectors and the nature of the bifurcation. Two dissipation mechanisms are taken into account for the waves and the mean flow: one linked to the effect of viscosity in the bulk, described by the dimensionless coefficient $\Lambda_1$ for the mean flow and $\alpha_2$ for the waves, and one for the dissipation at the lateral walls, described by the dimensionless coefficient $\Lambda_2$ for the mean flow and $\alpha_1$ for the waves. We have solved analytically this model, and we have obtained the value of the threshold as a function of the dimensionless parameters. The nature of the dissipation of the waves changes only quantitatively the threshold. In contrast, the nature of the dissipation of the mean flow plays a major role in the bifurcation: the bifurcation is supercritical for large values of $\Lambda_1$ and subcritical for small values of $\Lambda_1$.

It is remarkable that the model of Plumb and McEwan can be solved analytically. However, as shown in this article, the resolution is long. On the other hand, the numerical resolution of this model is simple and rapid: we think that it would be a good example for advanced undergraduate students.


We have measured simultaneously the waves and the mean flow using PIV. We have shown an example of feedback of the mean flow on the waves. The experimental results are consistent with the mechanism of the QBO: the attenuation of a wave (and its momentum transfer to the mean flow) depends on the relative sign of the mean flow and the horizontal phase velocity of the wave. In particular, at the special time where the mean flow has the same sign everywhere, one of the waves is filtered out and only the other wave is measured at the bottom.

We have also shown an example of a mean flow as a function of time for a fixed forcing. The example we have chosen correspond to a subcritical bifurcation: there is no mean flow during a very long time, and then it appears due to an unknown perturbation. The space-time diagram also shows that the oscillation of the mean flow is very stable, and last more that $9$ periods. This is significantly more than the duration of the oscillation of the mean flow reported in previous works~\cite{Plumb1978,Otobe1998}. This is due to the use of a pump to continuously inject fluid at the bottom and remove fluid at the top, which keep the stratification linear even close to the membranes and the injection of waves efficient.

The experimental mean flow at a given height can be well fitted by a sinusoidal fit. This is consistent with the weakly nonlinear approach used in the analytical model and is used to obtain the relative phase and the amplitude of the mean flow as a function of height $z$. The amplitude profile of the mean flow is qualitatively similar for the experimental, analytical and numerical results. In both cases, the amplitude is non-monotonous. We recover the main characteristics of the mean flow profile: a phase difference of $\pi$ between the top and the bottom, consistent with the usual qualitative description of the QBO~\cite{Baldwin2001}. The variation of the amplitude of the mean flow with height is not monotonic, displaying a local minimum in theory, simulations and experiments.



We have shown a set of values of the period of the mean flow as a function of the forcing. The experiments and the simulations are in good qualitative agreement and show that this period decreases when the amplitude of the forcing increases.

We also discuss in detail one of the main result of~\cite{Semin2018}: the bifurcation is always a Hopf one: when a mean flow is generated, it is always oscillating at a finite frequency, but the transition can be subcritical (sudden, possible hysteresis) or supercritical (continuous). The experimental results for the bifurcation are in agreement with the analytical results, and we have been able to find two different configurations for which the bifurcation is subcritical (small value of $N$, NaCl) or supercritical (large value of $N$, MgCl$_2$).


One of the perspective for the analog experiments is to study the effect of the ratio of the dissipation length to the height of water, and to force the membrane with two frequencies. The goal is to reproduce experimentally the numerical results of~\cite{Chartrand2024} about the robustness of the periodic behavior of the QBO.

\break
Another perspective for the analog experiments is to investigate a less regular forcing. It is surely feasible to use a similar set-up with membrane, and apply an irregular forcing, or at least use only a fraction of the motors. Obtaining a QBO without membrane is very difficult and has not been achieved yet. No oscillation has been observed in numerical simulations of the experiments for random forcing in the case of perfectly smooth membrane~\cite{Wedi2006}. Even with a very careful optimisation of the forcing and of the stratification it is not clear that such an experiment will work.

The geophysics perspective is the use of the tools of nonlinear physics, in particular bifurcation theory, for the QBO in GCM (general circulation model). Theses tools are becoming (again) very popular in atmospheric science to study the effects of climate change, and the bifurcations are one of the types of the so-called ``tipping points''~\cite{Bathiany2018,Wunderling2023,IPCC_2023}. The QBO is expected to decrease in amplitude due to climate change~\cite{Anstey2022}. The experiments have shown that the bifurcation can be subcritical or supercritical, depending mainly on the effective dissipation form of the mean flow. A first question in GCM will be to identify the dissipation of the main flow, which is neither molecular viscosity nor friction on the wall like in the experiment, but probably involves turbulence and maybe the coupling of the velocity field with other fields. An even more relevant question that can be addressed using GCM is the nature of the bifurcation: will the QBO disappear, and if yes continuously or suddenly? A first approach would be to plot the amplitude and period as a function of the CO$_2$ concentration or the wave forcing and search for the variations of the mean flow close to the onset of instability.




\appendix
\section{Fabrication of the seeding beads}\label{Fabrication of the seeding beads}

To make the seeding beads, we use the following ingredients:
\begin{itemize}
\item \emph{Paraffin wax:} We used paraffin wax with a solidification point at $\unit{52-54}{\celsius}$, supplier VWR, reference 26157.291

\item Aluminum oxide (alumina) $\unit{3}{\micro\meter}$ powder. Supplier Lam Plan reference A003 P105K

\item \emph{Hand washing liquid:} We use Green Lime Hand Washing Liquid 500 ml from the brand ``Arbre Vert''
\end{itemize}

It is very probably not necessary to use exactly the same ingredients to obtain beads, we just mention the ones that we use. We are not even sure that the hand washing liquid is necessary.


The protocol that we use is given below.

We mix $\unit{31}{\gram}$ of paraffin, $\unit{6}{\gram}$ of alumina and one droplet of washing-up liquid in a borosilicate glass beaker with a teflon stirring bar. Due to the size of alumina particles, we use a FFP2 mask for safety.

The beaker should be placed in a boiling water bath (``bain-marie'') and stirred. To do this, we put a crystallizing dish with water on a magnetic hotplate stirrer. The beaker with wax was placed inside the water, close enough to the surface of the hotplate so that the bar can stir the mixture.

When the wax is liquid, we poured the mixture in a plastic weighing boat to cool it down. It is possible to use a spoon to take the alumina which remains at the bottom of the beaker.

Then we let the alumina particles sediment. We started to stir again when the wax is very viscous but not yet solid.

When the wax is cold but still malleable by hand, we made a sphere with the mixture. We waited then one night so that the mixture is hard.

Then we used a chopper to make the particles (Moulinette Mini Chopper from the brand Moulinex).


It is important to stop after a few seconds (5-20) to avoid any melting of the wax.

The last step is to sieve the particles. The amount of beads around $\unit{1}{\milli\meter}$ is about $\unit{5-10}{\gram}$.



\section{Linear and weakly nonlinear solution}\label{Linear and weakly nonlinear solution}


Equations~\eqref{eq:dudtPlumb} and~\eqref{eq:fluxPlumb} are investigated analytically close to the threshold in the following.

We will determine the linear threshold, the eigenmodes and where the bifurcation is super-critical or sub-critical in the parameter space.

This will be achieved using a weakly non-linear expansion up to the third order in the amplitude of the mean flow $u$.

At this order, the operator $- \partial D / \partial z$ writes (see equation~\eqref{eq:fluxPlumb}) :
\begin{equation}\label{eq:dDdz3order}
- \dpart{D}{z}= F(L u+ \Phi (u))
\end{equation}
where $L$ is the linear operator:
\begin{equation}\label{eq:operatorL}
L u = - \dpart{}{z}\left(\exp(-z) \ell(u) \right)= \dpart{}{z}\left(2 (2\alpha_1+4\alpha_2) \exp(-z) \int_0^z u \ud z\right)
\end{equation}
$\ell (u)$ is defined by
\begin{equation}
\ell (u)= -2 (2\alpha_1+4\alpha_2) \int_0^z u \ud z'
\end{equation}
$\Phi$ is the non-linear operator.
\begin{equation}
\Phi(u)= - \dpart{}{z}\left(\exp(-z) \xi(u) \right)
\end{equation}
and $ \xi (u)$ is given by
\begin{multline}\label{eq:exprexiu}
\xi (u) = -2 (4\alpha_1+20\alpha_2) \int_0^z u^{3} \ud z' \\
+ 2(2\alpha_1+4\alpha_2) (3\alpha_1+10\alpha_2) \left(\int_0^z u \ud z' \right) \left(\int_0^z u^{2} \ud z' \right)
-\frac{1}{3}(2\alpha_1+4\alpha_2)^3 \left(\int_0^z u \ud z' \right)^3
\end{multline}


\subsection{Linear order}\label{Linear order}


At linear order, the model~\eqref{eq:dudtPlumb} reduces to:
\begin{equation}\label{eq:modelPlumblinear1}
\dpart{u}{t}=F \dpart{}{z}\left(2 (2\alpha_1+4\alpha_2) \exp(-z) \int_0^z u \ud z\right) +\Lambda_1 \ddpart{u}{z}-\Lambda_2 u
\end{equation}

We introduce the variable $\psi$ like in~\cite{Plumb1977}:
\begin{equation}\label{eq:defpsi}
\psi=\int_0^z u \ud z
\end{equation}

Then~\eqref{eq:modelPlumblinear1} becomes:
\begin{equation}\label{eq:modelPlumblinear2}
\frac{\partial^2 \psi}{\partial t \partial z}=F \dpart{}{z}\left(2 (2\alpha_1+4\alpha_2) \exp(-z) \psi \right) +\Lambda_1 \frac{\partial^3 \psi}{\partial z^3}-\Lambda_2 \dpart{\psi}{z}
\end{equation}


Using the linearity of the equations, we use complex variables, so that the solution is the real part of the corresponding complex variable. To simplify the notations, we use the same variable for the complex and real variables. We are considering eigenmodes of~\eqref{eq:modelPlumblinear2} so that:
\begin{equation}
\dpart{\psi}{t}= \mu \psi
\end{equation}
$\mu \in \bbC$. At the threshold, the real part of the growth rate vanishes and thus $\mu \in i \bbR$.

We then integrate~\eqref{eq:modelPlumblinear2}:
\begin{equation}\label{eq:modelPlumblinear3}
\mu \psi=F \times 2 (2\alpha_1+4\alpha_2) \exp(-z) \psi +\Lambda_1 \ddpart{\psi}{z} -\Lambda_2 \psi + K_1
\end{equation}
where $K_1$ is an unknown complex number. This equation can be rewritten:
\begin{equation}\label{eq:modelPlumblinear4}
\ddpart{\psi}{z}+ (a \exp(-z)-b) \psi= K,
\end{equation}
with:
\begin{align}
a &amp;= \frac{F \times 2 (2\alpha_1+4\alpha_2) }{\Lambda_1}, \label{eq:defa}
\\
b &amp;= \frac{\mu + \Lambda_2}{\Lambda_1}, \label{eq:defb}
\end{align}
and $K$ is unknown.

We first solve the homogeneous equation:
\begin{equation}\label{eq:modelPlumblinear5}
\ddpart{\psi}{z}+ (a \exp(-z)-b) \psi= 0.
\end{equation}
The solutions are:
\begin{equation}
\psi=J\left(\pm 2\sqrt{b}, 2 \sqrt{a} \exp(-z/2)\right),
\end{equation}
where $J$ is the Bessel function of the first kind, of order $ \pm 2\sqrt{b}$ and of argument $2 \sqrt{a} \exp(-z/2)$.


The function $J(- 2\sqrt{b}, 2 \sqrt{a} \exp(-z/2))$ does not have a finite limit when $z\rightarrow + \infty$, and we disregard this solution due to the boundary condition (see equation~\eqref{eq:boundarycondition2}).

The solution of~\eqref{eq:modelPlumblinear5} is thus:
\begin{equation}\label{eq:sollinhomog}
\psi_0=J\left(2\sqrt{b}, 2 \sqrt{a} \exp(-z/2)\right),
\end{equation}
up to a multiplicative constant, which is irrelevant for the linear problem.


We then solve~\eqref{eq:modelPlumblinear4}, and we search a particular solution using the method of the variation of parameter:
\begin{equation}
\psi(z)=f_1(z) \psi_0 (z).
\end{equation}


Using equation~\eqref{eq:modelPlumblinear4} and that $\psi_0$ is solution of the homogeneous equation, we obtain:
\begin{equation}\label{eq:modelPlumblinear7}
2 \deriv{f_1}{z} \deriv{\psi_0}{z}+ f_1 \dderiv{\psi_0}{z} = K.
\end{equation}
We multiply the previous equation by $\psi_0$:
\begin{equation}
\deriv{}{z} \left(\psi_0^2 \deriv{f_1}{z} \right)= K \psi_0,
\end{equation}
which gives:
\begin{equation}
\psi_0^2 \deriv{f_1}{z} = \int_{+ \infty}^z K \psi_0 (z') \ud z',
\end{equation}
assuming that $f_1$, $\deriv{f_1}{z}$ vanish when $z \rightarrow + \infty$ and using that $\psi_0$ vanishes when $z \rightarrow + \infty$.

We write $f_1=K f$, so that we have:
\begin{equation}\label{eq:valuef}
f(z)= \int_{z_0}^z \frac{\ud z' \int_{+ \infty}^{z'} \psi_0 (z'') \ud z''}{\left(\psi_0(z')\right)^2},
\end{equation}
where $z_0$ is an arbitrary positive constant ; the final solution does not depend on the value chosen for $z_0$.

So the solution of~\eqref{eq:modelPlumblinear4} is thus:
\begin{equation}\label{eq:solgene}
\psi(z)=\psi_0 (z)\left(1+K f(z)\right),
\end{equation}
where $\psi_0$ is given by~\eqref{eq:sollinhomog} and $f$ by~\eqref{eq:valuef}.


We can also write the solution for $u$:
\begin{equation}\label{eq:ufctionpsi0}
u=\dpart{}{z} \left(\psi_0 (z)\left(1+K f(z)\right) \right),
\end{equation}

Using~\eqref{eq:boundarycondition2}, \eqref{eq:defpsi}, we have the following boundary conditions:
\begin{align}
&amp; \psi (0)=0, \label{eq:BCpsi0}
\\
&amp; \deriv{\psi}{z}(0) = 0, \label{eq:BCdpsi0}
\\
&amp; \deriv{\psi}{z}(+ \infty)= 0, \label{eq:BCdpsiinfty}
\end{align}



Using~\eqref{eq:BCpsi0} in~\eqref{eq:solgene} gives:
\begin{equation}
0=\psi_0 (0)(1+K f(0)).
\end{equation}
So that:
\begin{equation}\label{eq:valueofK}
K=\frac{- 1}{f(0)}= \frac{- 1}{ \int_{z_0}^0 \frac{\int_{+ \infty}^{z} \psi_0 (z') \ud z'}{(\psi_0(z))^2} \ud z}.
\end{equation}

Using~\eqref{eq:BCdpsi0} in~\eqref{eq:solgene} gives:
\begin{equation}
0 = \psi_0'(0)\left(1+K f(0)\right)+\psi_0 (0)\left(K f'(0)\right)= \psi_0 (0)\left(K f'(0)\right).
\end{equation}
We thus have $f'(0)=0$, i.e.:
\begin{equation}
\int_0^{+ \infty} \psi_0(z) \ud z=0
\end{equation}
i.e.:
\begin{equation}\label{eq:linthresh2}
\int_0^{+ \infty} J\left(2\sqrt{b}, 2 \sqrt{a} \exp(-z/2)\right) \ud z=0
\end{equation}
This equation allows us to compute the threshold of the instability.

Using the new variable
\begin{equation}
v=\exp(-z/2),
\end{equation}
equation~\eqref{eq:linthresh2} gives:
\begin{equation}\label{eq:linthresh3}
\int_0^{1} \frac{J\left(2\sqrt{b}, 2 \sqrt{a} v\right)}{v} \ud v=0.
\end{equation}


Using the following property of Bessel functions:
\begin{equation}
\frac{2 \alpha}{y} J\left(\alpha,y\right)=J\left(\alpha-1,y\right)+J\left(\alpha+1,y\right),
\end{equation}
we finally obtain equation~\eqref{eq:linthresh4text}.



\subsection{Non-linear order}


Using equation~\eqref{eq:dDdz3order}, equation~\eqref{eq:dudtPlumb} is considered up to the third order in the amplitude of $u$:
\begin{equation}\label{eq:dudtPlumborder3}
\dpart{u}{t}=F(L u+ \Phi (u))+\Lambda_1 \ddpart{u}{z}-\Lambda_2 u
\end{equation}

More precisely, we note
\begin{equation}
F= F_0 + \varepsilon F_1
\end{equation}
where $F_0$ is the value at the instability threshold ($\Lambda_1$ and $\Lambda_2$ being fixed) and $\varepsilon$ a small parameter (dimensionless distance to the threshold).

We introduce a slow timescale $T=\varepsilon t$ which correspond to the slow variation of the amplitude of the oscillation. Using the multiple scale expansion, $\dpart{u}{t}$ should be replaced by :
\begin{equation}\label{eq:dudtexpansion}
\dpart{u}{t}+\dpart{u}{T} \dpart{T}{t}=\dpart{u}{t}+ \varepsilon \dpart{u}{T}
\end{equation}

We choose the expansion for $u$:
\begin{equation}\label{eq:expansionofu}
u=\sqrt{\varepsilon} \left(u_0 + \varepsilon u_1+ \varepsilon^2 u_2+\dots\right)
\end{equation}



At order $\varepsilon^{1/2}$, equation~\eqref{eq:dudtPlumborder3} gives:
\begin{equation}
\dpart{u_0}{t}= F_0 L u_0 +\Lambda_1 \ddpart{u_0}{z}-\Lambda_2 u_0
\end{equation}
This is the linear order (see equation~\eqref{eq:modelPlumblinear1}), which has already been solved. The solution is proportional to
\begin{equation}\label{eq:modeneutre}
u=\dpart{\psi}{z}
\end{equation}
where $\psi$ is given by equations~\eqref{eq:solgene}, \eqref{eq:valueofK}, \eqref{eq:valuef}, and~\eqref{eq:sollinhomog}.

We write $u_0$ as a product of the eigenmode of the linear problem $u$ with a slowly varying amplitude:
\begin{equation}\label{eq:defmodeneutre}
u_0(t,T,z)=\frac{1}{2} \left(A(T) u(t,z) + A^*(T) u^*(t,z) \right),
\end{equation}
where $u$ is given by equation~\eqref{eq:modeneutre}

At order $\varepsilon^{3/2}$, equation~\eqref{eq:dudtPlumborder3} gives:
\begin{equation}
\dpart{u_1}{t} +\dpart{u_0}{T} = F_0 L u_1+ F_1 L u_0+ \Phi(u_0) +\Lambda_1 \ddpart{u_1}{z}-\Lambda_2 u_1
\end{equation}
which can be rewritten:
\begin{equation}\label{eq:Lu1Hu0}
\underbrace{\dpart{u_1}{t}- F_0 L u_1- \Lambda_1 \ddpart{u_1}{z}+\Lambda_2 u_1}_{\clL u_1}=\underbrace{-\dpart{u_0}{T}+ F_1 L u_0+ \Phi(u_0)}_{\clH u_0}
\end{equation}


The next step is to calculate the adjoint operator of $\clL $. We introduce the (hermitian) scalar product:
\begin{equation}
\langle a|b \rangle=\lim_{T\,\rightarrow\,+ \infty} \frac{2}{T} \int_{-T}^{T} \int_0^{+\infty} a^* b \ud t \ud z'
\end{equation}
where $a^*$ is the complex conjugate of $a$.

The fields verify the following properties: they do not diverge when $t \rightarrow \pm \infty$ and they vanish when $x=0$ or $+ \infty$.

After integration by parts we obtain for the adjoint operator:
\begin{equation}
\clL ^\dagger c=+\dpart{c}{t}+ F_0 \times 2 \left(2 \alpha_1+4\alpha_2\right)\int_0^{z} \exp(-z') \dpart{c}{z} \ud z' + \Lambda_1 \ddpart{c}{z} -\Lambda_2 c
\end{equation}

We now search a vector $q$ in the kernel of $ \clL ^\dagger$:
\begin{equation}
\clL ^\dagger q=0
\end{equation}
We thus have:
\begin{equation}
\dpart{}{z} \left(\clL ^\dagger q \right) = 0
\end{equation}

We define $\varphi$ by:
\begin{equation}
\varphi=\dpart{q}{z}
\end{equation}

We search modes proportional to $\exp{(\kappa t)}$ so that
\begin{equation}
\dpart{\varphi}{t}= \kappa \varphi
\end{equation}
We thus obtain:
\begin{equation}
\ddpart{\varphi}{z} + \left(a \exp(-z)-B\right) \varphi = 0
\end{equation}
with:
\begin{align}
\begin{split}
a &amp;= \frac{F_0 \times 2 \left(2 \alpha_1+4\alpha_2\right)}{\Lambda_1}\\
B &amp;= \frac{-\kappa +\Lambda_2}{\Lambda_1}
\end{split}
\end{align}
This equation is nearly identical to the one solved for the direct problem (equation~\eqref{eq:modelPlumblinear4}). We thus have:
\begin{equation}
\varphi = J \left(2 \sqrt{B}, 2 \sqrt{a} \exp(-z/2)\right)
\end{equation}
so that we have finally obtained an expression of an element of the kernel of the adjoint of the linear operator:
\begin{equation}
q=\int_0^z J \left(2 \sqrt{B}, 2 \sqrt{a} \exp\left(-z'/2\right)\right) \ud z'
\end{equation}

From equation~\eqref{eq:Lu1Hu0}, the solvability condition then writes:
\begin{equation}\label{eq:Lu1Hu0proj}
\left\langle q \,\middle|\, \clL  u_1 \right\rangle = \left\langle q\,\middle|\, -\dpart{u_0}{T}+ F_1 L u_0+ \Phi(u_0)\right\rangle =0
\end{equation}

Using the expression of $u_0$, see eq. \eqref{eq:defmodeneutre}, we obtain the amplitude equation:
\begin{equation}\label{eq:Lu1Hu0proj2}
\dpart{A}{T}= A F_1 \frac{\left\langle q\,\middle|\, L u\right\rangle}{\left\langle q\,\middle|\,u\right\rangle}+ |A|^2 A \frac{\left\langle q\,\middle|\,\Phi\left(\frac{1}{2} \left(u(t,z) + u^*(t,z)\right)\right)\right\rangle}{\left\langle q\,\middle|\,u\right\rangle},
\end{equation}
from which we obtain the expression of the coefficients $\alpha$ and $\beta$ as defined in Eq.~\eqref{eq:Lu1Hu0proj2text}:
\begin{align}
\alpha&amp;=\frac{\left\langle q\,\middle|\,L u\right\rangle}{\left\langle q\,\middle|\,u\right\rangle}\,,\\
\beta&amp;=\frac{\left\langle q\,\middle|\,\Phi\left(\frac{1}{2} \left(u(t,z) + u^*(t,z)\right)\right)\right\rangle}{\left\langle q\,\middle|\,u\right\rangle}\,.
\end{align}

The different coefficients are given by
\begin{equation}
\begin{split}
\left\langle q\,\middle|\,u\right\rangle &amp;= \int_0^{+\infty} \left(\int_0^{z} J \left(2 \sqrt{B}, 2 \sqrt{a} \exp\left(-z'/2\right)\right) \ud z' \right)^* \left(\frac{1}{2} \dpart{\psi}{z} \right) \ud z, \\
&amp;= - \frac{1}{2} \int_0^{+\infty} J \left(2 \sqrt{B}, 2 \sqrt{a} \exp\left(-z'/2\right)\right) ^* \psi \ud z.
\end{split}
\end{equation}
where $\psi$ is given by eq~\eqref{eq:solgene}, and
\begin{align}
\begin{split}
\frac{1}{2(2 \alpha_1+4 \alpha_2)} \left\langle q\,\middle|\,L u\right\rangle &amp;= \int_0^{+\infty} \left(\int_0^{z} J \left(2 \sqrt{B}, 2 \sqrt{a} \exp\left(-z'/2\right)\right) \ud z' \right)^* \left(\frac{1}{2} \dpart{}{z}\left(\exp(-z) \psi(z) \right) \right) \ud z\\
&amp;= - \frac{1}{2} \int_0^{+\infty} J \left(2 \sqrt{B}, 2 \sqrt{a} \exp\left(-z'/2\right)\right) ^* \exp(-z) \psi \ud z
\end{split}
\end{align}
and for the non-linear term by
\begin{equation}
\left\langle q\,\middle|\,\Phi\left(\frac{1}{2} \left(u(t,z) + u^*(t,z)\right)\right)\right\rangle=\left\langle q\,\middle|\,\dpart{}{}\left(\exp(-z) \xi(u)\right) \right\rangle = Q_1+ Q_2+ Q_3,
\end{equation}
with:
\begin{align*}
Q_1=2 (4\alpha_1+20\alpha_2)\int_0^{+\infty} \ud z J \left(2 \sqrt{B}, 2 \sqrt{a} \exp(-z/2)\right)^* \exp(-z)\int_0^z \frac{3}{8} \left|u(z')\right|^2 u(z') \ud z', 
\end{align*}
\begin{multline*}
Q_2=- 2(2\alpha_1+4\alpha_2) (3\alpha_1+10\alpha_2) \int_0^{+\infty} \ud z  J \left(2 \sqrt{B}, 2 \sqrt{a} \exp(-z/2)\right)^* \exp(-z) \\
\frac{1}{8} \left(\left(\int_0^z u(z')^* \ud z'\right) \left(\int_0^z u(z')^2 \ud z'\right) +2 \left(\int_0^z u(z') \ud z'\right) \left(\int_0^z |u(z')|^2 \ud z'\right)  \right), 
\end{multline*}
\begin{align*}
Q_3= \frac{1}{3}(2\alpha_1+4\alpha_2)^3 \int_0^{+\infty} \ud z J \left(2 \sqrt{B}, 2 \sqrt{a} \exp(-z/2)\right)^* \exp(-z) \frac{3}{8} |\psi(z)|^2 \psi(z). \\
\end{align*}

\section*{Acknowledgements}

We gratefully acknowledge J. da Silva Quintas, C. Goncalves, E. Nicolau, C. Herrmann, and L.~Bonnet for their technical support. We thank our collaborators S. Fauve, N. Garroum and G.~Facchini for discussions and with whom~\cite{Semin2018} or~\cite{Semin2016} were written.

\section*{Declaration of interests}
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

%\bibliographystyle{crplain}
%\bibliography{Semin}

\printbibliography
\end{document}
</pre></body></html>