Generation of zonal flows through symmetry breaking of statistical homogeneity

In geophysical and plasma contexts, zonal flows (ZFs) are well known to arise out of turbulence. We elucidate the transition from homogeneous turbulence without ZFs to inhomogeneous turbulence with steady ZFs. Starting from the equation for barotropic flow on a β plane, we employ both the quasilinear approximation and a statistical average, which retains a great deal of the qualitative behavior of the full system. Within the resulting framework known as CE2, we extend recent understanding of the symmetry-breaking zonostrophic instability and show that it is an example of a Type Is?> instability within the pattern formation literature. The broken symmetry is statistical homogeneity. Near the bifurcation point, the slow dynamics of CE2 are governed by a well-known amplitude equation. The important features of this amplitude equation, and therefore of the CE2 system, are multiple. First, the ZF wavelength is not unique. In an idealized, infinite system, there is a continuous band of ZF wavelengths that allow a nonlinear equilibrium. Second, of these wavelengths, only those within a smaller subband are stable. Unstable wavelengths must evolve to reach a stable wavelength; this process manifests as merging jets. These behaviors are shown numerically to hold in the CE2 system. We also conclude that the stability of the equilibria near the bifurcation point, which is governed by the Eckhaus instability, is independent of the Rayleigh–Kuo criterion.


Introduction
Zonal flows (ZFs) are azimuthally symmetric, banded, shear flows that can be spontaneously generated from turbulence. Their prominent geophysical manifestations have inspired study for decades [1][2][3]. ZFs have also been discovered in magnetized plasmas, and these flows are believed to play an important role in regulation of turbulence and turbulent transport [4,5]. Further cementing their universality, ZFs as well as zonal magnetic fields have been seen in astrophysical simulations, specifically of accretion disc turbulence driven by the magnetorotational instability [6,7].
ZFs remain incompletely understood, even regarding the basic question of the jet width. Various attempts have been made to relate the jet width or spacing to length scales that emerge from the vorticity equation by heuristically balancing the magnitudes of the Rossby wave term and the nonlinear (NL) advection term. Those scales include the Rhines scale β = L U ( / ) R 1/2 , where U is the rms velocity and β is the northward gradient of the Coriolis parameter [1], and a related transitional scale, ε β = ε L ( / ) 3 1/5 [8][9][10], where ε is the energy input. In simple situations, the total energy can be estimated as ε μ = E /2 , where μ is the friction. By using = U E 1 2 2 , one can write the Rhines scale in terms of external parameters, which often gives good agreement with the ZF length scale. A Rhines-like length scale is also obtained from arguments based on potential vorticity staircases [11,12]. Despite these advances, there has not yet been a theory of how one would determine the ZF length scale systematically from the equations of motion, nor even a framework of what such a theory might look like.
The question of jet length scale is clearly related to the phenomenon of merging jets. Coalescence of two or more jets is ubiquitous in numerical simulations [13][14][15]. The merging process typically occurs during the initial transient period before a statistically steady state is reached, although in some cases occasional merging and branching may persist. Merging occurs as part of a dynamical process through which the ZF reaches its preferred length scale, but there has been limited theoretical understanding of it so far [16,17].
Our present work addresses these questions in the context of stochastically forced barotropic flow on a β plane, a model of fluid turbulence in a rotating system [8]. This rather simple model excludes a myriad of physics that arises in realistic situations. The simplicity is an advantage as it allows us to focus on the elements it does retain, in particular the emergence of steady ZFs driven by turbulence. It is known that the barotropic vorticity equation is mathematically very similar to the generalized Hasegawa-Mima equation [18,19] for electrostatic potential, a model of magnetized plasma turbulence in the presence of a background density gradient. The methodology used here is equally applicable to the Hasegawa-Mima equation [20].
Our paper extends the work of Srinivasan and Young [15] in providing a firm analytic understanding of ZF generation and equilibrium. Our study employs the quasilinear (QL) approximation, which neglects certain nonlinear interactions. The use of this approximation is motivated by the fact that it retains much of the generic ZF behavior. A statistical average yields the CE2 (second-order cumulant) framework. CE2 has been used in many studies of the dynamics of ZF-turbulence interactions [21][22][23][24][25][26][27][28]. CE2 is mathematically equivalent to SSST or S3T (Stochastic Structural Stability Theory), which is the name used by Farrell, Ioannou, Bakas, and Constantinou [21][22][23][24][25]. Within CE2, it was discovered that homogeneous turbulence can undergo a symmetry-breaking instability leading to ZF generation [21], which was then calculated in detail analytically [15]. This instability has been named zonostrophic instability (ZI). The mechanism of ZI has been studied by Bakas and Ioannou [29].
Using CE2, we show from the statistical point of view that the bifurcation that generates steady ZFs is an example of what is called a Type I s instability within the pattern formation literature [30][31][32][33][34]. Therefore, many of the general properties of pattern-forming systems apply to these ZFs. Near the point of bifurcation the system obeys a classic amplitude equation, sometimes called the real Ginzburg-Landau equation. Two important analytic results about ZFs follow from the properties of the amplitude equation. First, the ZF wavelength is not unique. Indeed, in an idealized, infinite system, any wavelength within a certain continuous band corresponds to a steady-state solution. Nonuniqueness of the number of jets or the jet wavelength has in fact been observed numerically in direct numerical simulations [35] and in CE2 simulations [21]. Second, of these wavelengths, only those within a smaller subband are linearly stable. Unstable wavelengths must evolve to reach a stable wavelength. For short (long) wavelength unstable jets, this process manifests as merging (branching) jets.
A fundamental dimensionless parameter controlling the ZF dynamics is γ ε β μ = − 1/4 1/2 5/4 [35]. This parameter is related to the zonostrophy index = β ε R L L / R by γ ≈ β R 5 [36]. Also, if one were to modify the definition of the small parameter α of [37] such that the normalization length scale were the Rhines scale L R rather than the size of the domain, then γ α = −1 . Since it has been shown were in [37] that it is α (and α 1/2 ) that naturally appear in the normalized equations of motion, we opt to use γ instead of β R as the descriptive parameter. This paper is structured as follows. Section 2 introduces a phenomenological model of the bifurcation. As a zero-dimensional model, it orients the reader before the plunge into the full problem with its complexity of spatial dependence. Section 3 introduces the fundamental equation of motion, the (equivalent) barotropic vorticity equation, followed by the QL approximation and the statistical formulation known as CE2. A review of ZI is provided in section 4. Then, in section 5 (with some details in the appendix), we perform a full bifurcation analysis into the regime of nonlinearly interacting eddies and ZFs. A number of important phenomena appear even though the bifurcation calculation is strictly valid only in a limited domain of parameter space. To build on this understanding, we numerically solve for the CE2 steady states in section 6 and their stability in section 7. Section 8 explores the ZF wave number selection problem. Finally, section 9 contains discussion of the results and our conclusions.

Phenomenological bifurcation model of zonostrophic instability
A zero-dimensional phenomenological model illustrates some of the key features of ZI and the bifurcation to a state with ZFs. The system is a variant of another treatment which models the appearance of shear flows in the L-H transition in plasmas [38]. However, the model we present more closely mirrors the structure and behavior of the CE2 equations. The model includes three interacting degrees of freedom: the homogeneous, or spatial average, part of the fluctuation covariance W h ; the inhomogeneous, or deviation from the spatial average, part of the fluctuation covariance W i ; and the ZF amplitude (not covariance) z. Both W i and z may be positive or negative. The model is given by The structure of the model reflects that of the CE2 equations (6) in several ways. To affect the homogeneous part of the turbulence, the ZF interacts only with the inhomogeneous part. Similarly, the ZF interacts with the homogeneous part to affect the inhomogeneous part. (CE2 also contains an interaction between the ZFs and the inhomogeneous part to affect the inhomogeneous part; this is neglected here, as in ZI analysis.) Finally, it is the inhomogeneous part of the turbulence that is responsible for driving steady ZFs. The above model neglects the eddy self-nonlinearities as in CE2. The appearance of the same coefficient α inẆ h and inż reflects the conservation by the NL interactions of an energy-like quantity We take all of the coefficients μ ν α η F , , , , to be positive. The model allows a homogeneous equilibrium, in which forcing F is balanced by dissipation μ, and for which W i and z are zero. The homogeneous equilibrium is unstable if Increasing the forcing or decreasing the dissipation tends to make the homogeneous equilibrium more zonostrophically unstable, which is characteristic of the more rigorous analysis.
When the homogeneous equilibrium goes unstable, it connects to a stable inhomogeneous equilibrium at . Furthermore, there are actually two symmetric solutions, with either sign of z and W i . There is thus a supercritical pitchfork bifurcation; this feature is also present in the complete model, but the discrete → − z z symmetry becomes a continuous symmetry associated with translational invariance.
The model demonstrates some of the qualitative features of ZI, although in simplifying it we have tossed out spatial dependence. Spatial dependence makes the problem both immensely more complicated and immensely more interesting. The CE2 equations contain the full spatial dependence and will be introduced in section 3.4.

Barotropic vorticity equation
Our starting point is the (equivalent) barotropic vorticity equation on a β plane. We focus on results for the infinite deformation radius case, though we allow for finite deformation radius (the case of finite deformation radius is treated in a plasma context in [20]). The equation is 10 3 , and h = 4. The covariance of the forcing is taken to be πε δ , and zero otherwise. We take = k 1 f , δ = k 1 8 , and ε = 1. Periodic boundary conditions are used. During the transient period, merging jets are observed, while in the late time, a statistically steady state is reached with steady jets. Although the jets may not be completely unwavering, this is a useful idealization that we make.

Quasilinear approximation
We restrict ourselves to the QL approximation of this system. To obtain the QL equations, one performs an eddy-mean decomposition, given by decomposing all fields into a zonal mean and a deviation from the zonal mean, then neglects the eddy-eddy nonlinearities within the eddy equation [15]. Put another way, the QL approximation neglects all triad interactions except those containing a ZF mode. The QL system is given by

Motivation for using the quasilinear approximation
Srinivasan and Young [15] have shown that the QL system exhibits many of the same basic zonal jet features as the full NL system, including the formation of stable jets and merging jets. With periodic boundary conditions, the equation of motion enjoys translational symmetry in both the x and y directions. As a parameter is varied, the simulations suggest a spontaneous breaking of statistical homogeneity in the y direction. At large μ (small γ), the NL system in figure 1(a) and the QL system in figure 1(d) do not exhibit steady ZFs, so the behavior is statistically homogeneous. At small μ (large γ), both the NL system in figure 1(b) and the QL system in figure 1(e) do exhibit steady ZFs, implying a breaking of statistical homogeneity in the y direction. We also observe that, in both the NL and QL systems, simulations that differ only in initial conditions and realizations of the random forcing can display different numbers of jets (figures 1(b), (c), (e), (f)). In addition to these features, both NL and QL exhibit merging jets, evident in figures 1(c), (f). Our motivation in adopting the QL approximation is not because we believe it to be quantitatively correct, but rather because the QL system apparently retains the necessary ingredients that lead to the rich behavior of ZF formation. The QL system may provide insight into the more realistic models, and the advantage, of course, is that the QL system is far more tractable. The phenomena described above will all be explained analytically within the QL approximation.
Separate from our motivations for using the QL approximation, Bouchet et al have argued that in the regime of large γ the flow becomes predominantly zonal and the QL approximation becomes rigorously valid [37]. Our present study examines the regime in parameter space where γ is not asymptotically large, for it is in this regime where ZFs are born at low amplitudes from turbulence. We study values of γ ≈ 7, which are close to the transition to the 'zonostrophic regime' at ≈ β R 2 [36], but we caution against direct quantitative comparison because that investigation did not use the QL approximation whereas we do.
The QL approximation has also been used by Herring in the study of thermal convection [39], where the only NL interaction retained was between a horizontally-averaged temperature and the fluctuating temperature and velocity; NL interactions between the fluctuating quantities were discarded. At large Rayleigh number, this approximation was able to reproduce some of the qualitative features observed in experiments.

CE2 equation
The eddy quantity ′ w fluctuates rapidly in space and time. Averaging over these turbulent fluctuations enables one to work with smoothly varying functions. Such statistical approaches provide one path to gaining physical insight. Sometimes statistical turbulence theories strive for quantitative accuracy, which requires rather complicated methods [40], but we eschew those methods here because they are not required for investigation of the QL system.
We consider an average of the QL system (4). Derivations can be found in [22,26] though we follow [15] because there are advantages to its formulation. The full derivation can be found there, but we give a brief overview of the procedure. One defines the two-point, one-time correlation function of vorticity using a zonal average as x L x where L x is some averaging length, the integration is over the sum coordinate = + x x x ( ) 1 2 1 2 , and the difference coordinate = − x x x 1 2 is held fixed. The correlation function Ψ of streamfunction can be defined similarly. One finds an evolution equation for W by taking a time derivative of (5), substituting the expression for ∂w t from (4a), and performing the average. Under an ergodic assumption, the zonal average is equivalent to a statistical ensemble average, and the stochastic forcing can be averaged to a deterministic quantity. Then one performs a linear coordinate transformation to the sum and difference variables = + y y y ( ) is the covariance of the external forcing, and D h is the hyperviscosity operator, given by It can be shown from the definitions that W and Ψ are related by y y y y y y The use of the sum and difference coordinates x y y , , allows the structure of the theory and especially of the bifurcation to be more easily understood than in the original coordinates. In the new coordinates x and y represent two-point separations and y represents the two-point average position. If the turbulence were homogeneous, there would be no dependence on y .
The only assumption necessary for CE2 to be an exact description of the QL model is ergodicity in the zonal (x) direction, such that a zonal average is equivalent to an ensemble average. No other assumptions are required because the QL model neglects the NL eddy-eddy term that would give rise to a closure problem. Alternatively, instead of the QL-based derivation, CE2 can be regarded as a truncated statistical closure of the NL model [21,22,[26][27][28]. However, for present purposes we prefer the former interpretation.
CE2, like the QL system, exhibits merging jets [21]. Since CE2 is deterministic, if the system approaches a stable steady state then merging and branching of jets can only occur transiently. Once the stable equilibrium is reached, the system is stuck there and no more dynamical behavior can occur. Non-steady attractors have been seen within CE2 as well, although these do not appear common [22,41]. Additionally, if the QL system is not fully ergodic, then CE2 is not an exact description of it and dynamical behavior like merging or branching can persist even in a statistically steady state [22,37]. Though ergodicity is often a useful idealization, lack of complete ergodicity is to be expected in any physical system.
The CE2 equations exhibit important symmetries of translation and reflection, given by δ → + y y y a , In other words, if W x y y t U y t { ( , , ), ( , )} is a solution, then so are where δy is some constant translational shift. The symmetry (9d), dubbed the exchange symmetry, is always obeyed by the correlation function such that

Zonostrophic instability
In this section we review ZI, for which substantial understanding has been recently obtained [15,24]. A homogeneous, steady-state solution of the CE2 equations always exists, arising from a simple balance between forcing and dissipation. This solution is This equilibrium is linearly stable in a certain regime of parameters. As a control parameter ρ is varied, this homogeneous state becomes zonostrophically unstable [15,21]. Physically, ZI occurs when dissipation is overcome by the mutually reinforcing processes of eddy tilting by ZFs and production of Reynolds stress forces by tilted eddies. The instability eigenmode consists of perturbations spatially periodic in y with zero real frequency [15], so ZI arises as a Type I s instability [30] of homogeneous turbulence.
To calculate the dispersion relation corresponding to ZI, one considers perturbations about the equilibrium (10a), (10b). Because the equilibrium is independent of y and t, the y and t dependence of the perturbations can be Fourier transformed. The fields are written as t qy i where q is the ZF wave number. It is convenient to Fourier transform in both x and y as well.
One is able to obtain a single NL equation for the eigenvalue λ by solving for δW k k ( , ) x y in terms of δU, then substituting into the equation for δU [15]. The dispersion relation is 1 , x Some algebraic manipulation shows that Λ Λ = − + − [15]. One can calculate λ q ( ) numerically, plots of which can be seen in figure 4 of [15]. The equilibrium is unstable to ZI if λ > q ( ) 0 anywhere. The equilibrium transitions from being ZI stable to unstable at a critical value ρ c called the instability threshold. At ρ c there is a single marginally stable mode at wave number q c (and at −q c ). Beyond ρ c there is a band of unstable modes.

Beyond zonostrophic instability: bifurcation and the amplitude equation
The existence of ZI indicates that a homogeneous equilibrium without zonal flow is unstable. Perturbations to this equilibrium grow exponentially, with wave number dependencies and growth rates that can be calculated. However, the ZI calculation alone does not predict how the system saturates.
To understand the behavior in the regime of nonlinearly interacting eddies and ZFs, we turn to a bifurcation analysis. Near the instability threshold, the distance from the threshold serves as a small parameter to facilitate analytic progress. It has been demonstrated numerically that the bifurcation is supercritical [21], which we confirm with our analytical calculations in the appendix. Thus only lowest-order terms in the bifurcation analysis are needed to provide saturation of the instability.
The bifurcation analysis follows a standard procedure and involves a multiscale perturbation expansion about the threshold [30,31]. A normalized parameter is defined as If we denote u as the state vector relative to the homogeneous equilibrium, i.e., , then the expansion proceeds as 1/2 1 2

3/2 3
At first order, one q y i c that undergoes bifurcation, and A is its amplitude. The amplitude is an envelope that slowly varies in space and time. The slow variation represents the effect of the infinity of wave numbers nearby q c that also go unstable when ϵ > 0. The goal is to determine A, as then u 1 will be fully specified. Here, one determines a PDE for A as a solvability condition at third order in the perturbation expansion. One eventually finds where the c i are order unity, real, positive constants to be calculated. Equation (16) is referred to as the amplitude equation, or sometimes as the real Ginzburg-Landau equation.
It turns out that, in order to understand the qualitative behavior of A, one does not need to carry out the calculation of the c i explicitly [31]. This is because the translation and reflection symmetries (9a) constrain the lowest-order PDE for A to consist generically of the form in (16). The qualitative behavior of (16) is independent of its parameters since they can be transformed to unity by simple rescaling.
But even if the qualitative behavior is understood, it is still worthwhile to carry out the calculation of the coefficients c i . First, computing these and verifying the results numerically provides a concrete check on our overall understanding. Second, the perturbation solution is convenient for certain numerical methods where it is useful to start with a good approximation to the true solution. In the appendix, we derive (16) and obtain expressions for the c i . This calculation has also been carried out independently [42]. To verify our results, we compare the analytic growth rate found from (16) with that from the exact dispersion relation (12). Similarly, the analytic ZF amplitude found from (16) is compared with that from solving the ideal states numerically as in section 6. The results are shown in figure 2 and are in excellent agreement.
The amplitude equation (16) is well understood [30][31][32], and much of its qualitative behavior is seen generically in pattern-forming systems. First, with all the parameters c i and ϵ set to unity, a steady-state solution exists for any wave number within the continuous band solutions with < k 2 1 3 are stable; those with > k 2 1 3 suffer the Eckhaus instability [31]. This is demonstrated in figure 3, where an unstable solution that has been slightly perturbed undergoes merging behavior until a stable wave number is reached.
The CE2 system is described by this bifurcation, so near the threshold, and more generally, CE2 exhibits solutions existing with a range of ZF wave numbers having a certain stability region. In the following sections we calculate the equilibria and stability of CE2.

Calculation of ideal states
We proceed to find the steady-state solutions of (6). In the context of an infinite domain with no boundaries, these solutions are referred to as ideal states. Let q denote the fundamental ZF wave number of an ideal state. For a given q, we solve the time-independent form of (6) directly. Since we are able to select q and determine the ZF wavelength directly, this method differs from the finite-spatial-domain techniques of [21,22,28]. Within a finite spatial domain, the wave numbers takes discrete values. The dominant ZF mode is not preselected and is typically not the lowest mode because the system evolves self-consistently to find a solution. Our infinitedomain technique, which allows the selection of the ZF wavelength, is advantageous for understanding the global NL dynamics. We can solve for both stable and unstable solutions as we continuously vary q and therefore can easily determine stability boundaries.
There is an important side effect to our choice of the dominant ZF wave number. If the dominant ZF wave number does not occupy the lowest mode, subharmonics of the dominant mode can be excited. Our method involves setting the lowest ZF mode to be the dominant one. This requires fewer resolved modes. It also excludes subharmonics (although that is not a limitation in principle).
An ideal state is represented as a Fourier-Galerkin series with coefficients to be determined [31,33,34]. We expand as follows:   Reprinted with permission from [20]. Copyright 2013, AIP Publishing LLC.
While the periodicity in y is desired, the correlation function should decay in x and y; periodicity in x and y is a consequence of using the convenient Fourier basis. Thus, a and b, unlike q, are numerical parameters. They represent the spectral resolution of the correlation function and should be small enough to obtain an accurate solution.
Because the CE2 equations have translational symmetry in y , there is an infinite number of solutions, all equivalent, corresponding to displacements along y . In order to obtain a wellposed numerical problem, one must restrict the set of solutions. To this end, we again look to the symmetries (9a). The CE2 symmetries allow us to seek a solution for which In other words, we choose the origin of y such that the reflection symmetries hold for the solution itself. We find that such solutions do exist. It turns out that this restriction does not uniquely specify the solution, as shifting a solution by a half wavelength δ π = y q / yields a distinct but equivalent solution. Still, this restriction is sufficient to make the problem wellposed numerically. The above constraints, along with the conditions that U y ( ) and | W x y y ( , ) are real, force U p to be real, and p p Furthermore we take = U 0 0 , as that would merely represent a static uniform velocity. With these constraints, the number of independent, real degrees of freedom is . Aside from the previous statements, there is still no guarantee that there is a physically unique solution. Indeed, in the zonostrophically unstable regime, once the above ansatz with a specific q has been substituted there are at least two solutions: the equilibrium with ZFs, and the unstable homogeneous solution without ZFs. In some instances we also find other unstable solutions, which may be unphysical artifacts of the numerical discretization. We obtain a system of NL algebraic equations for the coefficients W mnp and U p by substituting the Galerkin series into (6) and projecting onto the basis functions. To demonstrate the projection for (6a), let ϕ = e e e mnp max nby pqy i i i . We project (6a) onto ϕ rst by operating with . The other terms of (6a), as well as (6b), are handled similarly. In total, we generate as many equations as there are coefficients.
The system of NL algebraic equations is solved with a Newtonʼs method [43]. In practice we separate W mnp into real and imaginary parts. The Jacobian matrix is sparse and is easy to specify analytically. We note that because the ZF equation is linear, it is possible to eliminate the ZF degrees of freedom analytically. This is avoided, however, because the reduction of only P degrees of freedom is negligible and this step incurs the major disadvantage of making the Jacobian no longer sparse. A Newtonʼs method requires a good initial guess. An accurate initial guess near the instability threshold is provided by the bifurcation calculation (see the appendix).
To find other solutions we employ simple numerical continuation, where the solution at one value of a parameter is used as the initial guess for the solution at the next value of the parameter. Figure 4 shows the ZF amplitude coefficients U p as functions of q at μ = 0.21 and μ = 0.19. Near the instability threshold, ideal states exist at all q for which the homogeneous equilibrium is zonostrophically unstable (between the two lines labeled N in figure 4(a)). Farther from threshold, there is a region of q where the ideal state solution disappears (between the lines N and D in figure 4(b); see also figure 5).
The computational method as described above works very well near the threshold μ = 0.237 c (γ = 6.02 c ). However, far from the threshold, for μ < 0.12 (γ ≳ 14.2), the numerical method breaks down. This appears to be related to the existence of multiple solutions at a given parameter value, of which some are unphysical or unstable. Far from threshold the Newtonʼs method seems to inevitably get stuck on one of these undesirable solutions. However, the timeevolving simulations previously mentioned do not have these problems because the CE2 equations are statistically realizable and will only approach physical, stable solutions.

Stability of the ideal states
We now turn to calculate the stability of the ideal states. Ideal state stability, which concerns the inhomogeneous equilibria, is distinct from ZI, which is a property of the homogeneous equilibriun. Both types of instabilities can be described within the CE2 formalism. We consider perturbations δ | W x y y t ( , , ) and δU y t ( , ) about an inhomogeneous equilibrium W U { , }. The linearized CE2 equations are   Since the underlying equilibrium is periodic in y , the perturbations can be expanded as a Bloch state [31,33]: t Qy mnp mnp max nby pqy where Q is the Bloch wave number and can be taken to lie within the first Brillouin zone − < ⩽ q Q q 1 2 1 2 . We do not use a Q x or Q y associated with the periodicities in x and y because, as previously mentioned, those periodicities are artificial. A separate symmetry argument also suggests taking Q x and Q y to be zero. Specifically, the correlation function exchange symmetry, which we continue to enforce in the perturbations, requires δ δ | = − − | W x y y W x y y ( , ) ( , ). This forces Q x to be either zero or a 1 2 and Q y to be either zero or b  1 complex coefficients from the δU p . Equation (21a), (21b) is projected onto the basis functions in the same way as in the ideal state calculation. The projection results in a linear system at each Q for the coefficients δW mnp and δU p ; this determines an eigenvalue problem for σ. The equilibrium is unstable if for any Q there are any eigenvalues with σ > Re 0. The eigenvalues are computed using Arnoldi iteration. The resulting equations contain ZI as a special case, for which the equilibrium is the homogeneous one and Q takes on the role of the wave number q.
It is possible to show the following two symmetries regarding the eigenvalue, which result from the symmetry of the ideal state equilibrium. First, for arbitrary Q, suppose δ δ ′ ( ) is an eigenvector with eigenvalue σ. Then the vector δ δ ′ ( ) (2) is also an eigenvector, with eigenvalue σ * , and This guarantees that every complex eigenvalue comes in a conjugate pair. Second, suppose at some Q that δ δ ′ ( ) Thus, for determining stability one needs to check only ⩽ ⩽ Q q 0 1 2 because the eigenvalues for negative Q are symmetric.
The stability diagram is shown in figure 5. To vary γ, we change μ and hold other parameters fixed at their previous values. The stable ideal states exist inside of the marginal stability curves marked E, L 1 , and R 1 , which represent different instabilities. The Eckhaus instability (E) is a long-wavelength universal instability, present even in the amplitude equation (16). The L 1 and R 1 curves represent novel short-wavelength instabilities. The zonal jets are spontaneously generated by ZI for γ γ > = 6.02 c . For γ < < 6.02 14.2, the stability curve is consistent with the dominant ZF wave number observed in QL simulations. For γ > 14.2, we could not calculate the stability diagram with this approach due to the aforementioned numerical issues of finding the steady state.
This stability calculation differs from the Rayleigh-Kuo criterion [21,22]. That criterion, which states that a necessary condition for inviscid barotropic instability is that β − ′′ U y ( ) changes sign, is derived assuming laminar flow with a fixed mean zonal velocity profile U(y) [44,45]. Such an assumption is relevant to the case where external heating and frictional forces directly determine the laminar flow profile. In practice, simulations of shallow or two-dimensional systems do typically satisfy the Rayleigh-Kuo sufficient condition for stability [2]. Some real flows are observed to have β − ″ U y ( ) change sign, such as Jupiterʼs zonal wind. (However, when 3D effects are taken into account, an appropriately generalized stability criterion may be satisfied [46][47][48][49].) Even in 2D, however, the Kuo criterion is not necessarily the only important governing principle. In the present case, the ZF is directly driven by turbulence. Fundamentally, there is no reason a priori to neglect the turbulence in a stability calculation. Especially in the regime we have focused on, just above threshold where the ZF is weak relative to the turbulence, a stability criterion that is based on laminar flow without turbulence is inapplicable.
One consequence of retaining the turbulence in an equilibrium is that in a linear stability analysis the ZF profile is allowed to change, as in (21b) [21]. These ZF perturbations are absent in the Rayleigh-Kuo criterion. In fact, near threshold the instability eigenvectors are found numerically to have a nontrivial component in the ZF perturbation. This numerical result is also evident from analytic consideration of the Eckhaus instability in the amplitude equation (16).
It is also important to note that because the Rayleigh-Kuo criterion is based on a linearization around laminar flow, the turbulent self-interaction is ignored. This interaction is exactly what is neglected by the QL approximation. Therefore, the Kuo criterion may be found just as well from the QL approximation. In other words, the QL approximation cannot be blamed on its face for the inapplicability of the Rayleigh-Kuo criterion. Figure 6(a) shows, at μ = 0.15 (γ = 10.7) and q = 0.6, a comparison of β and ″ U . For the same value of μ, figure 6(b) shows the minimum value of β − ″ U at each q along with the actual stability boundaries that have been computed directly. At large q, the curve approaches β because the ZF amplitude goes to zero. For this value of μ, at no q does β − ″ U change sign.
From the above theoretical considerations and numerical evidence, we conclude that in the QL or CE2 system near threshold, the operative instabilities determining the ZF length scale are independent of the Kuo criterion. This conclusion has been offered previously [21]; we have provided further analysis to support it.

Wave number selection
As evident from figure 5, we are presented with the theoretical quandary of having a wide range of allowed, stable solutions and yet a narrow preferred region where QL realizations tend to appear. This is common to pattern-forming systems, and the problem of wave number selection is difficult [31]. The Rhines wave number One might naturally inquire as to whether the preferred mode is the fastest growing mode in the ZI about the homogeneous equilibrium. This does not appear to be the case away from the threshold at larger γ [15,21], as seen in figure 5. There is, however, a plausible scenario that emerges which may explain the merging of jets often observed in the beginning stages of simulations, especially those which initialize everything at low amplitudes. At large γ, it appears that the fastest growing mode may be to the right of the stability region. In a simulation, the turbulence quickly comes to a quasi-equilibrium on a short time scale and begins to drive the ZF. The growing ZF mode cannot stably saturate, for its wavelength is too small to coexist with the turbulence. As the system evolves through the subsequent instability to drive the jets toward larger wavelength, a Hovmöller visualization displays merging jets.
Another possibility is that some kind of variational principle applies. The amplitude equation, by which CE2 is governed near threshold, is a gradient system. The ideal states of varying wave number q have varying values of the effective free energy. However, the minimum of the effective free energy is not necessarily dynamically preferred [50]. In any case, away from threshold CE2 is not a gradient system and there is no rigorous theoretical basis for expecting variational behavior to occur. From a different perspective, variational principles for certain 2D turbulent systems have long been discussed theoretically. Some of these principles are based on the nonlinearly conserved quadratic quantities, the energy and the enstrophy. For instance, in freely decaying turbulence where viscosity provides the dissipation, the enstrophy is expected to decay more quickly than the energy. One might expect the decaying turbulence to reach a state of minimum enstrophy subject to the constraint of constant energy. Other principles exist based on minimum dissipation or maximum entropy or entropy production [51]. Although these principles do not directly apply to the damped, driven CE2 system, they at least motivate a numerical exploration to try to discover any correlation between the preferred wave numbers and other properties.
As a simple starting point for our exploration, we examine the energy and enstrophy of the ideal states. Plots of the energy and enstrophy, for both the total and just the eddies, are shown in figure 7 for μ = 0.15. For each quantity a distinct minimum is present. We find at each μ the minimum of all four quantities; the resulting curves are shown in figure 5. While the minima of the total energy and total enstrophy are consistent with the QL realizations, there is no clear indication that either is especially preferred. On the other hand, the accessible regime investigated here is not too far from threshold, so this is not in the asymptotic regime of large γ.

Discussion
The averaging procedure to obtain the CE2 equations from the QL equations merits further discussion. Under appropriate assumptions, which always include some kind of ergodicity assumption, multiple choices of average will lead to the same final equations. For instance, zonal [15], short-time [24], and coarse-graining [52] averages have been discussed. The ergodicity assumption allows one to transform the average over the random forcing into a deterministic quantity. One can also discuss things in terms of an ensemble average, in which case an assumption of statistical homogeneity in the zonal (x) direction is made, but inhomogeneity is allowed in the nonzonal (y) direction. In this case, ergodicity is not required in order to derive the CE2 equations, but it becomes necessary if one wants to interpret the solutions of the equations as having anything to do with the behavior of an individual realization.
When using the ensemble average, Kraichnan pointed out in the context of thermal convection that the definition of the statistical ensemble is somewhat subtle for the situation of spontaneous symmetry breaking [53]. Because of the translational symmetry, the zonal jets have no preferred location and are presumably equally likely to form with any particular phase. One choice of the statistical ensemble encompasses all possible realizations consistent with the prescribed parameters, in which case the ensemble itself is statistically homogeneous and any ensemble-averaged quantity must be homogeneous also. Therefore the average yields zero mean ZF (and then the ZF must be described as a fluctuation), despite the fact that each individual realization has a nonzero ZF. This was the procedure followed in [19]. Another possibility is that the ensemble might consist only of the realizations for which the zonal jets have a particular phase. The latter interpretation is the one that yields the CE2 equations identical to those obtained by zonal averaging. With the former ensemble, the ergodic assumption is invalid, since an individual realization is no longer mixing throughout the full set of realizations of this ensemble. This is consistent with the fact that the ensemble-averaged behavior is not equivalent to the behavior of an individual realization.
Quantitatively, CE2 is not expected to be accurate in all situations (although it may be in some cases [52]). Not only are the eddy self-nonlinearities entirely ignored, which neglects the traditional cascades, but CE2 involves only one-time correlation functions rather than the more general two-time functions. The lack of time-history information means that most of the effects of wave propagation are discarded [54]. To achieve greater quantitative accuracy, the effect of eddy self-nonlinearities must be retained in some way. One way this could be done is by incorporating the spectrum from NL simulations [25]. Another way is to go beyond CE2 and use third-order cumulants in a CE3 framework [28]. Yet another alternative is to use inhomogeneous second-order statistical closures such as the direct interaction approximation (DIA) [53]. Progress developing inhomogeneous closures has been limited. A simpler DIA variant called the quasi-diagonal DIA exists [55][56][57], but it approximates the interaction between the mean field and the fluctuation. That approximation would affect the stability properties of the ZF in ways currently unknown. Additionally, an inhomogeneous Markovianized closure exists in the test-field model [58], but it is not statistically realizable in the presence of waves [59,60]. A homogeneous realizable test-field model exists [60], but as of yet there is no version that is both realizable and inhomogeneous. More work along these lines needs to be done.
There is an important situation for which it is crucial to have the eddy self-nonlinearities retained in a statistical sense, at least to describe the bifurcation from homogeneous to inhomogeneous turbulence: a system with linear instability. The oceans are baroclinically unstable and linear instabilities in plasmas are common. Some plasma models have clearly demonstrated the symmetry-breaking bifurcation of ZF generation [61]. In order to describe this transition, a model must have an equilibrium of homogeneous turbulence to begin with. But in a QL or CE2 description, if no ZFs are present then there are no NL interactions, and it is impossible for a linear instability to saturate. With linear instability present, a QL description has no homogeneous equilibrium. We also point out that it is not difficult to modify the phenomenological model in section 2 to include an eddy self-nonlinearity.
It is important to understand to how our analysis, which has been carried out on an infinite β plane, pertains to a more geophysically realistic setting such as the surface of a rotating sphere. In that case, because of the variation of the Coriolis parameter, the turbulence is always inhomogeneous with respect to latitude and there is no latitudinal symmetry to be broken. However, there can still be spontaneous formation of a mean field, in which case the ZF continues to play the role of an order parameter. As some control parameter is varied, a statistical steady state without ZFs can go unstable in a bifurcation at which ZFs are generated. This scenario appears to have been observed in numerical simulations when the rotation rate is increased from zero [62]. Recently CE2 on the rotating sphere has been studied and ZFs have been observed within that framework [26,27]. One avenue for future investigation is to use CE2 to study ZI on the sphere. Qualitative insight could be gained into the structure of the unstable eigenfunction, including the direction of the equatorial jet.
In conclusion, we have examined in detail the transition from statistically homogeneous turbulence without ZFs to statistically inhomogeneous turbulence with steady ZFs. Our study has been carried out within the QL approximation using the CE2 statistical framework, which appears to contain much of the essential physics of ZF generation yet is amenable to analysis. We believe, and numerical evidence supports, that many of the features identified and examined here, including the spontaneous symmetry breaking of statistical homogeneity and nonuniqueness of the ZF wave number, are present in the full system as well as the QL system. We have shown analytically within CE2 that multiple steady-state solutions exist with varying ZF wavelengths, and that only those wavelengths within a certain subband are linearly stable. Unstable wavelengths must evolve to reach a stable wavelength, and this process manifests as merging jets. The simplest realization of the phenomenon of merging jets appears in the wellknown amplitude equation (16), which governs the CE2 system near threshold.
Here we derive the amplitude equation (16) directly from the CE2 equations (6) and verify the results numerically. First, we review the procedure for the perturbation expansion [31]. Then we fill in the algebraic details. We limit ourselves to quadratic nonlinearity. Let ϕ be an abstract vector, Λ be a linear operator, N be a symmetric, bilinear operator, and F be external forcing. Any of Λ, N, and F may depend explicitly on the small parameter ϵ. states that u 1 is an eigenvector with a zero eigenvalue. Then u 1 can be a linear combination of null eigenvectors with a to-be-determined amplitude. The reality condition on u restricts the form to be where ∼ r e q y i c (and its complex conjugate) are the right null eigenvectors. Here, q c is the critical wave number that first goes unstable as the parameter ϵ crosses zero. Given an inner product · · ( , ), then associated with the right null eigenvector is a left null eigenvector l of L 0 , such that = l L u ( , ) 0 0 for any u. The y dependence of l will also be e q y i c . The amplitude A will be determined by nonlinearities occurring at higher order.
At ϵ O ( ), we first note that = L u 0 1 1 automatically. This is because q c is marginally stable at the instability threshold: given a dispersion relation λ ϵ q ( , ) as a function of wavenumber q and control parameter ϵ, then both λ or ± e q y 2i c due to the quadratic nonlinearity, this solvability condition is always satisfied. Thus, given that a solution exists, one may write u 2 as a linear combination of homogeneous and particular solutions: We have introduced another unknown parameter A 2 , but we will not need it in order to solve for A. At for the same reason that = L u 0 1 1 . Upon writing the solvability condition, one finds that several terms vanish, leaving = + l L u l N u u 0 ( , ) ( , 2 ( , )). (A.10) This is the desired partial differential equation which determines the amplitude A. Observe that one never explicitly needs L 1 or N 1 . We now apply this procedure to (6). For simplicity, we set the viscosity to zero, take infinite deformation radius, and cross the instability threshold by varying the strength of the forcing (rather than by varying the friction as in the main text); modification for other scenarios is obvious. Let the forcing be given by ϵ = + F x y F x y ( , ) ( The symmetrized operator is then given by (A.14) un un We now introduce the slow space and time scales. One subtlety that was not mentioned in the general procedure described above is that the u i may need to be expanded in ϵ. This occurs for two reasons. First, the ± U y y ( ) 1 2 terms lead to Second, the eigenvector r itself contains the differential operator ∂ y , which must be expanded in the multiple-scales procedure. It is extremely convenient to introduce these expansions at the outset so as to keep all of the ϵ expansion in a single place. This procedure is even more motivated when we absorb these extra terms into L 1 and L 2 , for these terms are necessary in order to satisfy = L u 0 . To introduce our convenient shortcut, first recall that since N 1 is never needed, we only need to perform this within L. Then, for the places where ± U y y ( ) Then, letting ∂ → q i y (which we can do because we will only need to perform this expansion on a term e q y i c and not other harmonics), we see that To implement this, one can set, in L, where the "∂ q " means that the derivative acts only on r(q), not on the e qy i part of u 1 . We need only make this replacement in W, not U, because the U component of the right null eigenvector does not contain any derivatives ∂ y . One can verify that this shortcut gives the same results as if one proceeded more straightforwardly.
The problem is most conveniently expressed in terms of the Fourier transform of the difference variables x y , . We use the convention U y x y x y x y Here,Û is the Fourier transform of U. Our expressions for U y ( ) will always consist of periodic exponentials, and soÛ contains delta functions and the convolution integral can be immediately performed. We also have defined The q dependence of r W and l W is now suppressed except for where it matters in (A.45); they should be evaluated at = q q c .
At ϵ O ( ), we need to solve the particular solution of (A.9). Take an ansatz    (16). The coefficients c i involve integrals over the forcing spectrum, which is here presented in a form where the wave numbers are shifted, i.e., contain terms like − k q y 1 2 . It is also possible to shift the integration variable so all integrals contain just the unshifted forcing F k k ( , ) x y 0 , after which the q derivatives in c 2 can be explicitly computed [42]. It is also possible to obtain c 0 , c 1 , and c 2 , which govern the linear behavior, via the alternate and much simpler route of using the analytic dispersion relation (12). The dispersion relation can be put into the form λ ϵ = D q ( , , ) 0. The conditions of the instability threshold require = D q (0, 0, ) 0