# Recovering Quasi-Biennial Oscillations from Chaos

###### Abstract

The Quasi-Biennial Oscillation (QBO) is understood to result from wave-mean-
flow interactions, but the reasons for its relative stability remain a subject
of ongoing debate. In addition, consensus has yet to be reached regarding the
respective roles of different equatorial wave types in shaping the QBO’s
characteristics. Here we employ Holton-Lindzen-Plumb’s quasilinear model to
shed light on the robustness of periodic behavior in the presence of multiple
wave forcings. A comprehensive examination of the various dynamical regimes in
this model reveals that increased vertical wave propagation at higher
altitudes favors periodicity. In the case of single standing wave forcing,
enhanced vertical propagation is controlled by the wave attenuation length
scale. The occurrence of non-periodic states at high forcing amplitudes is
explained by the excitation of high vertical unstable modes. Increasing the
attenuation length scale prevents the emergence of such modes. When multiple
wave forcing is considered, the mean flow generated by a dominant primary wave
facilitates greater vertical propagation of a perturbation wave. Raising the
altitude where most of the wave damping occurs favors periodicity by
preventing the development of secondary jets responsible for the aperiodic
behavior. This mechanism underscores the potential role of internal gravity
waves in supporting the periodicity of a QBO primarily driven by planetary
waves.

##  1 Introduction

The Quasi-Biennial Oscillation (QBO) of equatorial winds, once considered a
hallmark of low-frequency periodic phenomena emerging from the complex
interplay of waves and turbulence (Baldwin et al., 2001), has exhibited a
surprising departure from its regularity in recent years. Notably,
interruptions and the emergence of a new vertical mode in 2016 and 2020, have
raised questions about the long-standing assumption of its steady periodicity
(e.g., Newman et al., 2016; Osprey et al., 2016). These disruptions, coupled
with similar observations of irregular wind reversals on other planets
(Fletcher et al., 2017), have cast uncertainty on the robustness of this
equatorial oscillation on Earth (Anstey et al., 2022).

The possibility of non-periodic states in equatorial wind reversals,
particularly when wave forcing is intensified, has been documented in various
idealized studies. These include simulations based on a 1D quasilinear model
(Kim and MacGregor, 2001; Renaud et al., 2019; Léard et al., 2020), 2D direct
numerical simulations of stratified fluids subject to lower boundary forcing
(Wedi and Smolarkiewicz, 2006; Couston et al., 2018; Renaud et al., 2019), and
3D general circulation simulations of gas giants and brown dwarfs (Showman et
al., 2019). These findings highlight the absence of a definitive explanation
for the remarkable periodicity of the QBO on Earth, prompting the question of
why the atmosphere appears to be finely tuned to parameters that correspond to
periodic solutions.

In this study, building upon previous work by Saravanan (1990) and Léard et
al. (2020), we delve into the role of multiple wave forcing and critical
layers in promoting periodicity. We use the quasilinear model originally
developed by Holton, Lindzen, and Plumb (Holton and Lindzen, 1972; Plumb,
1977). This model serves as an invaluable tool for exploring the dynamics of
equatorial wind reversals, replicating fundamental aspects of wave–mean-flow
interactions, while allowing the exploration of wide regions of the parameter
space.

Within the framework of the quasilinear model, several mechanisms have been
identified that can support the restoration of periodic states. A first
mechanism is due to tropical upwelling, represented as a constant-in-time,
bottom-intensified vertical transport term in zonally averaged zonal momentum
equations (Saravanan, 1990). This term effectively shifts the range of
observable periodic states to higher wave forcing values in the parameter
space (Rajendran et al., 2016). Upwelling’s significance in the altitude range
connecting the upper troposphere and the lower stratosphere in Earth’s
atmosphere is now widely acknowledged (Match and Fueglistaler, 2019, 2021).

A second mechanism is due to the introduction of a forcing term with
prescribed low-frequency oscillations such as a seasonal cycle (Read and
Castrejón-Pita, 2012). Such periodic forcing can result from various factors,
including seasonal modulations of upwelling (Rajendran et al., 2016), coupling
with semi-annual oscillations in the upper stratosphere (Dunkerton and Delisi,
1997), or seasonal variations of extratropical Rossby wave momentum flux
(Bardet et al., 2022).

A third mechanism was introduced by Léard et al. (2020), suggesting
periodicity recovery through the redistribution of forcing associated with a
monochromatic wave over a broad spectrum of frequencies via multimodal forcing
at the stratosphere’s base. They hypothesize that high-frequency waves, which
tend to yield more periodic oscillations for a given forcing, become dominant
in governing mean-flow oscillations beyond a specific threshold. However,
since this forcing redistribution simultaneously reduces the amplitude of the
primary wave while increasing the amplitude of high-frequency waves, the role
of waves other than the high-frequency ones remains challenging to decipher.

Here, we uncover two additional mechanisms that favor periodicity recovery in
the absence of a seasonal cycle and tropical upwelling. Firstly, we emphasize
the critical role played by a geometrical parameter that compares the depth of
the stratosphere to the typical wave attenuation length scale. We demonstrate
that an increase in the attenuation length scale promotes the resurgence of
periodic regimes in cases with monochromatic wave forcing. Secondly, we
explore the role of multimodal forcing by considering a scenario with two
pairs of waves: a fixed background standing wave and a perturbation standing
wave, each having distinct wavenumbers and frequencies. Results show
substantial regions of parameter space wherein the perturbation drives the
system to return to a periodic state, exhibiting a wind structure reminiscent
of the Quasi-Biennial Oscillation. Periodicity recovery occurs whether or not
waves with higher frequencies govern zonal wind oscillations. Our analysis
points to a synchronization mechanism, which coordinates the descent rate of
critical layers of both types of waves, preventing the growth of higher
vertical modes and explaining the periodicity recovery.

In Section 2, we introduce the quasilinear model. Section 3 presents the model
results, followed by a concluding discussion in Section 4.

##  2 Model

The quasilinear model provides a description of the coupled interaction
between internal gravity waves and zonal winds in a stratified fluid (Plumb,
1977; Vallis, 2017). In this model, the stratosphere is represented as a
horizontally averaged one-dimensional vertical profile of linearly stratified
fluid with a constant buoyancy frequency, denoted as $N$. The model tracks the
evolution of zonally averaged zonal winds, $\overline{u}$, under the influence
of waves or deviations to the zonal mean, represented as $u^{\prime}$. Forcing
from a wave generator at the lower boundary mimics the impact of the
tropopause on stratospheric dynamics. The fluid is confined within a layer of
depth $H$, and the mean flow is subject to stress-free and no-flow boundary
conditions at the top and bottom, respectively. Additionally, a radiation
boundary is applied to the waves at the top boundary.

In this study, our primary focus is on a scenario in which the wave generator
emits a maximum of two pairs of symmetric and counterpropagating waves. The
first pair, indexed as $b$, is specifically designed to establish the
background state of the QBO. Our objective is to investigate how this
background state is perturbed by the introduction of another pair, referred to
as $p$. The model describes the coupled interaction between the mean flow and
the waves as

$\displaystyle\partial_{t}\overline{u}$
$\displaystyle=\nu\partial_{zz}\overline{u}-\sum_{i}\partial_{z}\overline{u^{\prime}w^{\prime}}_{i}\,,\quad
i\in\\{b,p\\}\,,$ (1) $\displaystyle\overline{u^{\prime}w^{\prime}}_{i}$
$\displaystyle=F_{i}\exp\left(-\frac{1}{h_{i}}\int^{z}_{0}\frac{c_{i}^{2}}{(\overline{u}-c_{i})^{2}}\,\mathrm{d}z\right)-F_{i}\exp\left(-\frac{1}{h_{i}}\int^{z}_{0}\frac{c_{i}^{2}}{(\overline{u}+c_{i})^{2}}\,\mathrm{d}z\right)\,,$
(2)

where Eq. (1) describes the time evolution of the zonal winds forced by the
waves’ Reynolds stresses, $\overline{u^{\prime}w^{\prime}}_{i}$, and Eq. (2)
describes how these Reynolds stresses are in turn modulated by the zonal
winds. Within this framework, the parameters include $\nu$ representing the
eddy viscosity, $\pm F_{i}$ denoting the imposed wave momentum flux for the
westward $(-)$ and eastward $(+)$ propagating waves, $\pm c_{i}$ representing
the waves’ horizontal phase speeds, and $h_{i}$ describing the typical
vertical attenuation length scale of the wave field. The model is complemented
by two boundary conditions111Our parameterization of wave momentum flux
assumes implicitly a radiation condition for waves in the upper stratosphere,
such that wave momentum leaks out of the domain. Changing these boundary
conditions would change the expression of $\overline{u^{\prime}w^{\prime}}$
and details of the phase diagram in the case $\widehat{h}>1$, but not the
mechanisms described in this study. for the zonal winds:

$\overline{u}=0\,\,\,\text{at}\,\,\,z=0\,,\quad\partial_{z}\overline{u}=0\,\,\,\text{at}\,\,\,z=H\,.$
(3)

This model can be derived from the nonrotating primitive equations describing
the motion of a stratified fluid on the equatorial plane, combined with a
quasilinear approach. The derivation consists of three main steps. First, the
dynamics are linearized around a prescribed mean flow and the wave field is
computed assuming that the mean flow is frozen-in-time. Second, the Reynolds
stresses induced by this wave field are computed in the large-time limit,
assuming a scale separation between the vertical oscillations of the wave
field and the vertical variations of the winds. The final step involves
parameterizing the Reynolds stresses using a WKB (Wentzel-Kramers-Brillouin)
approximation to evolve the mean flow. Detailed discussions of the underlying
assumptions leading to this model can be found in Renaud and Venaille (2020).

Equation (2) for the Reynolds stresses is based on the assumption that wave
field dissipation is primarily caused by radiative damping, which occurs at a
constant rate denoted as $\gamma$. This parameter determines the attenuation
length scale, denoted as $h_{i}$, for a wave with a given frequency
$\omega_{i}$ and wavenumber $k_{i}$ through the relationship:

$h_{i}=\frac{\omega_{i}^{2}}{\gamma Nk_{i}}\,.$ (4)

For a given wave, fixing the parameter couple $(\omega_{i},k_{i})$ is
equivalent to fixing $(c_{i},h_{i})$, as the phase speed is given by
$c_{i}=\omega_{i}/k_{i}$. This relationship holds under the condition that the
buoyancy frequency and damping rate remain steady in time and consistent
across the vertical axis. The attenuation length scale, in combination with
the forcing amplitude $F_{i}$ and the phase speed $c_{i}$, provides an
estimate for the characteristic timescale of mean flow reversals for a single
monochromatic forcing, which can be expressed as $\tau_{i}=c_{i}h_{i}/F_{i}$
(Vallis, 2017).

We emphasize that the parameter $h_{i}$ characterizes the typical exponential
decay scale of the wave forcing profile in the absence of a mean flow.
However, in the presence of a mean flow, the local attenuation length scale at
a given altitude is determined by replacing $\omega_{i}$ with the Doppler-
shifted frequency $\omega_{i}-k_{i}\overline{u}(z)$ in Equation (4). This
modification accounts for the influence of the mean flow on the attenuation of
the wave field, which in turns modifies vertical profile of the wave forcing.
In the case of a constant mean flow, a shallower decay scale occurs when the
phase speed $c_{i}$ shares the same sign as the mean flow $\overline{u}$, and
conversely when $c_{i}$ has the opposite sign to $\overline{u}$. In general, a
sheared mean flow leads to significant deviations from an exponential wave
forcing profile.

A crucial effect of this Doppler-shift arises when the magnitude of the mean
flow $\overline{u}$ becomes comparable to $c_{i}$. This leads to a situation
where the local attenuation length scale approaches zero, effectively causing
complete absorption of wave momentum into the mean flow at this particular
altitude, thus preventing any wave forcing above it. This phenomenon is
referred to as a critical layer. The altitude at which these critical layers
can occur has significant implications for the dynamical behavior of the
oscillation.

In this study, we explore two distinct approaches to raise this altitude: (i)
in the case of monochromatic forcing, by increasing the attenuation length
scale $h_{i}$, and (ii) in the case of wave superposition, by elevating the
mean flow $\overline{u}$ to a value greater than the phase speed $c_{i}$. In
both cases, increasing the altitude at which critical layers can occur
promotes periodicity.

Note that the critical layer absorption phenomenon takes place irrespective of
the magnitude of the damping rate and the specific nature of the dissipation
mechanism. While the expression of the attenuation length scale may vary when
considering attenuation mechanisms other than radiative damping (such as
viscous dissipation or turbulent viscosity accounting for wave-breaking
effects), the qualitative nature of the solution is robust.

Equations (1) and (2) are solved on a uniform grid with $200$ levels, using a
3rd order Adams-Bashforth time stepping scheme with a timestep $\Delta
t=10^{-4}\tau$, where $\tau$ is the characteristic timescale for mean-flow
reversal222In the case of multi-modal forcing, the timestep is set by the
smallest of the two characteristic timescales.. An implicit diffusion scheme
is used for the eddy-viscous term. Injection of momentum for a given wave is
tapered above a value of $0.95c_{i}$.

##  3 Results

###  3.1 Periodicity recovery for increasing attenuation length scale

We begin by revisiting the case of a monochromatic standing wave forcing,
following the framework initially proposed by Plumb (1977). In this context,
we eliminate perturbations of the background state by setting
$\overline{u^{\prime}w^{\prime}}_{p}$ and $F_{p}$ to zero in Eqs (1) and (2).
For simplicity, we temporarily drop the index $b$ of the background state, as
we are dealing with a single pair of waves. This context reduces the model’s
dynamics to two dimensionless control parameters:

$Re=\frac{Fh}{\nu c}\,,\quad\widehat{h}=\frac{h}{H}.$ (5)

Figure 1: Bifurcation diagram for a monochromatic standing wave forcing. (a)
Growth rate of the linearized system for $\widehat{h}=1/3.5$. Black lines with
starred indices: oscillatory modes ($\mathrm{Im}(\sigma)\neq 0$). Orange
lines: non-oscillatory modes ($\mathrm{Im}(\sigma)=0$). Indices refer to the
vertical structure of the modes. (b) Dynamical regimes of the quasilinear
model in $(\widehat{h},Re)$ space, using Poincaré map metric. Black lines
indicate the appearance of an unstable oscillatory mode in the linearized
system, orange lines correspond to non-oscillatory ones. Green dots correspond
to those in the panel (a).

Figure 1b presents the bifurcation diagram obtained by varying of both
dimensionless control parameters. Bifurcation diagrams are constructed using
Poincaré maps, where values of $\overline{u}$ near the top of the domain that
intersect with $\overline{u}=0$ near the bottom are grouped into bins. The
shading in the Figure corresponds to the total number of bins populated in a
given experiment. Experiments with only one bin indicate a stable resting
state, those with two bins represent a QBO-like periodic oscillation, while
more than two bins signify the progression toward aperiodicity and chaotic
behavior (Renaud et al., 2019). The transition from a resting state to a QBO-
like oscillation is referred to as the first bifurcation, whereas the
transition from a QBO-like oscillation to higher oscillatory modes is termed
the second bifurcation.

In standard textbook treatments, the assumption is often made of a semi-
infinite stratosphere, where $\widehat{h}\rightarrow 0$, corresponding to the
left of Fig. 1b. Remarkably, for any fixed value of the forcing parameter
$Re$, the system transitions towards more periodic states as $\widehat{h}$
becomes sufficiently large. This implies that, for a finite domain depth, $H$,
modulations of wind reversal regimes can be achieved solely by varying the
attenuation length scale $h$, which is determined by the wave properties
originating in the upper troposphere, regardless of their forcing amplitude.

Focusing on the structure of bifurcation diagrams and considering the entire
range of $Re$ as $\widehat{h}$ increases, a new dynamical regime emerges,
which is not present in the semi-infinite assumption ($\widehat{h}\rightarrow
0$). This regime is characterized by bistability, represented by the light
grey shading in Figure 1b. Within this bistable regime, two steady non-zero
zonal wind solutions coexist, effectively breaking the east-west system
symmetry imposed by the wave generator. The introduction of this new dynamical
regime has a significant impact on the behavior of the first bifurcation,
which occurs between resting states and QBO-like oscillations. This transition
is marked by a shift from a supercritical Hopf bifurcation to a subcritical
one as $\widehat{h}$ surpasses a critical threshold, denoted as
$\widehat{h}_{c}$. For values of $\widehat{h}$ below $\widehat{h}_{c}$, the
system undergoes a first bifurcation from a state of rest to a QBO-like
oscillation. However, when $\widehat{h}$ exceeds $\widehat{h}_{c}$, the system
first transitions to a steady non-zero zonal wind profile before generating
QBO-like regimes. The system exhibits hysteresis behavior as $Re$ is varied
continuously across this regime. Similar transitions from supercritical to
subcritical Hopf bifurcations have been previously reported by Semin et al.
(2018) when varying other parameters of the quasilinear model, such as the
strength of an additional linear drag, see also Yoden and Holton (1988) in the
case of asymmetric wave forcing.

####  3.1.1 Predicting bifurcations with the linearized system

To gain insight into the dynamics explaining the bifurcation diagram of Fig.
1b, we linearize the model around a state of rest, following the approach of
Semin et al. (2018). Using nondimensional parameters:

$U=\frac{\overline{u}}{c}\,,\quad Z=\frac{z}{H}\,,\quad T=\frac{tF}{ch}\,,$
(6)

and assuming $U\ll 1$, we obtain, at the lowest order, the rescaled linear
system:

$\partial_{T}U=\frac{\widehat{h}^{2}}{Re}\partial_{ZZ}U+4\mathrm{e}^{-Z/\widehat{h}}\left(U-\frac{1}{\widehat{h}}\int_{0}^{Z}U\,\mathrm{d}Z^{\prime}\right)\,,$
(7)

with boundary conditions

$U=0\,\,\,\text{at}\,\,\,Z=0\,,\quad\partial_{Z}U=0\,\,\,\text{at}\,\,\,Z=1\,.$
(8)

This system can be turned into an eigenvalue problem,
$\sigma\widehat{U}=L\hat{U}$, by writing $U(z,t)=\widehat{U}\mathrm{e}^{\sigma
t}$. The eigenspectrum is spanned by solving numerically the integro-
differential operator $L$, corresponding to the right-hand side of Eq. (7). We
first focus on the classical scenario of a semi-infinite stratosphere,
characterized by $\widehat{h}\ll 1$, for which the real part of the eigenvalue
$\sigma$ is plotted in Fig. 1a as a function of $Re$. In this Figure, black
signifies that a mode is oscillatory ($\mathrm{Im}(\sigma)\neq 0$), while
orange signifies a mode that is not oscillatory ($\mathrm{Im}(\sigma)=0$).
Vertical modes are numbered based on the maximum number of zeros or nodes,
excluding boundary conditions (see Supplemental Materials). The point at which
an eigenvalue changes sign marks the onset of unstable exponential growth
(green dots in Fig. 1a). Within this semi-infinite domain context, the
threshold of the first bifurcation in Fig. 1b corresponds to the instability
of the first oscillating mode, denoted as 1⋆. Remarkably, the threshold of the
second and third bifurcations, as predicted in the quasilinear simulations,
closely aligns with the instability of the third and fifth oscillating mode,
referred to as 5⋆ and 9⋆. We conclude that the transition to chaos corresponds
to the excitation of unstable modes with an increasing number of alternating
horizontal jets in the vertical direction as the Reynolds number is increased.

The onset of instability of each modes is illustrated as black and orange
lines in the bifurcation diagram of Figure 1b. As $\widehat{h}$ increases, all
unstable eigenvalues become purely real, as the wave operator $L$
progressively becomes proportional to a Laplacian. This is why any black
branch (oscillatory modes) split into two orange branches (non-oscillatory)
modes as $\widehat{h}$ is increased in 1b. Consequently, the first QBO
bifurcation at large $\widehat{h}$ is different than in the semi-infinite case
($\widehat{h}\ll 1$). For values of $\widehat{h}$ larger than a critical
threshold denoted $\widehat{h}_{c}$, the emergence of the first unstable
eigenvalue corresponds to a pitchfork bifurcation of fixed points (from a
state of rest to bistability). The threshold of this first bifurcation
corresponds to the instability of the first stationary mode, denoted as 0 in
Fig. 1b. This bistable configuration becomes unstable at higher $Re$, with the
emergence of a limit cycle around the two fixed points. The transition occurs
before the second stationary unstable eigenvalue and is not captured by the
linearized dynamics around a state of rest, since this would instead require a
linearization around the fixed points. Surprisingly, the threshold of the
second bifurcation, originally captured by the instability of the oscillatory-
mode-5⋆ at $\widehat{h}\ll 1$, coincides with the instability of the
stationary-mode-4 as $\widehat{h}$ increases.

It is intriguing that the instability of the eigenmodes of the linearized
system is capable of capturing the emergence of successive bifurcations beyond
the first one given that they are computed around a state of rest.
Additionally, the connection between the instability of certain modes and
their impact on bifurcations, while others have limited effect, is also
somewhat puzzling. To fully understand these selection rules and explain the
predictive ability of the linearized dynamics around a state of rest for
subsequent bifurcations, a comprehensive Floquet analysis would be required.
However, conducting such an analysis falls beyond the scope of this paper.

Dimensional analysis helps explain why increasing $\widehat{h}$ promotes
periodicity. In a semi-infinite stratosphere where $\widehat{h}\rightarrow 0$,
$Re$ is the sole dimensionless parameter governing the system’s dynamics
Vallis (2017). When $\widehat{h}$ is large, a direct inspection of the
equation (7) shows that the appearance of new unstable eigenmodes is
controlled by $Re/\widehat{h}^{2}$. This amounts to defining a new Reynolds
number based on the domain height $H$ rather than on the attenuation length
scale $h$ and on the reduced forcing $F/\widehat{h}$ rather than $F$. In the
limit $\widehat{h}\gg 1$, bifurcations occur at finite values of
$Re/\widehat{h}^{2}$, which causes the range of $Re$ between each bifurcation
to widen as $\widehat{h}$ increases, as observed in Figure 1b.

We have seen that the emergence of aperiodic states involves the excitation of
higher vertical modes with multiple jets. During their nonlinear evolution,
some of these jets grow until reaching the upper-bound set by the critical
layers $\overline{u}=c$. In the case of a semi-infinite domain, the altitude
at which these critical layers form scales roughly with $h$. As $\widehat{h}$
becomes comparable to $1$ or larger, this altitude scales with the
stratosphere height $H$. More generally, for a given $Re$, increasing $h$
implies reduced wave damping, higher altitude wave propagation, and increased
vertical homogeneity, concomitantly with periodicity recovery. We hypothesize
that this vertical homogeneity in wave forcing promotes periodicity as it
hinders the development of higher vertical modes. This hypothesis will be
further substantiated in the case of multimodal forcing in Section 33.3.

###  3.2 A new metric to characterize dynamical regimes using vertical modes

As noted above, the onset of aperiodicity is inherently associated with the
instability of higher vertical modes of the linearized system, with the number
of alternating jets increasing as the forcing parameter grows. This excitation
of higher vertical modes with increasing forcing is manifested in the
quasilinear model, as illustrated using Hovmöller plots of the zonal winds in
Fig. 2a-d. Periodic oscillations are characterized by a maximum of one
downward propagating flow reversal in the vertical at any instant (or
1-node)333Boundary conditions are excluded from this definition.. In contrast,
strongly forced aperiodic regimes exhibit time periods with multiple downward
propagating nodes. Figure 2e depicts a bifurcation diagram of the histograms
employing Poincaré maps, for which the information is condensed as a number of
bins in Figure 1b. This is contrasted with Fig. 2f, where the ratio of the
number of flow reversals compared to that of the topmost grid point is plotted
as a function of altitude. As the forcing is increased, faster flow reversals
are observed at the bottom of the stratosphere compared to the top. These
higher vertical structures can be interpreted as an entanglement of multiple
oscillators with increasing frequency toward the ground (see Fig. 2d).
Importantly, each major transition in the bifurcation diagram of Fig. 2e is
closely linked to a regime transition in the vertical modes in Fig. 2f.

Figure 2: (a)-(d) Hovmöller diagrams for the wind speed profile at different
Reynolds numbers, for $\widehat{h}=1/3.5$. (e) Corresponding bifurcation
diagram computed with Poincaré metric. (f) Same bifurcation diagram computed
with the new metric $\mathcal{N}_{0}$, defined as the ratio between the number
of flow reversals on every level to that of the upper one.

Based on this observation, we propose a new metric for diagnosing QBO-like
regimes that focuses on the vertical structure of the oscillations, rather
than employing histograms of Poincaré maps (e.g., Renaud et al., 2019; Léard
et al., 2020). For a given total simulated time, this new metric computes the
ratio of the number of zeros counted in the bottom boundary layer (grey region
in Fig. 2f) compared to that counted near the top. This ratio is loosely
related to a vertical mode as it yields a time-averaged number of nodes in the
column, excluding periods without zonal wind reversal from this average. The
numbering of each vertical mode, explained in Supplemental Materials (Fig.
S2), is based on the maximum number of zeros in the column, excluding boundary
conditions. For oscillatory modes, we consider the maximum through a complete
oscillation.

Hereafter, states with a ratio of one, meaning a single node or mode-1, will
be referred to as QBO-like states. Notably, this metric is better suited to
describe the actual Earth QBO, which exhibits fluctuations around its average
period while remaining trapped in a mode-1 state prior to 2016. Yet, under
this new metric, both periodic and quasiperiodic solutions with mode-1
vertical structure are considered QBO-like. This results in a slight shift in
the value of the second bifurcation compared to the Poincaré metric,
corresponding to the quasiperiodic region (Fig. 2e,f).

This minor limitation is outweighed by a major advantage compared to Poincaré
sections: it exhibits a monotonic increase with both $Re$ and $h$. For
example, quasiperiodic regimes are often wedged between periodic and
frequency-locked regimes when increasing $Re$ (Fig. 2e). This is diagnosed as
a local maximum using Poincaré sections, whereas a monotonic increase is
observed using vertical modes. The monotonicity of the metric greatly helps in
interpreting the response of the system to a given change in parameter. In
particular, if increasing a parameter reduces the number of nodes, it can be
interpreted as favoring periodicity. Conversely, if the number of nodes
increases with an increasing parameter, it suggests a tendency towards
aperiodicity or more complex dynamics. This idea will greatly facilitate the
interpretation of the results in the following section, which involve
perturbations to the background QBO.

###  3.3 Periodicity recovery by perturbations to an aperiodic background
state

We now shift the focus to the case of a superposition of multiple
counterpropagating waves, which is arguably more realistic than monochromatic
wave forcing. As discussed above in the Model section, we consider, for
simplicity, the case of only two pairs of counterpropagating waves introduced
in Equations (1) and (2). Adding a pair of counterpropagating waves with
amplitude $F_{p}>0$, frequency $\omega_{p}$ and wavenumber $k_{p}$ implies
three new nondimensional parameters:

$\frac{F_{p}}{F_{b}}\,,\quad\frac{c_{p}}{c_{b}}\,,\quad\frac{h_{p}}{h_{b}}.$
(9)

Having already discussed in depth the role of the attenuation length scale in
the monochromatic case, we now constrain ourselves to the case $h_{p}=h_{b}\ll
H$ in the following discussion. This restriction does not compromise the
generality of the proposed mechanism for periodicity restoration. Below, we
describe how the introduction of a second wave exerts a stabilizing influence
on the system by preventing the excitation of higher vertical modes
responsible for aperiodicity.

In order to disentangle the role of the perturbation from that of the base
state, we set $F_{b}$ to a fixed value characterized by a mode-3 oscillation,
beyond the second bifurcation, which corresponds to the solution shown at Fig.
2c. We then conduct a series of experiments where we incrementally increase
the perturbation $F_{p}$ while keeping $F_{b}$ constant — thus increasing the
total forcing — for a wide range of phase speed ratios $c_{p}/c_{b}$. One
might intuitively anticipate that as the total forcing increases, the system
would transition towards higher vertical modes and chaos, as noted by Renaud
et al. (2019). Paradoxically, our results, depicted in Fig. 3a, reveal a
parameter space featuring instead two prominent regions where the flow
transitions back to a periodic QBO-like oscillation: one for $c_{p}<c_{b}$,
denoted as QBO(-), and another for $c_{p}>c_{b}$, denoted as QBO(+). Apart
from a small region in the upper left corner of Fig. 3a ($c_{p}\ll c_{b}$),
where higher vertical modes with chaotic reversals are excited, no significant
changes are observed elsewhere.

Figure 3: (a) Bifurcation diagram for the case of multiple wave forcing. See
text and Figure 2 for a definition of the metric $\mathcal{N}_{0}$. (b)
Maximal wind speed recorded in each regime. (c) Frequency of wind reversals
computed from a time series of the fastest oscillations near the lower
boundary in each regime. In each panel, thick black contours indicate isolines
of $\tau_{p}/\tau_{b}$.

To understand the primary drivers in these regions, we consider the maximum
velocity of the mean flow, displayed in Fig. 3b. When the background wave with
phase speed $c_{b}$ governs the system, the maximum speed of the mean flow is
$\overline{u}_{max}\approx c_{b}$. In most of the parameter space, regimes are
dominated by the background wave. There are two notable exceptions. Firstly,
the system undergoes an abrupt transition when $c_{p}\ll c_{b}$, where the
mean flow decelerates suddenly (Fig. 3b) as it becomes governed by the
perturbation associated with chaotic high vertical modes (Fig. 3a). Secondly,
a smooth transition occurs for $c_{p}>c_{b}$, where the mean flow gradually
accelerates as it increasingly falls under the influence of the perturbation.
This smooth transition for $c_{p}>c_{b}$ aligns with a decrease in the
frequency of wind reversals (as indicated in Fig. 3c) and is closely linked to
the QBO(+) region of periodicity recovery observed in Fig. 3a. This contrasts
with the region QBO(-), where the frequency of wind reversals increases while
no change in mean flow is observed. Below, we focus individually on each of
the periodic regions.

In the QBO(-) region, where the background wave governs the system, two types
of critical levels can emerge: (i) critical levels associated with the
background waves, which occur when the maximum wind speed approaches the phase
speed of the background waves ($\overline{u}_{\max}\approx\pm c_{b}$); and
(ii) critical levels linked to the perturbation waves, which occur when the
mean wind speed matches the horizontal phase speed of the perturbation waves
($\overline{u}=\pm c_{p}$).

To gain insight into the role of each type of critical level, we conduct an
experiment depicted in Fig. 4, where a perturbation is introduced at $t=0$
from the reference experiment using monochromatic forcing.

Figure 4: Transition from mode-3 to QBO-like oscillations with the addition of
a perturbation wave at $t=0$, corresponding to the orange star in Figure 3.
(a) Hovmöller diagram. Black and yellow contours trace the downward
propagation of critical layers of the perturbation and the background waves,
respectively. (b) and (c) Instantaneous vertical profiles of the mean flow
(green line) under monochromatic wave forcing, with blue and red dots
representing westward and eastward Reynolds stress derivatives, respectively
(dot size is proportional to amplitude). (d) and (e) Same for multiple wave
forcing, where Reynolds stresses are decomposed into background wave (at
$\overline{u}=\pm c_{b}$) and perturbation wave (at $\overline{u}=\pm c_{p}$)
contributions.

After just a few oscillation cycles following the introduction of the
perturbation wave, a QBO-like cycle is recovered (Fig. 4a), corresponding to
the experiment marked with an orange star in the Figure 3. In both the
monochromatic and multimodal cases, we analyze two snapshots typical of the
westward phase of the mean winds ($\overline{u}<0$), as well as their
associated wave forcing. The first snapshot is taken soon after the previous
wind reversal, while the second snapshot is right before the next wind
reversal (as shown in Fig. 4).

In both cases, the forcing from the background waves is similar at both times
(Fig. 4b-e). However, in the multimodal case, an additional eastward forcing
from the perturbation is initially introduced in the upper part of the
stratosphere (see $t_{3}$ in Fig. 4e). This high-altitude eastward forcing is
a consequence of the Doppler shift’s impact on the local attenuation length
scale444The Doppler-shifted attenuation length scale of the perturbation in
the presence of a mean flow, where $\overline{u}\sim c_{b}$ is much larger
than the phase speed $c_{p}\ll c_{b}$, scales as $(c_{b}/c_{p})^{2}h$ (as
explained in Section 2)., which allows the perturbation wave with an eastward
phase speed to propagate upward with minimal attenuation until it reaches the
upper part of the stratosphere, where the winds are much weaker. This
efficient transport leads to high-altitude deposition of eastward momentum
from the perturbation and to the development of a fast downward propagating
critical layer that closely follows the zero-wind line (Fig. 4d,e). This fast-
moving critical layer prevents the excitation of secondary westward jets or
higher vertical modes responsible for aperiodicity, by smoothing the vertical
wind profile (compare Fig. 4c,e).

Periodicity recovery results from a synchronization mechanism between the fast
descent of the perturbation-driven zero-wind line and the background-driven
mean-wind oscillation. Synchronization occurs when the perturbation-driven
critical layers are sufficiently faster than the background-driven critical
layers to prevent the growth of secondary fronts forced by the background
wave. The typical time scale associated with the descent of background-driven
critical layers is $\tau_{b}=c_{b}h_{b}/F_{b}$, while that of the
perturbation-driven critical layers is $\tau_{p}=c_{p}h_{p}/F_{p}$, assuming
that all the momentum of the perturbation wave is efficiently used to move the
location where $\overline{u}=c_{p}$ downward. A periodic state is thus
expected when $\tau_{p}<\tau_{b}$.

In each panel of Fig. 3, lines representing the ratio
$\tau_{p}/\tau_{b}=0.1,1,10$ are displayed. The region QBO(-) is well
represented by the area between $\tau_{p}/\tau_{b}=0.1$ and $1$, as expected
from the above analysis. In this region a periodic behavior with a QBO-like
vertical structure emerges from the combination of perturbed and background
wave forcing, which would independently lead to nonperiodic behavior. This
synchronization mechanism ceases to be effective when the critical levels
associated with the perturbation are fast enough to govern the oscillation
itself, which becomes chaotic at $\tau_{p}/\tau_{b}\approx 0.1$.

The occurrence of the region QBO(+) coinciding with the line of
$\tau_{p}/\tau_{b}=10$ suggests that a similar synchronization mechanism is at
play in this second regime of periodicity recovery. However, in this case, the
roles of the background and perturbation are somewhat reversed: critical
levels associated with the background forcing occur progressively at higher
altitudes as the maximum mean flow gradually accelerates when the oscillation
falls under the influence of the perturbation (see Fig. S3). Additionally, the
transfer of momentum from the perturbation to the background state is much
less efficient than in QBO(-) region due to a compensation effect between the
wave momentum flux induced by the eastward-moving and westward-moving
perturbations, which prevents these perturbations from reaching critical
levels. An extreme case is observed when $c_{p}/c_{b}\gg 1$, where the
perturbation has virtually no effect on the wind reversals. This is evident in
Fig. 3, where QBO properties remain unchanged beyond $c_{p}=10c_{b}$. This
occurs because the wave momentum flux induced by the eastward-moving
perturbation becomes nearly equal to the contribution from the westward-moving
perturbation when $\overline{u}/c_{p}\ll 1$, as described by Eq. (2). For
finite values of $c_{p}/c_{b}>1$, the primary effect of the perturbation is to
shift the maxima of the wind profile to values greater than the critical level
$c_{b}$. Assuming that most of the descent of these maxima remains dominated
by contributions from the background wave forcing, the period of the QBO can
be estimated as $\overline{u}_{\mathrm{max}}h_{b}/F_{b}$, with
$\overline{u}_{\mathrm{max}}$ being an increasing function of $F_{p}$. This
explains the decrease in frequency associated with the QBO(+) region in Fig.
3c.

In summary, perturbation waves with phase speeds $c_{p}<c_{b}$ lead to a
faster descent of the QBO front to the bottom of the stratosphere, resulting
in a shorter period of wind reversals. On the other hand, waves with
$c_{p}>c_{b}$ amplify the amplitude of the QBO front, leading to an extended
period of reversals. In both cases, the presence of a wave with a slower phase
speed prevents the growth of unstable modes of the larger wave due to
resonance between two timescales: any local extremum in the velocity profile
is shifted downward before its amplitude reaches a critical level. This
effectively prevents the excitation of higher vertical modes responsible for
aperiodicity and explains the observed periodicity recovery regions
characterized by a vertical mode-1 akin to the QBO. This phenomenon is all the
more surprising as it emerges from the combination of perturbed and background
wave forcing that independently would lead to nonperiodic reversals.

##  4 Conclusion

Building upon the insights gained from the classical monochromatic scenario,
we have proposed a physical mechanism for intrinsic synchronization in
internal-wave–mean-flow interactions with multiple-wave forcing. The key
concept is that raising the altitude at which wave momentum is deposited
promotes periodic behavior. This can be achieved by modifying the attenuation
length scale in the monochromatic case or through multimodal wave forcing,
where the mean flow generated by the faster wave significantly alters the
Doppler-shifted attenuation length scale of the slower wave. To maintain
simplicity, we focused on a scenario involving only two pairs of
counterpropagating waves, which lays the groundwork for interpreting more
complex configurations involving a wider distribution of waves generated in
the upper troposphere.

Using a quasilinear model similar to ours, forced with a Gaussian frequency
spectrum, Léard et al. (2020) found that increasing the variance of the
distribution, while maintaining a constant total forcing, favors periodic
behavior. The resulting QBO-like periodic regime exhibited significantly
higher velocities compared to the chaotic monochromatic forcing, leading the
authors to suggest that high-frequency waves, given a fixed horizontal wave
number, played a crucial role in shaping QBO characteristics. Translating
their observation into our specific framework, we associate this recovery of
periodicity to the region denoted as QBO(+), for which a similar zonal winds
amplification is observed. Importantly, our study unveiled another region of
periodicity recovery, QBO(-), where the phase speed of the perturbation wave
is smaller than the background wave phase speed. In contrast to the
periodicity recovery in the QBO(+) region, the maximum zonal wind remains
unchanged in QBO(-) since the mean flow is governed by the background wave in
this region. Thus, beyond modulations in the forcing frequency, our study
highlights the critical influence of perturbation phase speed variations,
achieved by simultaneous adjustments in frequency and wavenumber.

In the context of a more realistic wave representation, the respective roles
of the different classes of equatorial waves remains a subject of debate
(e.g., Pahlavan et al., 2021; Holt et al., 2022). While planetary-scale waves
are widely acknowledged as a central driver of the oscillation, previous
numerical studies of the atmospheric QBO have emphasized the need to
parameterize the effects of small-scale internal gravity waves in order to
induce a QBO-like regime (Giorgetta et al., 2006; Lott and Guez, 2013).

In the framework of our study, the primary standing wave serve as a
rudimentary representation of the influence exerted by planetary Rossby,
Yanai, and Kelvin waves, while the perturbation standing wave can be
understood as the result of smaller-scale and higher-frequency internal waves
generated by localized convection events.

Our study demonstrates the crucial role of the zonal phase speed ratio,
denoted as $c_{p}/c_{b}$, which characterizes an important relationship
between background and perturbation waves. The phase speed of convectively
coupled equatorial Kelvin waves, estimated at about $15$ m/s, serves as a
reference value for $c_{b}$. These waves exhibit a vertical wavelength of
$\pi/H_{t}$, where $H_{t}$ represents the tropospheric height based on the
equatorial Rossby radius of deformation, $\sqrt{c_{b}/\beta}$. We assume that
perturbations arise from localized convection events that generate internal
gravity waves with similar vertical length scales, namely $\pi/H_{t}$.
However, their horizontal length scales are significantly smaller than the
equatorial radius of deformation. This results in a perturbation phase speed
with an absolute value close to $c_{b}$. Nonetheless, the projection of this
phase speed in the zonal direction is on average smaller than $c_{b}$ in the
case of oblique ray propagation with respect to the equator (Kim et al.,
2023). Synchronization induced by such oblique waves would then correspond to
the region QBO(-) highlighted in our study.

To establish this interpretation on firmer ground, we will need to generalize
our study to a 3D configuration on the stratified equatorial beta plane (Plumb
and Bell, 1982). In this expanded model, the attenuation length scale for
Rossby and Yanai waves will differ from the scales used in our study, which
concentrated solely on internal gravity waves. It should be reminded that we
also neglected additional effects such as coupling with an external low
frequency oscillator (Read and Castrejón-Pita, 2012), stochasticity
(Saravanan, 1990; Wedi and Smolarkiewicz, 2006), wave-wave interactions
(Couston et al., 2018), transient behavior of the wave field (Dunkerton,
1981). These effects are currently under investigation for efficient
parameterizations of unresolved internal gravity waves in climate model (e.g.,
Achatz et al., 2023). Ensuring the robustness of the periodicity recovery
mechanism described in our study will require thorough testing against these
influences.

###### Acknowledgements.

We acknowledge Cerasela Calugaru and the whole team from PSMN/CBP at ENS de
Lyon for computing facilities. We warmly thank Nicolas Perez for his valuable
help in performing the linear analysis with the Dedalus solver. We also thank
Louis Couston and Jason Reneuve for useful discussions on wind reversals with
multiple wave forcing. Louis-Philippe Nadeau is partly supported through NSERC
Award RGPIN-2022-04306. The model used for this study is available on GitHub
at https://github.com/xavier-chartrand/qbo1d.git. Model outputs are archived
on PSMN clusters and remain available upon request from XC.

## References

  * Achatz et al. (2023) Achatz, U., Y. H. Kim, and G. S. Voelker, 2023: Multi-Scale Dynamics of the Interaction Between Waves and Mean Flows: From Nonlinear WKB Theory to Gravity-Wave Parameterizations in Weather and Climate Models. arXiv preprint arXiv:2310.07334.
  * Anstey et al. (2022) Anstey, J. A., and Coauthors, 2022: Impacts, processes and projections of the quasi-biennial oscillation. Nature Reviews Earth & Environment, 3 (9), 588–603.
  * Baldwin et al. (2001) Baldwin, M. P., and Coauthors, 2001: The quasi-biennial oscillation. Reviews of Geophysics, 39 (2), 179–229.
  * Bardet et al. (2022) Bardet, D., A. Spiga, and S. Guerlet, 2022: Joint evolution of equatorial oscillation and interhemispheric circulation in Saturn’s stratosphere. Nature Astronomy, 6 (7), 804–811.
  * Couston et al. (2018) Couston, L. A., D. Lecoanet, B. Favier, and M. Le Bars, 2018: Order out of chaos: slowly reversing mean flows emerge from turbulently generated internal waves. Physical review letters, 120 (24), 244 505.
  * Dunkerton (1981) Dunkerton, T. J., 1981: Wave transience in a compressible atmosphere. Part I: Transient internal wave, mean-flow interaction. Journal of Atmospheric Sciences, 38 (2), 281–297.
  * Dunkerton and Delisi (1997) Dunkerton, T. J., and D. P. Delisi, 1997: Interaction of the quasi-biennial oscillation and stratopause semiannual oscillation. Journal of Geophysical Research: Atmospheres, 102 (D22), 26 107–26 116.
  * Fletcher et al. (2017) Fletcher, L. N., and Coauthors, 2017: Disruption of Saturn’s quasi-periodic equatorial oscillation by the great northern storm. Nature Astronomy, 1 (11), 765–770.
  * Giorgetta et al. (2006) Giorgetta, M. A., E. Manzini, E. Roeckner, M. Esch, and L. Bengtsson, 2006: Climatology and forcing of the quasi-biennial oscillation in the MAECHAM5 model. Journal of Climate, 19 (16), 3882–3901.
  * Holt et al. (2022) Holt, L. A., and Coauthors, 2022: An evaluation of tropical waves and wave forcing of the QBO in the QBOi models. Quarterly Journal of the Royal Meteorological Society, 148 (744), 1541–1567.
  * Holton and Lindzen (1972) Holton, J. R., and R. S. Lindzen, 1972: An updated theory for the quasi-biennial cycle of the tropical stratosphere. Journal of the Atmospheric Sciences, 29 (6), 1076–1080.
  * Kim and MacGregor (2001) Kim, E. j., and K. B. MacGregor, 2001: Gravity wave-driven flows in the solar tachocline. The Astrophysical Journal, 556 (2), L117.
  * Kim et al. (2023) Kim, Y. H., G. S. Voelker, G. Bölöni, G. Zängl, and U. Achatz, 2023: Crucial Role of Obliquely Propagating Gravity Waves in Tropical Stratospheric Circulation. arXiv preprint arXiv:2309.15301.
  * Léard et al. (2020) Léard, P., D. Lecoanet, and M. Le Bars, 2020: Multimodal Excitation to Model the Quasibiennial Oscillation. Physical review letters, 125 (23), 234 501.
  * Lott and Guez (2013) Lott, F. G., and L. Guez, 2013: A stochastic parameterization of the gravity waves due to convection and its impact on the equatorial stratosphere. Journal of Geophysical Research: Atmospheres, 118 (16), 8897–8909.
  * Match and Fueglistaler (2019) Match, A., and S. Fueglistaler, 2019: The buffer zone of the quasi-biennial oscillation. Journal of the Atmospheric Sciences, 76 (11), 3553–3567.
  * Match and Fueglistaler (2021) Match, A., and S. Fueglistaler, 2021: Large internal variability dominates over global warming signal in observed lower stratospheric QBO amplitude. Journal of Climate, 34 (24), 9823–9836.
  * Newman et al. (2016) Newman, P. A., L. Coy, S. Pawson, and L. R. Lait, 2016: The anomalous change in the QBO in 2015–2016. Geophysical Research Letters, 43 (16), 8791–8797.
  * Osprey et al. (2016) Osprey, S. M., N. Butchart, J. R. Knight, A. A. Scaife, K. Hamilton, J. A. Anstey, V. Schenzinger, and C. Zhang, 2016: An unexpected disruption of the atmospheric quasi-biennial oscillation. Science, 353 (6306), 1424–1427.
  * Pahlavan et al. (2021) Pahlavan, H. A., J. M. Wallace, Q. Fu, and G. N. Kiladis, 2021: Revisiting the quasi-biennial oscillation as seen in ERA5. Part II: Evaluation of waves and wave forcing. Journal of the Atmospheric Sciences, 78 (3), 693–707.
  * Plumb (1977) Plumb, R. A., 1977: The interaction of two internal waves with the mean flow: Implications for the theory of the quasi-biennial oscillation. Journal of the Atmospheric Sciences, 34 (12), 1847–1858.
  * Plumb and Bell (1982) Plumb, R. A., and R. C. Bell, 1982: A model of the quasi-biennial oscillation on an equatorial beta-plane. Quarterly Journal of the Royal Meteorological Society, 108 (456), 335–352.
  * Rajendran et al. (2016) Rajendran, K., I. M. Moroz, P. L. Read, and S. M. Osprey, 2016: Synchronisation of the equatorial QBO by the annual cycle in tropical upwelling in a warming climate. Quarterly Journal of the Royal Meteorological Society, 142 (695), 1111–1120.
  * Read and Castrejón-Pita (2012) Read, P. L., and A. A. Castrejón-Pita, 2012: Phase synchronization between stratospheric and tropospheric quasi-biennial and semi-annual oscillations. Quarterly Journal of the Royal Meteorological Society, 138 (666), 1338–1349.
  * Renaud et al. (2019) Renaud, A., L. P. Nadeau, and A. Venaille, 2019: Periodicity disruption of a model quasibiennial oscillation of equatorial winds. Physical review letters, 122 (21), 214 504.
  * Renaud and Venaille (2020) Renaud, A., and A. Venaille, 2020: On the Holton–Lindzen–Plumb model for mean flow reversals in stratified fluids. Quarterly Journal of the Royal Meteorological Society, 146 (732), 2981–2997.
  * Saravanan (1990) Saravanan, R., 1990: A multiwave model of the quasi-biennial oscillation. Journal of Atmospheric Sciences, 47 (21), 2465–2474.
  * Semin et al. (2018) Semin, B., N. Garroum, F. Pétrélis, and S. Fauve, 2018: Nonlinear saturation of the large scale flow in a laboratory model of the quasibiennial oscillation. Physical review letters, 121 (13), 134 502.
  * Showman et al. (2019) Showman, A. P., X. Tan, and X. Zhang, 2019: Atmospheric Circulation of Brown Dwarfs and Jupiter-and Saturn-like Planets: Zonal Jets, Long-term Variability, and QBO-type Oscillations. The Astrophysical Journal, 883 (1), 4.
  * Vallis (2017) Vallis, G. K., 2017: Atmospheric and oceanic fluid dynamics. Cambridge University Press.
  * Wedi and Smolarkiewicz (2006) Wedi, N. P., and P. K. Smolarkiewicz, 2006: Direct numerical simulation of the Plumb–McEwan laboratory analog of the QBO. Journal of the atmospheric sciences, 63 (12), 3226–3252.
  * Yoden and Holton (1988) Yoden, S., and J. R. Holton, 1988: A new look at equatorial quasi-biennial oscillation models. Journal of Atmospheric Sciences, 45 (19), 2703–2717.

Supplemental Materials to: Recovering Quasi-Biennial Oscillations from Chaos

## I. Computation of the eigenvalue problem

In the following, we present the approach to solve the eigenvalue problem
arising from the Holton-Lindzen-Plumb model when linearized around a state of
rest, along with a detailed picture of eigenvalues $\sigma$ and eigenmodes
$\widehat{U}$. Recalling our choice of dimensionless parameters, namely
$U=\overline{u}/c$, $Z=z/H$ and $T=tF/(ch)$, we search for solutions
$U\propto\widehat{U}(Z)\exp(\sigma T)$ that satisfy

$\sigma\widehat{U}=L\widehat{U}\qquad\text{where}\qquad
L\widehat{U}=\frac{\widehat{h}^{2}}{Re}\partialderivative[2]{z}\widehat{U}+4\mathrm{e}^{-Z/\widehat{h}}\quantity(\widehat{U}-\frac{1}{\widehat{h}}\int_{0}^{Z}\widehat{U}\,\differential
Z^{\prime})\,.$ (S1)

By defining $\widehat{Q}=\partial_{Z}\widehat{U}$ and
$\widehat{U}=\partial_{Z}\widehat{F}$, Eq. (S1) can be written as a set of
three first-order homogeneous equations:

$\matrixquantity(\partialderivative{Z}&\frac{Re}{\widehat{h}^{2}}\quantity(4\mathrm{e}^{-Z/\widehat{h}}-\sigma)&-\frac{4Re\mathrm{e}^{-Z/\widehat{h}}}{\widehat{h}^{3}}\\\
1&\partialderivative{Z}&0\\\ 0&1&\partialderivative{Z}\\\
)\matrixquantity(\widehat{Q}\\\ \widehat{U}\\\
\widehat{F})=\matrixquantity(0\\\ 0\\\ 0)\,.$ (S2)

The no-slip and stress-free boundary conditions are imposed at the bottom and
top, respectively

$\widehat{F}(0)=0\,,\quad\widehat{U}(0)=0\,,\quad\widehat{Q}(1)=0\,.$ (S3)

The eigenspectrum of the matrix operator on the LHS of Eq. (S2) is solved
numerically using Dedalus along with the Eigentools library (see Burns et al.
(2020); Oishi et al. (2021)). The solutions are computed on a Chebyshev
polynomial basis with $nz=128$ grid points. Note that the operator $L$ is not
self-adjoint, therefore the eigenspectrum cannot be directly considered as a
natural basis to linearly decompose any solutions of the quasilinear model.
Nevertheless, the properties of the eigenspectrum are helpful to understand
the appearance of various bifurcations and each type of oscillation resolved
by the quasilinear model.

Figure S1 compares bifurcation diagrams with solutions of the eigenvalue
problem for two different values of $\widehat{h}$. At small $\widehat{h}$, the
first unstable eigenmode is oscillatory, with a non-zero imaginary part
(panels a-b-c). This corresponds to a supercritical Hopf bifurcation for the
rest state in the quasi-linear system. At large $\widehat{h}$, the unstable
eigenmodes are purely real (panels d-e-f). The first bifurcation corresponds
to a pitchfork bifurcation leading to nontrivial stationary states, breaking
the east-west symmetry in the quasilinear system. QBO-like states emerge in
that case from secondary instabilities of those non-trivial stationary states,
which cannot be detected from the operator obtained by linearizing around a
state of rest. There is an intermediate case
$\widehat{h}\in[\widehat{h}_{c},\widehat{h}_{c2}]$, where the first and second
unstable eigenmodes are real, while the third unstable mode is oscillatory.
This corresponds to a subcritical Hopf bifurcation for the base state at rest.
The critical value $\widehat{h}=\widehat{h}_{c}$ thus corresponds to a
transition from a supercritical to a subcritical Hopf bifurcation.

The numbering of each vertical mode, shown in Fig. S2, is based on the maximum
number of zeros in the column, excluding boundary conditions. For oscillatory
modes, we consider the maximum through a complete oscillation.

## II. Periodicity recovery in the QBO(+) region

In this section, we show the periodicity recovery in the region QBO(+)
($c_{p}>c_{b}$) using an experiment similar to that shown in Fig. 4 (see Fig.
S3). For this case, we selected $c_{p}/c_{b}=4.4$ and $F_{p}/F_{b}=0.8$. A
synchronization mechanism analogous what’s observed in the QBO(-) occurs
between critical layers of the background and the perturbation waves albeit
with their role reversed. Here, the perturbation enhances the magnitude of the
westward winds at $t_{3}$, leading to an injection of momentum by the
background wave at higher altitudes. The rapid descent of the background
wave’s critical layer prevents the excitation of higher vertical modes,
responsible for aperiodicity. A major difference with the periodicity recovery
in the region QBO(-) is that the perturbation doesn’t entirely form a critical
layer in this case. This suggests that critical layers of the fastest waves
are not a necessary ingredient of this synchronization mechanism.

## References

Burns, K. J., G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, 2020:
Dedalus: A flexible framework for numerical simulations with spectral methods.
Physical Review Research, 2 (2), 023 068.
Oishi, J. S., K. J. Burns, S. E. Clark, E. H. Anders, B. P. Brown, G. M.
Vasil, and D. Lecoanet, 2021: eigentools: A Python package for studying
differential eigenvalue problems with an emphasis on robustness. Journal of
Open Source Software, 6 (62), 3079.

Fig. S 1: Comparison of quasilinear simulations with $\widehat{h}=1/3.5$ (left
column) and $\widehat{h}=1/1.25$ (right column), to their linearized version.
(a), (d) Poincaré bifurcation diagram of quasilinear simulations. (b), (e) the
real and (c), (f) the imaginary part of the eigenvalue $\sigma$. Fig. S 2:
Indexation of eigenmodes by the number of nodes (colored dots). Real part of
$\widehat{U}$ for oscillatory modes indexed 1⋆, 3⋆ and 5⋆ presenting 1, 3 and
5 nodes, respectively. The same indexation applies for stationary modes. Fig.
S 3: Same as Figure 4, but for the case $c_{p}>c_{b}$. Black contours indicate
critical layers of the background wave.