Magnetic Suppression of Zonal Flows on a Beta Plane

The system of Equation (1) supports two basic waves, the fast and slow magneto-Rossby waves, which are mixtures of the Rossby wave and the shear Alfvén wave. To derive the dispersion relations of the magneto-Rossby waves, we linearize the unforced equations of motion about $(\zeta ,A)=(0,0)$ and substitute perturbations of the form ${e}^{i{\boldsymbol{k}}\cdot {\boldsymbol{x}}-i\omega t}$. We obtain the dispersion relation

$\begin{eqnarray}{\omega }_{f,s} & = & \displaystyle \frac{{\omega }_{{\rm{R}}}}{2}\left(1-i\displaystyle \frac{(\nu +\eta ){k}^{2}}{{\omega }_{{\rm{R}}}}\right.\\ & & \left.\pm \sqrt{{\left[1-i\displaystyle \frac{(\nu -\eta ){k}^{2}}{{\omega }_{{\rm{R}}}}\right]}^{2}+\displaystyle \frac{4{\omega }_{{\rm{A}}}^{2}}{{\omega }_{{\rm{R}}}^{2}}}\right),\end{eqnarray} \tag{ 4 }$

where ${k}^{2}={k}_{x}^{2}+{k}_{y}^{2}$ and ${\omega }_{{\rm{R}}}\mathop{=}\limits^{\mathrm{def}}-\beta {k}_{x}/{k}^{2}$ and ${\omega }_{{\rm{A}}}\mathop{=}\limits^{\mathrm{def}}{k}_{x}{B}_{0}$ are the frequencies of the undamped Rossby and shear Alfvén waves, respectively. The fast wave ${\omega }_{f}$ takes the + sign, and the slow wave ${\omega }_{s}$ takes the − sign.

The eigenmodes can be obtained from the linearized magnetic equation as

$\begin{eqnarray}&&{\left[\begin{array}{l}\zeta \\ A\end{array}\right]}_{f,s}=\displaystyle \frac{1}{k\sqrt{| {\omega }_{f,s}+i\eta {k}^{2}{| }^{2}+{\omega }_{{\rm{A}}}^{2}}}\left[\begin{array}{c}{k}^{2}({\omega }_{f,s}+i\eta {k}^{2})\\ -{\omega }_{{\rm{A}}}\end{array}\right].\end{eqnarray} \tag{ 5 }$

For later convenience, the normalizing factor has been chosen such that the quantity ${k}^{2}(| \psi {| }^{2}+| A{| }^{2})={k}^{2}(| \zeta {| }^{2}/{k}^{4}+| A{| }^{2})$, which is equal to the mode energy (up to a factor of two), is unity.

We examine two limits to elucidate the physical nature of these waves. First, in the nondissipative limit where ν and η vanish, the frequencies are

$\begin{eqnarray}&&{\omega }_{f,s}=\displaystyle \frac{{\omega }_{{\rm{R}}}}{2}\left(1\pm \sqrt{1+\displaystyle \frac{4{\omega }_{{\rm{A}}}^{2}}{{\omega }_{{\rm{R}}}^{2}}}\right).\end{eqnarray} \tag{ 6 }$

For vanishing magnetic field the Rossby wave is recovered, while for strong magnetic field the shear Alfvén wave is recovered.

Second, in this paper we focus on the regime where in the length scales of interest, the Rossby wave is the fastest process, such that νk², $\eta {k}^{2}$, ${\omega }_{{\rm{A}}}\ll {\omega }_{{\rm{R}}}$. In this regime, Equation (4) reduces to

$\begin{eqnarray}&&{\omega }_{f}={\omega }_{{\rm{R}}}-i\nu {k}^{2},\end{eqnarray} \tag{ 7a }$

$\begin{eqnarray}&&{\omega }_{s}=-\displaystyle \frac{{\omega }_{{\rm{A}}}^{2}}{{\omega }_{{\rm{R}}}}-i\eta {k}^{2}.\end{eqnarray} \tag{ 7b }$

The fast wave is essentially the Rossby wave, while the slow wave involves both the magnetic field and beta effect. In this regime, the eigenmodes in Equation (5) simplify to

$\begin{eqnarray}&&{\left[\begin{array}{c}\zeta \\ A\end{array}\right]}_{f}=\displaystyle \frac{1}{k}\left[\begin{array}{c}{k}^{2}\\ -{\omega }_{{\rm{A}}}/{\omega }_{{\rm{R}}}\end{array}\right],\end{eqnarray} \tag{ 8a }$

$\begin{eqnarray}&&{\left[\begin{array}{c}\zeta \\ A\end{array}\right]}_{s}=\displaystyle \frac{1}{k}\left[\begin{array}{c}{k}^{2}{\omega }_{{\rm{A}}}/{\omega }_{{\rm{R}}}\\ 1\end{array}\right].\end{eqnarray} \tag{ 8b }$

In this regime, the fast wave is dominated by the vorticity component and the slow wave is dominated by the magnetic component.

A useful framework for addressing the dynamics of coherent flows embedded in and driven by turbulence involves studying the dynamics of the statistics of the flow fields (e.g., statistical moments). Rather than working directly with flow fields that rapidly vary in time and space, studying the behavior of dynamical equations for statistical quantities can provide qualitative insight of turbulence–mean flow interaction. However, forming statistically averaged equations of nonlinear systems inevitably runs into the closure problem, where an infinite hierarchy of moment equations is required to obtain a closed system. Thus, a turbulence closure is needed.

Here, we study the dynamics of the magnetized fluid in Equation (1) using the quasilinear second-order closure. This closure has proven useful in gaining analytic understanding and physical insight regarding coherent-structure formation in turbulent flows. In the quasilinear second-order closure, the eddy-mean flow interaction is accurately captured; indeed, this interaction is not approximated whatsoever. This particular closure comes (unfortunately) in the literature under two names: "S3T," which stands for Stochastic Structural Stability Theory (Farrell & Ioannou 2003) and "CE2," which stands for Cumulant Expansion at second order (Marston et al. 2008). Hereafter we refer to this closure as CE2.

We consider a decomposition of the flow fields into a coherent and an incoherent component. Here, we identify the coherent component with the zonal mean (denoted by over bar) and the incoherent component, or eddies, with the fluctuations about the zonal mean (denoted by prime), e.g.,

$\begin{eqnarray}&&\overline{\psi }(y,t)\mathop{=}\limits^{\mathrm{def}}\displaystyle \frac{1}{{L}_{x}}{\int }_{0}^{{L}_{x}}{dx}\,\psi ({\boldsymbol{x}},t),\end{eqnarray} \tag{ 9a }$

$\begin{eqnarray}&&\psi ^{\prime} ({\boldsymbol{x}},t)\mathop{=}\limits^{\mathrm{def}}\psi ({\boldsymbol{x}},t)-\overline{\psi }(y,t).\end{eqnarray} \tag{ 9b }$

The quasilinear approximation consists of neglecting the eddy–eddy nonlinearity in the eddy evolution equations while keeping the mean flow dynamics intact. Thus, from Equation (1), we obtain the quasilinear equations

$\begin{eqnarray}&&{\partial }_{t}\bar{u}=\overline{v^{\prime} \,\zeta ^{\prime} }-\overline{({\partial }_{x}A^{\prime} ){{\rm{\nabla }}}^{2}A^{\prime} }+\nu {\partial }_{y}^{2}\bar{u},\end{eqnarray} \tag{ 10a }$

$\begin{eqnarray}&&{\partial }_{t}\bar{A}=-{\partial }_{y}(\overline{v^{\prime} A^{\prime} })+\eta {\partial }_{y}^{2}\bar{A},\end{eqnarray} \tag{ 10b }$

$\begin{eqnarray}&&{\partial }_{t}\zeta ^{\prime} +\bar{u}{\partial }_{x}\zeta ^{\prime} +(\beta -{\partial }_{y}^{2}\bar{u})v^{\prime} \\ &&=\,-({B}_{0}+{\partial }_{y}\bar{A}){\partial }_{x}{{\rm{\nabla }}}^{2}A^{\prime} +({\partial }_{y}^{3}\bar{A})({\partial }_{x}A^{\prime} )+\nu {{\rm{\nabla }}}^{2}\zeta ^{\prime} +\xi ,\end{eqnarray} \tag{ 10c }$

$\begin{eqnarray}&&{\partial }_{t}A^{\prime} +\bar{u}{\partial }_{x}A^{\prime} =-({B}_{0}+{\partial }_{y}\bar{A})v^{\prime} +\eta {{\rm{\nabla }}}^{2}A^{\prime} .\end{eqnarray} \tag{ 10d }$

From the quasilinear equations above we can form the closed system for the evolution of the first and second flow cumulants. The first cumulants being the mean flow components,

$\begin{eqnarray}&&\bar{u}\ \mathrm{and}\ \bar{A},\end{eqnarray} \tag{ 11 }$

while the second cumulants are the same-time two-point eddy covariances:

$\begin{eqnarray}W & \mathop{=}\limits^{\mathrm{def}} & \overline{\zeta ^{\prime} ({{\boldsymbol{x}}}_{a},t)\zeta ^{\prime} ({{\boldsymbol{x}}}_{b},t)},\quad M\mathop{=}\limits^{\mathrm{def}}\overline{\zeta ^{\prime} ({{\boldsymbol{x}}}_{a},t)A^{\prime} ({{\boldsymbol{x}}}_{b},t)},\\ N & \mathop{=}\limits^{\mathrm{def}} & \overline{A^{\prime} ({{\boldsymbol{x}}}_{a},t)\zeta ^{\prime} ({{\boldsymbol{x}}}_{b},t)},\quad G\mathop{=}\limits^{\mathrm{def}}\overline{A^{\prime} ({{\boldsymbol{x}}}_{a},t)A^{\prime} ({{\boldsymbol{x}}}_{b},t)}.\end{eqnarray} \tag{ 12 }$

The stresses that appear in the mean flow Equations 10(a) and (b) are expressed in terms of the eddy covariances through

$\begin{eqnarray}&&\overline{v^{\prime} \zeta ^{\prime} }=\tfrac{1}{2}{[({\partial }_{{x}_{a}}{{\rm{\nabla }}}_{a}^{-2}+{\partial }_{{x}_{b}}{{\rm{\nabla }}}_{b}^{-2})W]}_{a=b},\end{eqnarray} \tag{ 13a }$

$\begin{eqnarray}&&\overline{({\partial }_{x}A^{\prime} ){{\rm{\nabla }}}^{2}A^{\prime} }=\tfrac{1}{2}{[({\partial }_{{x}_{a}}{{\rm{\nabla }}}_{b}^{2}+{\partial }_{{x}_{b}}{{\rm{\nabla }}}_{a}^{2})G]}_{a=b},\end{eqnarray} \tag{ 13b }$

$\begin{eqnarray}&&\overline{v^{\prime} A^{\prime} }=\tfrac{1}{2}{[{\partial }_{{x}_{a}}{{\rm{\nabla }}}_{a}^{-2}M+{\partial }_{{x}_{b}}{{\rm{\nabla }}}_{b}^{-2}N]}_{a=b},\end{eqnarray} \tag{ 13c }$

where the subscript a = b denotes that the function of ${{\boldsymbol{x}}}_{a}$ and ${{\boldsymbol{x}}}_{b}$ inside square brackets is transformed into a function of a single spatial coordinate by setting ${{\boldsymbol{x}}}_{a}={{\boldsymbol{x}}}_{b}={\boldsymbol{x}}$. Thus, the mean flow equations in the CE2 closure are exactly Equations 10(a) and (b) with the stresses given by Equation (13).

By manipulating Equations 10(c) and (d) and also using Equation 2(b) we obtain the evolution equations for the eddy covariances Equation (12):

$\begin{eqnarray}&&{\partial }_{t}W=({{ \mathcal L }}_{a}^{\zeta \zeta }+{{ \mathcal L }}_{b}^{\zeta \zeta })W\,+\,{{ \mathcal L }}_{a}^{\zeta A}N\,+{{ \mathcal L }}_{b}^{\zeta A}M+Q,\end{eqnarray} \tag{ 14a }$

$\begin{eqnarray}&&{\partial }_{t}M=({{ \mathcal L }}_{a}^{\zeta \zeta }+{{ \mathcal L }}_{b}^{{AA}})M+\,{{ \mathcal L }}_{a}^{\zeta A}G\,+{{ \mathcal L }}_{b}^{A\zeta }W,\end{eqnarray} \tag{ 14b }$

$\begin{eqnarray}&&{\partial }_{t}G=({{ \mathcal L }}_{a}^{{AA}}+{{ \mathcal L }}_{b}^{{AA}})G+{{ \mathcal L }}_{a}^{A\zeta }M+{{ \mathcal L }}_{b}^{A\zeta }N,\end{eqnarray} \tag{ 14c }$

where the ${ \mathcal L }$ operators depend on the mean flow fields, $\bar{u}$ and $\bar{A}$, and are given by

$\begin{eqnarray}&&{{ \mathcal L }}^{\zeta \zeta }\mathop{=}\limits^{\mathrm{def}}-\bar{u}\,{\partial }_{x}-(\beta -{\partial }_{y}^{2}\bar{u}){{\rm{\nabla }}}^{-2}{\partial }_{x}+\nu {{\rm{\nabla }}}^{2},\end{eqnarray} \tag{ 15a }$

$\begin{eqnarray}&&{{ \mathcal L }}^{\zeta A}\mathop{=}\limits^{\mathrm{def}}-({B}_{0}+{\partial }_{y}\bar{A}){{\rm{\nabla }}}^{2}{\partial }_{x}+({\partial }_{y}^{3}\bar{A}){\partial }_{x},\end{eqnarray} \tag{ 15b }$

$\begin{eqnarray}&&{{ \mathcal L }}^{A\zeta }\mathop{=}\limits^{\mathrm{def}}-({B}_{0}+{\partial }_{y}\bar{A}){{\rm{\nabla }}}^{-2}{\partial }_{x},\end{eqnarray} \tag{ 15c }$

$\begin{eqnarray}&&{{ \mathcal L }}^{{AA}}\mathop{=}\limits^{\mathrm{def}}-\bar{u}\,{\partial }_{x}+\eta {{\rm{\nabla }}}^{2}.\end{eqnarray} \tag{ 15d }$

In Equation (14), Q is the forcing covariance defined in Equation 2(b) and subscripts in the ${ \mathcal L }$ operators denote the variables on which the differential operators act and at which the mean flow fields are evaluated. We have assumed ergodicity to replace zonal averages over the random-forcing realizations with their ensemble averages. The evolution equation for mixed covariance N is redundant because of the symmetry

$\begin{eqnarray}&&M({{\boldsymbol{x}}}_{a},{{\boldsymbol{x}}}_{b},t)=N({{\boldsymbol{x}}}_{b},{{\boldsymbol{x}}}_{a},t).\end{eqnarray} \tag{ 16 }$

The evolution equation for N can be obtained from Equation 14(b) by exchanging $\zeta \leftrightarrow A$ in the superscripts of the ${ \mathcal L }$ operators together with exchanging $a\leftrightarrow b$ in the subscripts.

Note that only the quasilinear approximation in Equation (10) is enough produce the CE2 closure. Thus, a closure of the flow statistics at second order is exactly equivalent with the neglect of the eddy–eddy nonlinearity in the eddy dynamics.

The terms on the right-hand side of Equation 10(a) can be rewritten using integration by parts as

$\begin{eqnarray}&&\overline{v^{\prime} \,\zeta ^{\prime} }=-{\partial }_{y}\overline{u^{\prime} v^{\prime} },\end{eqnarray} \tag{ 17a }$

$\begin{eqnarray}&&\overline{({\partial }_{x}A^{\prime} ){{\rm{\nabla }}}^{2}A^{\prime} }=-{\partial }_{y}\overline{{B}_{x}^{{\prime} }{B}_{y}^{{\prime} }}.\end{eqnarray} \tag{ 17b }$

These identities allow the forces in Equation 10(a) to be written in the form of divergence-of-a-stress. Equation 17(a) is Taylor's identity that relates the vorticity flux with the Reynolds-stress divergence; Equation 17(b) is analogous to Equation 17(a) in providing an identity for the vorticity flux associated with the Maxwell stress. We will use either of the expressions in Equation 17(a) interchangeably and refer to them simply as the "Reynolds stress"; similarly we refer to either of the expressions in Equation 17(b) as the "Maxwell stress."

In summary, the CE2 equations consist of the evolution Equation (14) for the eddy covariances, and the evolution Equations 10(a)–(b) for the zonally averaged flow and magnetic potential (in which the stresses are given by Equation (13)).

Here, we demonstrate how the Kelvin–Orr shearing wave gives rise to the tendency for hydrodynamic fluctuations in the presence of a long-wavelength shear flow to transfer energy to the shear flow. The calculation was presented by Bakas & Ioannou (2013a), which we review here because we use the same techniques when we include a magnetic field in Section 3.2.

First, we consider the energetics of the mean flow. The zonally averaged momentum equation, ignoring magnetic fields and dissipation, is given by Equation 10(a): ${\partial }_{t}\bar{u}=-{\partial }_{y}\overline{u^{\prime} v^{\prime} }$. Multiplying by $\bar{u}$ and averaging over y, we obtain

$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{ZF}}}{{dt}}=\displaystyle \frac{1}{{L}_{y}}{\int }_{0}^{{L}_{y}}{dy}\,\overline{u^{\prime} v^{\prime} }{\partial }_{y}\overline{u},\end{eqnarray} \tag{ 18 }$

where ${E}_{\mathrm{ZF}}\mathop{=}\limits^{\mathrm{def}}\tfrac{1}{{L}_{y}}{\int }_{0}^{{L}_{y}}{dy}\,\tfrac{1}{2}{\overline{u}}^{2}$ is the spatially averaged energy density of the zonal flow and we have neglected boundary terms.

For the rest of this section, we consider the evolution of perturbation vorticity under the assumption of a fixed, linear shear flow $\bar{u}={Sy}$. Unlike a periodic $\overline{u}$, a linear flow appears incompatible with the neglect of boundary terms, a point which we will return to at the end of this section. The linearized evolution equation for vorticity is

$\begin{eqnarray}&&({\partial }_{t}+{Sy}{\partial }_{x})\zeta ^{\prime} =\nu {{\rm{\nabla }}}^{2}\zeta ^{\prime} .\end{eqnarray} \tag{ 19 }$

As we are interested in studying the emergence of zonal flows, we assume that $\bar{u}$ is very weak. This assumption implies that the shear S is very small, in a manner to be quantified later. We substitute an ansatz $\zeta ^{\prime} ({\boldsymbol{x}},t)=Z(t){e}^{i{\boldsymbol{k}}(t)\cdot {\boldsymbol{x}}}$. Requiring the coefficients of the terms linear in x and y to vanish, we see that ${{dk}}_{x}/{dt}=0$ and ${{dk}}_{y}/{dt}=-{{Sk}}_{x}$. Hence,

$\begin{eqnarray}&&{k}_{x}=\mathrm{constant},\end{eqnarray} \tag{ 20a }$

$\begin{eqnarray}&&{k}_{y}(t)={k}_{y0}-{{Sk}}_{x}t.\end{eqnarray} \tag{ 20b }$

The resulting equation for the amplitude Z can be solved, yielding

$\begin{eqnarray}&&\zeta ^{\prime} ={Z}_{0}{e}^{i[{k}_{x}x+{k}_{y}(t)y]}{e}^{-\nu {\int }_{0}^{t}d\tau k{(\tau )}^{2}},\end{eqnarray} \tag{ 21 }$

where $k{(t)}^{2}\mathop{=}\limits^{\mathrm{def}}{k}_{x}^{2}+{k}_{y}{(t)}^{2}$. Equation (21) describes the shearing wave. From Equation (21) we can compute $u^{\prime} \,={{ik}}_{y}(t)\zeta ^{\prime} /{k}^{2}(t)$ and $v^{\prime} =-{{ik}}_{x}\zeta ^{\prime} /{k}^{2}(t)$.

We next compute the net energy change of the mean flow due to a single wave that shears over and eventually dissipates. We combine the time-dependent shearing wave solution with our previous energetics calculations. We require the zonal average $\overline{u^{\prime} v^{\prime} }$, which is quadratic in wave fields. The wave fields are ultimately real, and accounting for their complex representation, we have

$\begin{eqnarray}&&\overline{u^{\prime} v^{\prime} }\to \displaystyle \frac{1}{2}\mathrm{Re}(u^{\prime} v{{\prime} }^{* })=-\displaystyle \frac{1}{2}\displaystyle \frac{{k}_{x}{k}_{y}(t)}{k{(t)}^{4}}| Z(t){| }^{2}.\end{eqnarray} \tag{ 22 }$

The change in ${E}_{\mathrm{ZF}}$ is obtained by integrating Equation (18) over the lifetime of the shearing wave:

$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{ZF}}={\int }_{0}^{\infty }{dt}\displaystyle \frac{{{dE}}_{\mathrm{ZF}}}{{dt}}=-\displaystyle \frac{1}{2}S{\int }_{0}^{\infty }{dt}\displaystyle \frac{{k}_{x}{k}_{y}(t)}{k{(t)}^{4}}| Z(t){| }^{2}.\end{eqnarray} \tag{ 23 }$

We used that ${\partial }_{y}\overline{u}=S$ (independent of y) and that the stress $\overline{u^{\prime} v^{\prime} }$ for an individual wave is also independent of y, and therefore the average over y does nothing.

Equation (23) shows that waves starting off with ${k}_{y}/{k}_{x}\gt 0$ (quadrants I and III in the ${\boldsymbol{k}}$ plane) will take energy from the zonal flow, while waves starting off with ${k}_{y}/{k}_{x}\lt 0$ (quadrants II and IV) will give energy to the zonal flow. The simplest form of the Kelvin–Orr shearing wave dynamics for growth of the shear flow arises from considering two waves at the same amplitude, with initial wavevectors $({k}_{x},{k}_{y0})$ and $(-{k}_{x},{k}_{y0})$. In isolation, one of the waves would grow in expense of the mean flow while the other would decay and give energy to the mean flow. The two waves must be considered together because the net leading order contribution to ${\rm{\Delta }}{E}_{\mathrm{ZF}}$ cancels out. We ignore interactions between waves, meaning that in the computation of the stress $\overline{u^{\prime} v^{\prime} }$, we ignore cross terms.

The total energy change of the zonal flow ${\rm{\Delta }}{E}_{\mathrm{ZF}}$ is the sum of that of the two waves individually, given by

$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{ZF}}=\displaystyle \frac{S}{2}{\int }_{0}^{\infty }{dt}\left[\displaystyle \frac{{k}_{x}{k}_{y+}(t)}{{k}_{+}{(t)}^{4}}| {Z}_{+}(t){| }^{2}-\displaystyle \frac{{k}_{x}{k}_{y-}(t)}{{k}_{-}{(t)}^{4}}| {Z}_{-}(t){| }^{2}\right].\end{eqnarray} \tag{ 24 }$

Here, a term with subscript ± stems from the wave with initial wavevector $(\mp {k}_{x},{k}_{y0})$, where ${k}_{y\pm }(t)\mathop{=}\limits^{\mathrm{def}}{k}_{y0}\pm {{Sk}}_{x}t$. We take the initial amplitudes ${Z}_{+}(0)={Z}_{-}(0)={Z}_{0}$. From Equation (21), we have ${| {Z}_{\pm }(t)| }^{2}={| {Z}_{0}| }^{2}{e}^{-2\nu {\int }_{0}^{t}d\tau {k}_{\pm }{(\tau )}^{2}}$. Assuming ${k}_{x}{St}/{k}_{y0}\,\ll \,1$, expanding to leading order in S, and dropping the 0 subscript on k_y0, we obtain

$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{ZF}}=\displaystyle \frac{{S}^{2}{k}_{x}^{2}| {Z}_{0}{| }^{2}}{4{\nu }^{2}{k}^{4}}\displaystyle \frac{{k}_{x}^{2}-5{k}_{y}^{2}}{{k}^{6}}.\end{eqnarray} \tag{ 25 }$

One immediate conclusion is that a pair of waves with wavevectors at a shallow enough angle to the k_x-axis tends to contribute energy to the mean flow, reinforcing it. The critical angle is given by $\tan ({\phi }_{\mathrm{crit}})=1/\sqrt{5}$, or ${\phi }_{\mathrm{crit}}\approx 24^\circ$.⁵ A pair of waves with an angle greater than ${\phi }_{\mathrm{crit}}$ draws energy from the mean flow, diminishing it.

We briefly comment on the use of periodic boundary conditions, infinite plane waves, and linear shear, which are mathematically convenient but could potentially raise some concern because of possible inconsistencies or physical subtleties. Within the literature, others have explored the use of more realistic profiles for the shear flow and perturbations, such as using wavepackets rather than infinite plane waves.

These more realistic profiles have not been found to fundamentally alter the direction of energetic transfer from those in simpler calculations. For instance, in a calculation involving perturbation growth in a finite-amplitude linear shear flow, Farrell (1987) used localized perturbation wavepackets and showed that similar conclusions about energetic changes are obtained as when infinite plane waves are used. Another calculation, more relevant to the present study as it is concerned with the growth of the mean flow, uses localized wavepackets and periodic, rather than linear, shear flow (Parker & Krommes 2018). That calculation found energy transfer to the shear flow, just as is found here.

We extend now the analysis of the Kelvin–Orr shearing wave to include magnetic fields. Again, we consider the energetics of the zonally averaged flow. We neglect $\overline{A}$, which is justified by the later numerical findings in Section 4.

Returning to Equation 10(a), retaining the magnetic fluctuations, and performing similar steps as in the previous subsection, we find the energetics of the mean zonal flow are now given by

$\begin{eqnarray}&&\displaystyle \frac{{{dE}}_{\mathrm{ZF}}}{{dt}}=\displaystyle \frac{1}{{L}_{y}}{\int }_{0}^{{L}_{y}}{dy}(\overline{u^{\prime} v^{\prime} }{\partial }_{y}\bar{u}-\overline{{B}_{x}^{{\prime} }{B}_{y}^{{\prime} }}{\partial }_{y}\bar{u}).\end{eqnarray} \tag{ 26 }$

We also need the generalization of the shearing wave that includes magnetic fields. With $\overline{A}=0$ and $\bar{u}={Sy}$, the linearized, non-forced equations for the perturbations $\zeta ^{\prime}$ and A' are

$\begin{eqnarray}&&({\partial }_{t}+{Sy}{\partial }_{x})\zeta ^{\prime} +\beta {\partial }_{x}\psi ^{\prime} =-{B}_{0}{\partial }_{x}{{\rm{\nabla }}}^{2}A^{\prime} +\nu {{\rm{\nabla }}}^{2}\zeta ^{\prime} ,\end{eqnarray} \tag{ 27a }$

$\begin{eqnarray}&&({\partial }_{t}+{Sy}{\partial }_{x})A^{\prime} =-{B}_{0}{\partial }_{x}\psi ^{\prime} +\eta {{\rm{\nabla }}}^{2}A^{\prime} .\end{eqnarray} \tag{ 27b }$

Assuming $\zeta ^{\prime} ({\boldsymbol{x}},t)=Z(t){e}^{i{\boldsymbol{k}}(t)\cdot {\boldsymbol{x}}}$ and $A^{\prime} ({\boldsymbol{x}},t)=a(t){e}^{i{\boldsymbol{k}}(t)\cdot {\boldsymbol{x}}}$, we find the same shearing dependence for ${\boldsymbol{k}}(t)$ as before (see Equation (20)). Then, we have

$\begin{eqnarray}&&\displaystyle \frac{{dZ}}{{dt}}=\displaystyle \frac{{{ik}}_{x}\beta }{k{(t)}^{2}}Z+{{ik}}_{x}{B}_{0}k{(t)}^{2}a-\nu k{(t)}^{2}Z,\end{eqnarray} \tag{ 28a }$

$\begin{eqnarray}&&\displaystyle \frac{{da}}{{dt}}=\displaystyle \frac{{{ik}}_{x}{B}_{0}}{k{(t)}^{2}}Z-\eta k{(t)}^{2}a.\end{eqnarray} \tag{ 28b }$

If k² did not depend on time, then these equations would be exactly the linearized equations without mean flow and would give rise to the fast and slow waves with frequencies ${\omega }_{f}$, ${\omega }_{s}$. In that case, the solution for any initial condition could be decomposed into the fast and slow eigenmode components. In particular, if the linear combination Z(0) and a(0) start off exactly in the fast eigenmode, then the time dependence of Z(t) and a(t) is given by $\exp (-i{\omega }_{f}t)$, where the imaginary part of ${\omega }_{f}$ determines the damping rate.

The shear flow complicates matters because k² now changes with time. However, when the shear is small, such that ${k}_{x}{St}/{k}_{y0}\ll 1$, k² remains nearly constant up through the decay time of the wave. Hence, the constant-k² solution of the previous paragraph is mostly retained. We expand k² to leading order in S. If a wave starts as an eigenmode, it will stay in that eigenmode to lowest order; the solution for Z is then given by

$\begin{eqnarray}&&Z(t)={Z}_{0}{e}^{-i\theta (t)}\exp \left[{\int }_{0}^{t}d\tau \ \mathrm{Im}\ \omega (\tau )\right],\end{eqnarray} \tag{ 29 }$

where $\theta (t)\mathop{=}\limits^{\mathrm{def}}{\int }_{0}^{t}d\tau \ \mathrm{Re}\ \omega (\tau )$ is some phase. An expression similar to Equation (29) also holds for a(t).

We now restrict ourselves to the parameter regime where νk², $\eta {k}^{2}$, ${\omega }_{{\rm{A}}}\ll {\omega }_{{\rm{R}}}$. The fast and slow frequencies ${\omega }_{f}$ and ${\omega }_{s}$ simplify to the expressions in Equation (7). In this limit, $\mathrm{Im}({\omega }_{f})=-\nu {k}^{2}$ and $\mathrm{Im}({\omega }_{s})=-\eta {k}^{2}$, and the wave damping behaves purely diffusively. When starting in the fast eigenmode, the solution for small shear is

$\begin{eqnarray}&&{Z}_{f}(t)={Z}_{0}{e}^{-i{\theta }_{f}(t)}\exp \left[-\nu {\int }_{0}^{t}d\tau \,k{(\tau )}^{2}\right],\end{eqnarray} \tag{ 30a }$

$\begin{eqnarray}&&{a}_{f}(t)={A}_{0}{e}^{-i{\theta }_{f}(t)}\exp \left[-\nu {\int }_{0}^{t}d\tau \,k{(\tau )}^{2}\right].\end{eqnarray} \tag{ 30b }$

The initial amplitudes Z₀ and A₀ are related by the eigenmode relation Equation (8). A similar expression exists for the slow wave, with ν replaced by η.

We use the shearing wave solution in Equation (30) to compute the energetic changes of the mean flow. For a single wave, the Reynolds stress $\overline{u^{\prime} v^{\prime} }$ is given by Equation (22); similarly, the Maxwell stress is given by

$\begin{eqnarray}&&\overline{{B}_{x}^{{\prime} }{B}_{y}^{\prime} }\to \displaystyle \frac{1}{2}\mathrm{Re}({B}_{x}^{{\prime} }{B}_{y}^{{\prime} * })=-\displaystyle \frac{1}{2}{k}_{x}{k}_{y}(t)| a(t){| }^{2}.\end{eqnarray} \tag{ 31 }$

Integrating over the lifetime of the wave, the net energy change in the mean flow due to a single wave shearing over is then

$\begin{eqnarray}&&{\rm{\Delta }}{E}_{\mathrm{ZF}}=-\displaystyle \frac{1}{2}{{Sk}}_{x}{\int }_{0}^{\infty }{dt}\left[\displaystyle \frac{{k}_{y}(t)}{k{(t)}^{4}}| Z(t){| }^{2}-{k}_{y}(t)| a(t){| }^{2}\right].\end{eqnarray} \tag{ 32 }$

We consider the effect of two (noninteracting) waves, with initial wavevectors $({k}_{x},{k}_{y0})$ and $(-{k}_{x},{k}_{y0})$. The procedure is much the same as in Section 3.1. Expanding to leading order in S, we obtain

$\begin{eqnarray}&&{({\rm{\Delta }}{E}_{\mathrm{ZF}})}_{f}=\displaystyle \frac{{S}^{2}{k}_{x}^{2}}{4{\nu }^{2}{k}^{4}}\left(\mathop{\underbrace{\displaystyle \frac{{k}_{x}^{2}-5{k}_{y}^{2}}{{k}^{6}}| {Z}_{0}{| }^{2}-\displaystyle \frac{{k}_{x}^{2}-{k}_{y}^{2}}{{k}^{2}}| {A}_{0}{| }^{2}}}\limits_{\mathop{=}\limits^{\mathrm{def}}J}\right).\end{eqnarray} \tag{ 33 }$

The corresponding expression for ${({\rm{\Delta }}{E}_{\mathrm{ZF}})}_{s}$ is identical with ν replaced by η.

Expression (33) generalizes the energy transfer to a weak mean flow due to Kelvin–Orr shearing wave dynamics to include magnetic fields. It is a major result of this paper. The term proportional to $| {Z}_{0}{| }^{2}$ stems from the Reynolds stress while the term proportional to $| {A}_{0}{| }^{2}$ comes from the Maxwell stress. We note that no explicit dependence on β or B₀ has yet appeared in ${\rm{\Delta }}{E}_{\mathrm{ZF}}$. Both β and B₀ have only an indirect effect on the size of the perturbations Z₀ and A₀.

Focusing on the wavevector dependence, we examine the quantity J inside the parentheses in Equation (33), which ${({\rm{\Delta }}{E}_{\mathrm{ZF}})}_{f,s}$ is proportional to. Substituting the energy-normalized eigenfunctions from Equation (8) and letting ${k}_{x}=k\cos \phi$ and ${k}_{y}=k\sin \phi$, we obtain

$\begin{eqnarray}&&{J}_{f}=({\cos }^{2}\phi -5{\sin }^{2}\phi )-\displaystyle \frac{{\omega }_{{\rm{A}}}^{2}}{{\omega }_{{\rm{R}}}^{2}}({\cos }^{2}\phi -{\sin }^{2}\phi ),\end{eqnarray} \tag{ 34a }$

$\begin{eqnarray}&&{J}_{s}=\displaystyle \frac{{\omega }_{{\rm{A}}}^{2}}{{\omega }_{{\rm{R}}}^{2}}({\cos }^{2}\phi -5{\sin }^{2}\phi )-({\cos }^{2}\phi -{\sin }^{2}\phi ),\end{eqnarray} \tag{ 34b }$

for the fast and slow wave, respectively.

We make several observations. First, for the fast wave, J_f = 0 determines the critical angle ϕ_crit that separates the waves that drive the mean flow from those that suppress it. As mentioned before, without magnetic field, ${\phi }_{\mathrm{crit}}\approx 24^\circ$. However, from Equation 34(a), we see that turning on the magnetic field causes the second term to become nonzero. Increasing the magnetic field reduces the critical angle, implying that now a smaller subset of fast-wave perturbations can contribute positively toward the growth of the shear flow. Figure 1(a) shows how the critical angle varies with the background magnetic field B₀. Waves with $\phi \lt {\phi }_{\mathrm{crit}}$ are of primary interest because these fast waves contribute positively to ${\rm{\Delta }}{E}_{\mathrm{ZF}}$, potentially driving strong growth of the mean flow. Second, J_s, like J_f, can be of either sign. However, for those waves with $\phi \lt {\phi }_{\mathrm{crit}}$, the slow wave opposes the mean flow, i.e., ${J}_{s}\lt 0$. Third, for $\phi \lt {\phi }_{\mathrm{crit}}$, the magnitude of J_f decreases as the magnetic field increases. The magnitude of J_s also somewhat decreases (see Figure 1(b)).

Zoom In Zoom Out Reset image size

The Kelvin–Orr shearing wave calculation does not capture the relative fraction of energy that resides in magnetic fluctuations compared with hydrodynamic fluctuations. Rather, the strength of these fluctuations, Z₀ and A₀, are taken here as given. In the parameter regime we have examined, magnetic fluctuations reside primarily in the slow wave and hydrodynamic fluctuations reside in the fast wave. Because Z₀ and A₀ are exogenous to this calculation, the use of energy-normalized eigenfunctions eases the interpretation of the physics by separating the effect of the wave from the amount of energy contained in each wave. Intuitively, and as we shall see later, as B₀ increases, more energy resides in the magnetic fluctuations and the slow wave more strongly suppresses the growth of zonal flow.

We have assumed an initial condition that starts off as a pure fast or slow wave, and calculated the effects of the two waves separately. Mathematically, this is equivalent to neglecting cross terms in the Reynolds and Maxwell stresses, which are quadratically nonlinear. From a physical point of view, this amounts to an assumption that the interaction between waves is negligible.

To summarize this section, we generalized the Kelvin–Orr shearing wave for a weak shear flow to include magnetic fields. We obtained Equation (33), one of the major results of this article, which describes a mean shear flow's energetic change due to a pair of shearing waves. Our calculation shows that magnetic fluctuations, through the slow magneto-Rossby wave, will oppose the growth of a mean shear flow. An additional effect is that a stronger B₀ also reduces the fast wave's contribution to driving a mean flow. We shall see later that the former is the dominant effect (see Figure 5(b) and surrounding discussion).

The Kelvin–Orr calculation is not a complete description because it does not close the loop and say how the zonal flow dynamically evolves. Furthermore, the computation is limited to long-wavelength shear flows. It also does not provide a growth rate. But it does give a clear physical picture of the effect of a weak shear flow on fluctuations, and shows, unambiguously, that a magnetic field opposes the growth of zonal flows. This simple calculation also quantitatively predicts which wavevectors contribute to driving or suppressing zonation.

The next section includes a more detailed and elaborate computation that is both dynamically consistent and also is not limited to long-wavelength mean flows. We shall see that the key conclusions of the wavenumber dependence of the Reynolds and Maxwell stress found in this simple Kelvin–Orr calculation (Equation (33)) are recovered from the more consistent calculation of the next section, in the appropriate asymptotic limit.

We present results from the ZI analysis. We consider a domain of size $2\pi \times 2\pi$, use parameter values $\beta =2$, $\nu \,=\eta ={10}^{-4}$, and take isotropic forcing centered about a total wavenumber k_f. That is:

$\begin{eqnarray}&&{\hat{Q}}_{{\boldsymbol{k}}}={Q}_{0}\,{e}^{-{(k-{k}_{f})}^{2}/(2\delta {k}_{f}^{2})},\end{eqnarray} \tag{ 35 }$

where

$\begin{eqnarray}&&{Q}_{0}=5\times {10}^{-5},{k}_{f}=12,\mathrm{and}\ \delta {k}_{f}=1.5.\end{eqnarray} \tag{ 36 }$

This forcing injects energy into hydrodynamic fluctuations at a rate $\epsilon ={\sum }_{{\boldsymbol{k}}}{\hat{Q}}_{{\boldsymbol{k}}}/(2{k}^{2})=4.81\,\times \,{10}^{-5}$. The forcing introduces a length scale ${k}_{f}^{-1}$ and a timescale ${(\epsilon {k}_{f}^{2})}^{-1/3}$.

For each q, there are multiple eigenmodes, each with its own ZI eigenvalue λ. Figure 2 shows the eigenvalue with maximum growth rate as a function of the mean flow wavenumber q for various values of the strength of the background magnetic field B₀ (normalized as $| {\omega }_{{\rm{A}}}/{\omega }_{{\rm{R}}}|$). As B₀ increases, the ZI is inhibited. This inhibition is also seen in Figure 3 in which the eigenvalue λ is shown as a function of the magnetic field strength for fixed mean-field wavenumber q.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When there is instability, the mean-flow components of the eigenfunction consists primarily of mean zonal jet $\delta \bar{u}$ rather than mean magnetic field $\delta \bar{A};$ see Figure 2(c). That the mean flow eigenfunction is dominated by $\delta \bar{u}$ is a general characteristic of the ZI of Equation (1), at least in all parameter ranges we have explored. The smallness of the mean magnetic component compared to the mean flow justifies our choice in the Kelvin–Orr calculation (Section 3.2) to use only a mean shear flow and to neglect a mean sheared magnetic field.

When the ZI is robustly strong—typically at low values of the magnetic field, $| {\omega }_{{\rm{A}}}/{\omega }_{{\rm{R}}}| \lesssim 0.25$—the eigenvalue is typically real. As the magnetic field becomes stronger, not only does the growth rate drop considerably, but also the eigenvalue becomes complex; this is seen in both Figures 2 and 3. While our ZI calculation is only linear and does not predict the final nonlinearly saturated state, the physics of a stationary (real eigenvalue) and translating (complex eigenvalue) mode can be quite different, and it is useful to distinguish between these cases. For instance, it is possible that the growing mode with real eigenvalue saturates into stationary zonal flows, while the mode with complex eigenvalue does not.

We can gain insight into how the ZI is inhibited by examining the Reynolds and Maxwell stresses for the eigenmodes.⁶ Recalling the zonally averaged momentum Equation 10(a), the Reynolds and Maxwell stresses are the fluctuation-driven terms that can drive or oppose the growth of the mean flow.

The perturbation equation for the mean flow eigenmode is described by

$\begin{eqnarray}&&(\lambda +\nu {q}^{2})\delta \bar{u}-\delta \overline{v^{\prime} \zeta ^{\prime} }+\delta \overline{({\partial }_{x}A^{\prime} ){{\rm{\nabla }}}^{2}A^{\prime} }=0,\end{eqnarray} \tag{ 37 }$

which comes directly from Equation 10(a). To a good approximation the above simplifies to

$\begin{eqnarray}&&\lambda +\nu {q}^{2}-\delta {\overline{v^{\prime} \zeta ^{\prime} }}^{u}+\delta {\overline{({\partial }_{x}A^{\prime} ){{\rm{\nabla }}}^{2}A^{\prime} }}^{u}=0.\end{eqnarray} \tag{ 38 }$

Here, the u superscript refers to the parts of the stresses associated with the perturbation mean flow $\delta \bar{u}$, neglecting the contribution associated with the perturbation mean magnetic field $\delta \bar{A}$. This decomposition of the stresses into components associated with $\delta \bar{u}$ and $\delta \bar{A}$ emerges from the instability calculation detailed in the Appendix. Because the mean magnetic component of the eigenfunction is small, $\delta \overline{v^{\prime} \zeta ^{\prime} }\,\approx {e}^{{iqy}}\,\delta {\overline{v^{\prime} \zeta ^{\prime} }}^{u}$ and $\delta \overline{({\partial }_{x}A^{\prime} ){{\rm{\nabla }}}^{2}A^{\prime} }\approx {e}^{{iqy}}\,\delta {\overline{({\partial }_{x}A^{\prime} ){{\rm{\nabla }}}^{2}A^{\prime} }}^{u}$.

Figure 4 shows the fluctuation stresses $\delta {\overline{v^{\prime} \zeta ^{\prime} }}^{u}$ and $\delta {\overline{({\partial }_{x}A^{\prime} ){{\rm{\nabla }}}^{2}A^{\prime} }}^{u}$. Panels (a)–(c) show the Reynolds stress and Maxwell stress for a marginally stable eigenmode $\lambda =0$ at three different values of the magnetic field. In this figure, positive values of the Reynolds stress reinforce the zonal flow, while positive values of the Maxwell stress oppose the zonal flow. At zero magnetic field (panel (a)), the Maxwell stress is zero and the Reynolds stress drives growth. At moderate magnetic field (panel (b)), the Maxwell stress is nonzero and opposes the zonal flow, but it does not have a significant effect because it is still considerably less than the Reynolds stress. At a large magnetic field (panel (c)), the Maxwell stress has grown such that it is almost exactly equal to the Reynolds stress. The Maxwell stress completely counteracts the driving effect of the Reynolds stress.

Zoom In Zoom Out Reset image size

It is also possible to take a closer look and examine the spectral decomposition of the Reynolds and Maxwell stresses. Considering marginally stable modes has been a useful way in earlier studies of ZI in unmagnetized fluids to understand which of the spectral components of the forcing contribute to ZI (Bakas & Ioannou 2013a; Bakas et al. 2015). Using analytic formulas derived in the course of the ZI calculation, we can extract the contribution of individual Fourier modes to the stresses. The procedure to obtain these analytic formulas is described in the Appendix, but the formulas themselves are not written explicitly because they are extremely complicated.

We can thus determine which fluctuation wavevectors tend to contribute positively or negatively toward the Reynolds and Maxwell stress. Figures 4(d)–(i) depict the spectral decomposition of marginally stable eigenmodes on a $(q,\phi )$ polar grid. For example, for the case with ${B}_{0}=0$, panel (d) implies that when a mean-flow perturbation $\delta \bar{u}$ with wavenumber $q/{k}_{f}=0.4$ is introduced in the flow, the forcing components ${\boldsymbol{k}}=(\pm {k}_{f},0)$ will induce Reynolds stresses with $\delta {\overline{v^{\prime} \zeta ^{\prime} }}^{u}\approx 0.08{(\epsilon {k}_{f}^{2})}^{1/3}\gt 0$ that tend to reinforce $\delta \bar{u}$, leading to instability.

We can see that for small values of the background magnetic field, the contribution of each component of the forcing to the Reynolds stresses remains mostly unchanged. In other words, panel (e) is mostly the same as panel (d). On the other hand, panel (h) shows the spectral decomposition of the Maxwell stress at a moderate magnetic field. For high values of the magnetic field, comparison of panels (f) and (i) shows that the cancellation between Reynolds stresses and Maxwell stresses occurs at each individual wavevector.

At this point, we can make close connection with the Kelvin–Orr shearing wave calculation presented in Section 3.2. The analytic formulas used for the spectral decompositions of the Reynolds and Maxwell stresses in Figure 4 can be asymptotically expanded in a limit relevant to the Kelvin–Orr shearing wave. The limit consistent with the Kelvin–Orr calculation is to take $\lambda \to 0$, small B₀, and small q. In this limit, the leading order terms for the Reynolds and Maxwell stresses are

$\begin{eqnarray}&&\delta {\overline{v^{\prime} \zeta ^{\prime} }}^{u}=\displaystyle \frac{2{k}_{x}^{2}{q}^{2}({k}_{x}^{2}-5{k}_{y}^{2})}{\nu {k}^{8}}{\hat{W}}_{{\boldsymbol{k}}}^{H},\end{eqnarray} \tag{ 39a }$

$\begin{eqnarray}&&\delta {\overline{({\partial }_{x}A^{\prime} ){{\rm{\nabla }}}^{2}A^{\prime} }}^{u}=\displaystyle \frac{2{k}_{x}^{2}{q}^{2}({k}_{x}^{2}-{k}_{y}^{2})}{\eta {k}^{4}}{\hat{G}}_{{\boldsymbol{k}}}^{H}.\end{eqnarray} \tag{ 39b }$

The parameter scalings for the Reynolds and Maxwell stresses are remarkably similar to those in Equation (33). (Recall that the first term on the right-hand side of Equation (33) arises from the Reynolds stress, and the second term from the Maxwell stress.) In particular, the wavevector dependence that determines positive versus negative contribution, $({k}_{x}^{2}-5{k}_{y}^{2})$ for the Reynolds stress and $({k}_{x}^{2}-{k}_{y}^{2})$ for the Maxwell stress, is exactly the same in the asymptotic limit of the ZI calculation and in the Kelvin–Orr calculation. The computer algebra system Mathematica was used both to derive the expressions for the stresses and to take the asymptotic limit.

Figure 4 shows, roughly, how small q must be (i.e., how long wavelength the zonal flow must be) for this asymptotic limit to be accurate. For example, we see that for $q/{k}_{f}\lesssim 0.2$, the Reynolds stresses are positive only for $\phi \lt 24^\circ ;$ similarly, the Maxwell stresses are positive for $\phi \lt 45^\circ$. For $q/{k}_{f}\gtrsim 0.2$, the constant-angle boundary (dash–dotted line) between positive and negative stresses is no longer accurate.

Figure 5 shows the balance between the Reynolds and Maxwell stresses as the resistivity η changes. As η changes from large to small, the Maxwell stress grows larger (Figure 5(a)). In Figure 5(b), we see that the Maxwell stress grows at the same rate as the overall level of magnetic fluctuations, as measured by the magnetic energy stored in the covariance G^H of the CE2 homogeneous equilibrium. At large η, for which magnetic fluctuations are suppressed, strong ZI occurs and the growth rate is about 0.4( k_f²)^(1/3). As η decreases and the level of magnetic fluctuations grows, eventually the Maxwell stress becomes comparable to the Reynolds stress, and the ZI is suppressed, with the growth rate weakening considerably. The eigenvalue λ and the stresses even become complex at $\eta ={10}^{-8}$, whereas these quantities are real for larger η.

Zoom In Zoom Out Reset image size

Figures 6(a)–(c) show the behavior of ZI on an $(\eta ,{B}_{0})$ grid. For each parameter value, a marker depicts whether the homogeneous equilibrium leads to growing, stationary ZF (ZI eigenvalue λ is real and positive, plus + signs), no growing ZF (λ is real and negative, circles ◦), or something indeterminate (λ is complex, often with positive real part, asterisks ∗). For these plots, only η and B₀ change while all other parameters are kept the same. Only a single ZF wavenumber q = 6 is used, which is typically close to the most unstable wavenumber. Figures 6(a) and (b) use the same parameters except the amplitude of the input forcing Q₀ is varied. Figure 6(c) uses a different value of ν.

Zoom In Zoom Out Reset image size

Up to some maximum B₀, the boundary in $(\eta ,{B}_{0})$ space between the growing, stationary zonal flow and the other behaviors is fitted well by a line $\eta /{B}_{0}^{2}=\mathrm{constant}$, which was also found by Tobias et al. (2007). The parameters of Figures 6(a) and (b) are chosen to match those of the simulations performed by Tobias et al. (2007), the results of which are summarized in Figure 7 (figure reproduced from paper by Tobias et al. 2007). However, we could not match the amplitude and spectral distribution of the input forcing exactly, as these values were not reported in detail. Despite an imperfect matching of forcings, there is nevertheless remarkable agreement between our findings, which result from examining only the ZI within a quasilinear theory, and the results from the fully nonlinear direct numerical simulations by Tobias et al. (2007). Part of the reason for this success is that within the ZI calculation, the details of the forcing turn out not that important. As we have argued in Sections 3 and 4, zonal jet appearance is controlled by the competition between the drive (Reynolds stresses) and suppressor (Maxwell stresses). The amplitude of the forcing, though, does not control this difference since both Reynolds and Maxwell stresses are proportional to the total energy input rate by the forcing. For example, compare Figures 6(a) and (b), which use the same input parameters except for a forcing strength that differs by two orders of magnitude. Qualitatively and quantitatively, the zonation boundary separating robust zonal flow growth (plus signs) from other behavior (circles and asterisks) changes little.

Zoom In Zoom Out Reset image size

Also shown in each plot of Figure 6 is a black contour, which depicts the curve ${({\omega }_{{\rm{A}}}^{2}/{\omega }_{{\rm{R}}}^{2})(1+{\Pr }_{{\rm{m}}})}^{2}/{\Pr }_{{\rm{m}}}=1$. To compute a single number for ${\omega }_{{\rm{A}}}^{2}/{\omega }_{{\rm{R}}}^{2}$, we use a characteristic wavenumber, which we take to be the forcing wavenumber k_f. In the regime ${\Pr }_{{\rm{m}}}\gg 1$, or $\nu \gg \eta$ (the bottom half of the curve), this curve reduces to $({\omega }_{{\rm{A}}}^{2}/{\omega }_{{\rm{R}}}^{2}){\Pr }_{{\rm{m}}}=1$. This equation recovers the observed scaling ${B}_{0}^{2}/\eta =\mathrm{constant}$, but also provides a value for the constant. As seen in Figure 6, this constant works remarkably well at disparate values of ν (separated by four orders of magnitude) at determining the $\eta /{B}_{0}^{2}$ boundary.

The parameter

$\begin{eqnarray}&&{\rm{\Upsilon }}\mathop{=}\limits^{\mathrm{def}}({\omega }_{{\rm{A}}}^{2}/{\omega }_{{\rm{R}}}^{2}){(1+{\Pr }_{{\rm{m}}})}^{2}/{\Pr }_{{\rm{m}}}\end{eqnarray} \tag{ 40 }$

is derived from the level of magnetic fluctuations in the homogeneous equilibrium G^H. The expression is given in Equation (44). A key parameter determining the homogeneous equilibrium is

$\begin{eqnarray}&&z\mathop{=}\limits^{\mathrm{def}}{\omega }_{{\rm{R}}}^{2}+{(\nu +\eta )}^{2}{k}^{4}+\displaystyle \frac{{(\nu +\eta )}^{2}}{\nu \eta }{\omega }_{{\rm{A}}}^{2}.\end{eqnarray} \tag{ 41 }$

In the regime of $\nu {k}^{2},\,\eta {k}^{2},\,{\omega }_{{\rm{A}}}\ll {\omega }_{{\rm{R}}}$, the middle term of z is negligible. The third term can be large or small compared to ${\omega }_{{\rm{R}}}^{2}$ because ${(\nu +\eta )}^{2}/\nu \eta ={(1+{\Pr }_{{\rm{m}}})}^{2}/{\Pr }_{{\rm{m}}}$ can be big if either of ν or η is much larger than the other. The critical parameter ${\rm{\Upsilon }}$ is the ratio of the third term to the first term. If ${\Pr }_{{\rm{m}}}$ is not too large or too small such that ${\rm{\Upsilon }}\ll 1$, then $z\approx {\omega }_{{\rm{R}}}^{2}$. Furthermore, if the additional assumption is made that ${\Pr }_{{\rm{m}}}\gg 1$ but still ${\rm{\Upsilon }}\ll 1$, the covariance of magnetic fluctuations becomes, from Equation 44(d),

$\begin{eqnarray}&&{\hat{G}}_{{\boldsymbol{k}}}^{H}=\displaystyle \frac{{\omega }_{{\rm{A}}}^{2}}{{\omega }_{{\rm{R}}}^{2}}\displaystyle \frac{{\hat{Q}}_{{\boldsymbol{k}}}}{2\eta {k}^{6}}.\end{eqnarray} \tag{ 42 }$

Hence, the covariance of magnetic fluctuations scales as B₀²/η, while in the same regime, the covariance of hydrodynamic fluctuations ${\hat{W}}^{H}$ is independent of both η and B₀. Thus, we have related parameters that determine the magnetic fluctuation level to the boundary of zonostrophic instability and found good agreement. The precise physics determining ${\rm{\Upsilon }}=1$ as a critical value (when ${\Pr }_{{\rm{m}}}\gg 1$) are not fully understood. However, the agreement between ${\rm{\Upsilon }}=1$ and the zonostrophic instability boundary is broadly consistent with the idea that magnetic fluctuations oppose zonostrophic instability, and hence suppress zonal flow.