2.1. Governing equations

We start with the two-dimensional vorticity equation on the $\beta$ plane:

(2.1)\begin{equation} \zeta_t-\psi_y\zeta_x+\psi_x\zeta_y+\beta\psi_x=F-\mu\zeta+\nu \nabla^2 \zeta,\quad \zeta=\nabla^2\psi. \end{equation}

Here the $x$- and $y$-directions are the longitudinal and latitudinal directions, respectively, $\psi$ is the stream function, $\zeta$ is the vorticity, $F$ is an external forcing which we will specify later, $\mu$ is the bottom drag and $\nu$ is the viscosity. The variables have been non-dimensionalised by characteristic length and time scales.

Despite the simplicity of (2.1), it is too difficult for an analysis of zonostrophic instability as it is nonlinear. We make the ‘quasilinear approximation’ to proceed analytically. We decompose the flow as

(2.2a,b)\begin{equation} \psi=\bar{\psi}(y,t)+\psi'(x,y,t),\quad \zeta=\bar{\zeta}(y,t)+\zeta'(x,y,t), \end{equation}

where the overline represents the zonal average,

(2.3a,b)\begin{equation} \bar{\psi}=\frac{k}{2{\rm \pi}}\int_0^{{2{\rm \pi}}/{k}}\psi\, \mathrm{d}x,\quad \bar{\zeta}=\frac{k}{2{\rm \pi}}\int_0^{{2{\rm \pi}}/{k}}\zeta\, \mathrm{d}x, \end{equation}

and $\psi '$ and $\zeta '$ are the fluctuating fields. Substituting (2.2a,b) into (2.1), we neglect all the nonlinear terms except those with mean components. Using $U=-\bar {\psi }_y$ to denote the mean zonal velocity, we derive the equation for the fluctuating vorticity,

(2.4)\begin{equation} \zeta'_t+U\zeta'_x+(\beta-U_{yy})\psi_x'=F-\mu\zeta'+\nu\nabla^2\zeta',\quad \zeta'=\nabla^2\psi'. \end{equation}

To find the evolution equation for $U$, we take the zonal average of (2.1) and find

(2.5)\begin{equation} U_t-(\overline{\psi'_x\psi'_y})_y={-}\mu U+\nu U_{yy}, \end{equation}

taking the forcing $F$ to have zero zonal average. Equations (2.4) and (2.5) constitute the approximate quasilinear system that we will use to study the zonostrophic instability analytically. We will also undertake numerical simulations to test the theory and the validity of the quasilinear approximation using the full nonlinear equation (2.1).

We consider the forcing

(2.6)\begin{equation} F=\sigma \hat{\xi}(t)\,\exp(\mathrm{i} kx)+\mathrm{c.c.}, \end{equation}

where $\sigma$ is the strength of the forcing, $k$ is the wavenumber in the $x$-direction, $\hat {\xi }$ is a complex Gaussian white noise and ‘c.c.’ represents the complex conjugate of the previous term (or terms). The expectation of the white noise has the properties

(2.7a–c)\begin{equation} \mathbb{E}[\hat{\xi}(t_1)\hat{\xi}^*(t_2)]=\delta(t_1-t_2),\quad \mathbb{E}[\hat{\xi}(t_1)\hat{\xi}(t_2)]=0,\quad \mathbb{E}[\hat{\xi}(t)]=0, \end{equation}

which indicate that the values of the white noise at two different times are independent (2.7a), it has zero expectation (2.7c) and its statistical properties are independent of time. Because of these properties, white noise forcing is a standard idealisation for stochastic differential equations. In our problem, we use it to drive waves with random amplitudes, as a very idealised model of turbulence. More details of the white noise and our numerical implementation are given in Appendix A.

The forcing (2.6) is the same as that used by Farrell & Ioannou (Reference Farrell and Ioannou2003,  Reference Farrell and Ioannou2007). It also has similarities with the approach of Srinivasan & Young (Reference Srinivasan and Young2012), who use a ‘ring forcing’, namely the stochastic driving of a ring of wavenumbers of given radius in Fourier space (for details, see Appendix B). The ring forcing is essentially isotropic, so the $\beta$-effect is necessary for the zonostrophic instability. This property of isotropy may be more relevant to modelling the formation of zonal jets on planets. Our forcing (2.6), on the other hand, has a single wavenumber in the $x$-direction whilst remaining uniform in the $y$-direction, and is therefore a ‘point forcing’ which is anisotropic. Although not as realistic as the isotropic ring forcing, it can still reveal key properties of zonostrophic instability. For MHD instabilities, because of the complexity of the analysis, we simplify to a point forcing as a first step, and leave any elaboration to an isotropic ring forcing for future research.

2.2. The ‘basic state’

We now consider the stability problem governed by (2.4) and (2.5). Because of the white-noise forcing, $\psi$ is stochastic, which complicates typical stability analyses. The usual method to remove stochasticity is to use a spatiotemporal correlation function which will evolve deterministically (Farrell & Ioannou Reference Farrell and Ioannou2003,  Reference Farrell and Ioannou2007; Srinivasan & Young Reference Srinivasan and Young2012). We take a different approach: we directly solve for $\psi$ in terms of the white noise and so retain its randomness. When we proceed to the equation for the mean flow, where quadratic terms of the fundamental waves will appear, we compute the expectation, taking advantage of the properties of the white noise given in (2.7a–c).

The forcing $F$ in (2.6) is independent of $y$ and generates waves. We consider a state in which there is zero zonal mean flow $U$ as the ‘basic state’, upon which zonostrophic instability develops. We thus take a solution of the form

(2.8)\begin{equation} \psi'=\psi_1=\skew{3}{\hat}{\psi}_1(t)\exp (\textrm{i} kx)+\mathrm{c.c.},\quad U=0. \end{equation}

Substitution into (2.4) yields

(2.9)\begin{equation} \frac{\mathrm{d}\skew{3}{\hat}{\psi}_1}{\mathrm{d}t}+ \left(-\mathrm{i}\frac{\beta}{k}+\mu+\nu k^2\right)\skew{3}{\hat}{\psi}_1={-}\frac{\sigma }{k^2}\hat{\xi}. \end{equation}

The solution for $\hat {\psi }_1$ is then

(2.10)\begin{equation} \skew{3}{\hat}{\psi}_1={-}\frac{\sigma}{k^2}\exp(\lambda_1 t)\int^t\hat{\xi}(\tau) \exp({-\lambda_1\tau})\,\mathrm{d}\tau, \end{equation}

where

(2.11)\begin{equation} \lambda_1={-}\mu-\nu k^2+\mathrm{i}\frac{\beta}{k} \end{equation}

is the eigenvalue of the homogeneous system (2.9). We have omitted the lower integral limit to relax the requirement on the initial condition, which has negligible effect on the long-time behaviour of $\hat {\psi }_1$ given the damping $\mathrm {Re}(\lambda _1)<0$.

The fluctuating flow $\hat {\psi }_1$ has the form of a damped Rossby wave, with its amplitude driven by the white noise. It is unsteady and stochastic, but statistically, $\hat {\psi }_1$ has zero expectation and its probability density will settle down to a steady distribution as $t\rightarrow \infty$, as can be checked by solving the corresponding Fokker–Planck equation.

2.3. Instability analysis

We then study the stability of the basic state in the statistical sense, adding now small perturbations $\psi _2$ and $U_2$:

(2.12a,b)\begin{equation} \psi'=\psi_1+\psi_2 + \cdots,\quad U=U_2 + \cdots . \end{equation}

Assuming that $\psi _2$ and $U_2$ are small, we substitute (2.12a,b) into (2.4) and (2.5), and then linearise terms involving $\psi _2$ and $U_2$, which gives

(2.13)$$\begin{gather} \zeta_{2,t}+\beta\psi_{2,x}+\mu\zeta_2-\nu\nabla^2\zeta_2={-}U_2 \zeta_{1,x}+U_{2,yy}\psi_{1,x}\quad \zeta_2=\nabla^2\psi_2, \end{gather}$$

(2.14)$$\begin{gather}U_{2,t}+\mu U_2-\nu U_{2,yy}=(\overline{\psi_{1,y}\psi_{2,x}+\psi_{2,y}\psi_{1,x}})_y. \end{gather}$$

We seek solutions in the form

(2.15a)\begin{align} U_2&=\hat{U}(t)\exp (\textrm{i} my) +\mathrm{c.c.}, \end{align}

(2.15b)\begin{align} \psi_2&=\skew{3}{\hat}{\psi}_{21}(t) \exp({\mathrm{i}kx+\mathrm{i}my})+ \skew{3}{\hat}{\psi}_{22}(t) \exp({-\mathrm{i}kx+\mathrm{i}my})+\mathrm{c.c.}, \end{align}

that is, the mean flow $U_2$ has a wavenumber $m$ in the transverse direction, and $\psi _2$ has a wavenumber combination of the basic wave and the mean flow. Substituting into (2.13) and (2.14), we obtain

(2.16)$$\begin{gather} \frac{\mathrm{d}\skew{3}{\hat}{\psi}_{21}}{\mathrm{d} t}-\lambda_2\skew{3}{\hat}{\psi}_{21}={-}\mathrm{i}k \varLambda \hat{U}\skew{3}{\hat}{\psi}_1, \end{gather}$$

(2.17)$$\begin{gather}\frac{\mathrm{d}\skew{3}{\hat}{\psi}_{22}}{\mathrm{d} t}-\lambda_2^*\skew{3}{\hat}{\psi}_{22}=\mathrm{i}k \varLambda \hat{U}\skew{3}{\hat}{\psi}_1^*, \end{gather}$$

(2.18)$$\begin{gather}\frac{\mathrm{d}\hat{U}}{\mathrm{d}t}+(\mu+m^2\nu)\hat{U}=\mathrm{i}m^2k(\skew{3}{\hat}{\psi}_{21} \skew{3}{\hat}{\psi}_1^*-\skew{3}{\hat}{\psi}_{22}\skew{3}{\hat}{\psi}_1), \end{gather}$$

where

(2.19a,b)\begin{equation} \lambda_2={-}\mu-\nu(k^2+m^2)+\frac{\mathrm{i}k\beta}{k^2+m^2},\quad \varLambda= \frac{k^2-m^2}{k^2+m^2}. \end{equation}

The solutions of (2.16) and (2.17) are

(2.20)$$\begin{gather} \skew{3}{\hat}{\psi}_{21}={-}\mathrm{i}k\varLambda \exp(\lambda_2 t) \int^t\hat{U}(\tau)\skew{3}{\hat}{\psi}_1(\tau) \exp({-\lambda_2 \tau})\,\mathrm{d}\tau, \end{gather}$$

(2.21)$$\begin{gather}\skew{3}{\hat}{\psi}_{22}=\mathrm{i}k\varLambda \exp({\lambda_2^* t}) \int^t\hat{U}(\tau)\skew{3}{\hat}{\psi}_1^*(\tau) \exp({-\lambda_2^* \tau})\,\mathrm{d}\tau. \end{gather}$$

We have omitted the lower integral limit again as we may ignore initial conditions in determining exponential growth rates for the zonostrophic instability.

We substitute (2.20) and (2.21) into (2.18), and then obtain an equation for $\hat {U}$:

(2.22) \begin{align} & \frac{\mathrm{d}\hat{U}}{\mathrm{d}t}+(\mu+\nu m^2)\hat{U} \nonumber\\ &\quad =m^2k^2\varLambda\int^t\hat{U}(\tau)\!\left\{\skew{3}{\hat}{\psi}_1(\tau) \skew{3}{\hat}{\psi}_1^*(t) \exp[{-\lambda_2(\tau-t)}]+\skew{3}{\hat}{\psi}_1^*(\tau) \skew{3}{\hat}{\psi}_1(t)\exp[{-\lambda_2^*(\tau-t)}]\right\}\mathrm{d}\tau. \end{align}

This equation describes the evolution of the mean flow driven by Rossby waves. The quantity $\hat {U}$ is again stochastic, but our interest is its expectation $\mathbb {E}(\hat {U})$. In particular, we look for the solution that grows exponentially, to indicate zonostrophic instability. Thus, we consider

(2.23)\begin{equation} \mathbb{E}[\hat{U}(t)]=\tilde{U}\exp (st), \end{equation}

and the main objective of the analysis is to solve for the growth rate $s$.

In order to proceed mathematically, we need to make the assumption that the mean flow $\hat {U}(\tau )$ and the quadratic term of the fundamental wave $\hat {\psi }_1^*(\tau )\hat {\psi }_1(t)$ are statistically uncorrelated, hence their expectations are separable,

(2.24)\begin{equation} \mathbb{E}[\hat{U}(\tau)\skew{3}{\hat}{\psi}_1(\tau)\skew{3}{\hat}{\psi}_1^*(t)]= \mathbb{E}[\hat{U}(\tau)] \mathbb{E}[\skew{3}{\hat}{\psi}_1(\tau)\skew{3}{\hat}{\psi}_1^*(t)], \end{equation}

thus the expectation of (2.22) becomes

(2.25) \begin{align} & \frac{\mathrm{d}\mathbb{E}(\hat{U})}{\mathrm{d}t}+(\mu+\nu m^2)\mathbb{E}(\hat{U}) \nonumber\\ &\quad =m^2k^2\varLambda\int^t\mathbb{E}[\hat{U}(\tau)]\left\{\mathbb{E}[\skew{3}{\hat}{\psi}_1(\tau) \skew{3}{\hat}{\psi}_1^*(t)] \exp[{-\lambda_2(\tau-t)}] \right.\nonumber\\ &\left.\qquad +\mathbb{E}[\skew{3}{\hat}{\psi}_1^*(\tau) \skew{3}{\hat}{\psi}_1(t)]\exp[{-\lambda_2^*(\tau-t)}]\right\}\,\mathrm{d}\tau. \end{align}

We do not have a proof for (2.24): it is an assumption that we will make to derive an analytical dispersion relation, which we will test through comparison with numerical simulations. But (2.24) is directly related to the assumption of zonal-mean ergodicity that has been widely used in mean-flow dynamics (e.g. Srinivasan & Young Reference Srinivasan and Young2012; Farrell & Ioannou Reference Farrell and Ioannou2003; Marston, Conover & Schneider Reference Marston, Conover and Schneider2008). This assumption states that the zonal-mean velocity of an individual realisation is equal to the ensemble average, or expectation

(2.26)\begin{equation} \hat{U}=\mathbb{E}[\hat{U}]; \end{equation}

in this case since $\hat {U}$ is no longer stochastic, (2.26) implies assumption (2.24). Indeed, we will show that our result of zonostrophic instability is consistent with the result of Srinivasan & Young (Reference Srinivasan and Young2012) based on the ergodic assumption. But unlike (2.26), our assumption (2.24) still retains stochasticity in the mean flow and is therefore a weaker assumption. We will refer to (2.26) as the ‘full ergodic assumption’ and (2.24) as the ‘partial ergodic assumption’ in what follows.

The expectation of the fundamental-wave term $\hat {\psi }_1^*(\tau )\hat {\psi }_1(t)$ in (2.25) may be computed explicitly: applying property (2.7a) to (2.10), we find that for $t>\tau$,

(2.27)\begin{align} \mathbb{E}[\skew{3}{\hat}{\psi}_1(\tau)\skew{3}{\hat}{\psi}_1^*(t)] &= \frac{\sigma^2}{k^4}\exp({\lambda_1^*t+\lambda_1\tau}) \int^{p=t}\int^{q=\tau}\mathbb{E}[\hat{\xi}^*(p)\hat{\xi}(q)]\exp({-\lambda_1^*p-\lambda_1q})\,\mathrm{d}q\,\mathrm{d}p \nonumber\\ &= \frac{\sigma^2}{k^4}\exp({\lambda_1^*t+\lambda_1\tau}) \int^{p=t}\int^{q=\tau}\delta(p-q)\exp({-\lambda_1^*p-\lambda_1q})\,\mathrm{d}q\,\mathrm{d}p \nonumber\\ &={-}\frac{\sigma^2}{k^4}\frac{1}{\lambda_1+\lambda_1^*} \exp[{\lambda_1^*(t-\tau)}]. \end{align}

Note that the lower integral limits have been neglected as their contribution is exponentially small for large $t$.

Substituting (2.23) and (2.27) into (2.25), we obtain the dispersion relation determining the growth rate $s$:

(2.28)\begin{equation} s+\mu+\nu m^2={-}\frac{m^2\sigma^2\varLambda}{k^2(\lambda_1+\lambda_1^*)} \left(\frac{1}{s-\lambda_1^*-\lambda_2}+\frac{1}{s-\lambda_1-\lambda_2^*}\right). \end{equation}

We rewrite the dispersion relation in terms of the original variables,

(2.29)\begin{align} s+\mu+\nu m^2=\frac{m^2\sigma^2(k^2-m^2)}{2 k^2(\mu +\nu k^2)} \frac{1}{(s+2\mu+2\nu k^2+\nu m^2)(k^2+m^2)+\mathrm{i}\beta m^2/k} +\mathrm{c.c.e.}s., \end{align}

where the expression ‘c.c.e.$s$.’ represents the previous terms with all quantities complex conjugated except $s$ (cf. (2.28)); this notation is helpful to give succinct equations in this paper, and we will use it repeatedly in the MHD case.

To reduce the number of independent quantities, we rescale according to

(2.30a–e)\begin{equation} s=s_\star \sigma^{{2}/{3}},\quad m=m_\star k, \quad \mu=\mu_\star \sigma^{{2}/{3}},\quad \nu=\nu_\star \sigma^{{2}/{3}}k^{{-}2},\quad \beta=\beta_\star k \sigma^{{2}/{3}}. \end{equation}

This corresponds to a non-dimensionalisation based on the forcing strength $\sigma$, the viscosity $\nu$, and the scale $k^{-1}$ of the forcing. The non-dimensional quantity obtained from the viscosity may be linked to a Grashof number given by $\mathrm {Gr} = \nu _\star ^{-2}$ and $\beta _\star$ is a non-dimensional measure of the strength of the $\beta$-effect on the same basis (see Childress, Kerswell & Gilbert Reference Childress, Kerswell and Gilbert2001; Durston & Gilbert Reference Durston and Gilbert2016). Under this rescaling the dispersion relation (2.29) becomes

(2.31)\begin{equation} s_\star{+}\mu_\star{+}\nu_\star m_\star^2=\frac{m_\star^2(1-m_\star^2)}{2(\mu_\star{+}\nu_\star)} \frac{1}{(s_\star{+}2\mu_\star{+}2\nu_\star{+}\nu_\star m_\star^2)(1+m_\star^2)+\mathrm{i} \beta_\star m_\star^2} +\mathrm{c.c.e.}s. \end{equation}

The rescaled parameters will be convenient for finding conditions for instability in the parameter space, and for deriving asymptotic expressions for the dispersion relation of MHD instabilities shown later on. However, when presenting general results, we will mainly use the original parameters which are more relevant to the physics.

2.5. Numerical simulation

To verify our theory for hydrodynamic zonostrophic instability, we undertake numerical simulations for the flow governed by (2.1), with further details in Appendix A. We use a pseudospectral method with a spatial discretisation of $256\times 256$ mesh points in the domain $[0, 2{\rm \pi} /k]\times [0, 2{\rm \pi} /m]$. For temporal discretisation, we use the Crank–Nicolson method, with the nonlinear terms advanced using Euler's method. Itô's interpretation is used for integration of the white noise. A very small time step of ${\rm \Delta} t=0.005$ is used for both the temporal evolution and discretisation for the white noise, to ensure a good approximation to the decorrelation property in (2.7a). We choose the parameters $\beta =5$, $\nu =10^{-4}$, $\mu =0$, $\sigma =0.05$, $k=16$ and $m=5$, corresponding to $\beta _\star =2.3$, $\nu _\star =0.19$, which have been used in figure 1. We take $\psi =10^{-4}\cos (my)$ at $t=0$ to render a very weak zonal flow $U=-\overline {\partial _y\psi }$ initially. Our main objective here is to verify the theory rather than aim for physical realism, hence we only consider one Fourier mode for the initial zonal flow to make comparison straightforward.

In figure 2, we show an example of the evolution of the vorticity field, alongside two snapshots of the zonal mean flow profiles. At earlier times $t\leq 80$, the vorticity field has the same pattern as the forcing $F$ (i.e. it is periodic in the $x$-direction and homogeneous in the $y$-direction). Then unstable zonal flows gradually grow stronger, causing bending in the $y$-direction visible at $t=100$. At $t=110$, we observe that nonlinear effects become significant, generating a localised zonal jet near $y=1$. The jet continues to grow stronger, and then saturates and undergoes nonlinear evolution. At $t=180$, another jet in the opposite direction is visible around $y=0.2$. Note that a simulation with a different white noise results in a different flow pattern, but with qualitatively the same features as seen in figure 2.

Figure 2. Numerical solution for the evolution of vorticity and zonal velocity without magnetic field: $\beta =5$; $k=16$; $\nu =10^{-4}$; $\mu =0$; $\sigma =0.05$; $m=5$; $\beta _\star =2.3$; $\nu _\star =0.19$; $m_\star =0.31$.

In figure 3, we plot the evolution of the zonal flow; figure 3(a) shows the density field of $U(y,t)$ in a Hovmöller diagram, giving the evolution of the mean flow with time on the horizontal axis (snapshots of $U(y,t)$ at $t=110$ and $t=180$ are shown in figure 2). The formation of the two jets is clearly seen. Note that unlike previous studies on zonal jets (e.g. Srinivasan & Young Reference Srinivasan and Young2012; Parker & Krommes Reference Parker and Krommes2013), we do not see jet-merging in our simulation. We think this may be due to the simple structure of the forcing that we use to drive the waves (cf. (2.6)).

Figure 3. Evolution of the zonal velocity $U$ for figure 2. Panel (a) is a Hovmöller diagram showing $U(y,t)$, and (b) shows the r.m.s. value of $U$. The solid line is the solution of the numerical simulation, and the red straight line is the exponential growth predicted by the dispersion relation (2.29), with growth rate $s=6.17\times 10^{-2}$.

In figure 3(b), we show the r.m.s. velocity with respect to $y$, i.e.

(2.33)\begin{equation} U_{rms}= \left({\frac{m}{2{\rm \pi}}\int_0^{{2{\rm \pi}}/{m}}U^2\,\mathrm{d}y}\right)^{1/2}. \end{equation}

We also plot the prediction of zonostrophic instability, showing good agreement between the theory and the simulation for $t\in [80,120]$ after some transient behaviour up to $t\simeq 80$. The good agreement justifies the theory, including the quasilinear approximation and the assumption (2.24) of neglecting wave–mean flow correlation, used in the analysis. To show the effect of different realisations of the forcing, in figure 4 we plot the evolution of $U_{rms}$ for five different examples of the driving white noise. The exponential growth of zonostrophic instability theory is plotted as a straight line. The $U_{rms}$ starts to grow at different times for different realisations, but the growth rates are all very close to the theoretical prediction.

Figure 4. Five realisations of $U_{rms}$ with different white-noise forcing; the parameters are the same as in figure 2. The red straight line represents the exponential growth predicted by the dispersion relation (2.29).

The simulation results allow us to make further comments on the full assumption of mean-flow ergodicity in (2.26), that is assuming the mean flow has no stochasticity. Our simulations show that this assumption is clearly not satisfied, since $U_{rms}$ is different for different realisations. But if our focus is on whether zonostrophic instability occurs or not, then this assumption seems to be reasonable, because all realisations have periods of growth at a similar rate as our prediction for the expectation, derived under the weaker partial assumption (2.24). The underlying reason for this property, i.e. that the mean flow exhibits similar growth rates in different realisations, however, remains elusive. Previous theoretical studies have indicated that full mean-flow ergodicity holds when the flow reaches a statistically steady state (Marston et al. Reference Marston, Conover and Schneider2008), when the drag $\mu$ is sufficiently strong (Bouchet et al. Reference Bouchet, Nardini and Tangarife2013), or when there is time scale or space scale separation between the mean flow and the waves (Bouchet et al. Reference Bouchet, Nardini and Tangarife2013; Durston & Gilbert Reference Durston and Gilbert2016). None of these conditions applies to our case: the flow is in a transient state during the instability, our simulations used zero drag $\mu =0$, and our mean flow and waves have similar length scales. More theoretical work is required to understand when full or partial mean-flow ergodicity holds.

3.1. Governing equations

Following the same methodology, we study the MHD zonostrophic instability – we include a constant magnetic field $B_0$ in the $x$-direction and study its effect. The original MHD equations for the two-dimensional flow are

(3.1)$$\begin{gather} \zeta_t-\psi_y\zeta_x+\psi_x\zeta_y+\beta\psi_x ={-}a_y j_x+a_x j_y+F-\mu\zeta+\nu \nabla^2 \zeta, \end{gather}$$

(3.2)$$\begin{gather}a_t-\psi_y a_x+\psi_x a_y =\eta \nabla^2 a, \end{gather}$$

(3.3a,b)$$\begin{gather}\zeta =\nabla^2\psi,\quad j =\nabla^2 a. \end{gather}$$

Here $a$ is the magnetic potential, i.e. the magnetic field is ${\boldsymbol {B}}=(-a_y, a_x, 0)$, and $j$ is the current density. We again apply the quasilinear approximation for this system. For the flow, we apply the same decomposition as (2.2a,b). For the magnetic field, we decompose the wave and mean by

(3.4a,b)\begin{equation} a=\bar{a}(y,t)+a'(x,y,t) + \cdots ,\quad j=\bar{j}(y,t)+j'(x,y,t) + \cdots , \end{equation}

and for the mean field we set

(3.5)\begin{equation} -\bar{a}_y=B_0+B(y,t), \end{equation}

where $B_0$ is the externally added constant mean field in the $x$-direction, while $B$ is a small variation of magnetic field averaged in the $x$-direction and caused by the flow, with $|B|\ll |B_0|$. For the analysis of zonostrophic instability, we again need the quasilinear approximation; we substitute (3.4a,b) and (3.5) into (3.1) and (3.2), and then only keep the nonlinear terms involving the flow velocity $U$ and field $B$. The governing equations for the waves are

(3.6)$$\begin{gather} \zeta'_t+U\zeta'_x+(\beta-U_{yy})\psi'_x-(B_0+B)j'_x+B_{yy} a'_x={-}\mu\zeta'+\nu \nabla^2 \zeta'+F, \end{gather}$$

(3.7)$$\begin{gather}a'_t+U a'_x-(B_0+B)\psi'_x=\eta \nabla^2 a', \end{gather}$$

and the mean flow and field evolve according to

(3.8)$$\begin{gather} U_t-\overline{(\psi_x'\psi_y'-a_x'a_y')_y}={-}\mu U+\nu U_{yy}, \end{gather}$$

(3.9)$$\begin{gather}B_t-\overline{(a'\psi'_x)}_{yy}=\eta B_{yy}. \end{gather}$$

Equations (3.6)–(3.9) constitute the quasilinear approximation for the MHD equations. Its validity will be tested by numerical simulation of the full equations (3.1)–(3.3a,b) later in § 3.6. The term $a_x'a_y'$ is often referred to as the Maxwell stress, which has the opposite sign but a similar structure to the Reynolds stress $\psi _x'\psi _y'$. Intuitively, one might expect that the Maxwell stress is the mechanism by which a magnetic field can inhibit zonal flows, but we need to solve the actual zonostrophic instability problem to reach a concrete conclusion.

3.2. Basic state

For the basic state, the forcing $F$ generates waves with zero mean flow and perturbation field, $U=B=0$. With the forcing (2.6), we consider solutions of the form

(3.10a,b)\begin{equation} \psi'=\psi_1=\skew{3}{\hat}{\psi}_1(t)\exp (\textrm{i} kx)+\mathrm{c.c.},\quad a'=a_1=\hat{a}_1(t)\exp (\textrm{i} kx)+\mathrm{c.c.}, \end{equation}

and we seek the amplitudes $\hat {\psi }_1(t)$ and $\hat {a}_1(t)$. Substituting (2.6) and (3.10a,b) into (3.6) and (3.7) yields two ordinary differential equations in $t$, which we write in vector form:

(3.11) \begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix} \skew{3}{\hat}{\psi}_1\\ \hat{a}_1\end{matrix}\right)=\left(\begin{matrix} \displaystyle{\dfrac{\mathrm{i}\beta}{k}}-\mu-\nu k^2 & \mathrm{i}kB_0\\ \mathrm{i}kB_0 & -\eta k^2\end{matrix}\right)\left(\begin{matrix} \skew{3}{\hat}{\psi}_1\\ \hat{a}_1\end{matrix}\right)+\left(\begin{matrix} \displaystyle{-\dfrac{\sigma}{k^2}}\hat{\xi}\\ 0\end{matrix}\right). \end{equation}

Their solutions are

(3.12)$$\begin{align} \skew{3}{\hat}{\psi}_1 &={-}\varPsi_+\exp({\lambda_{1+}t})\int^t \exp({-\lambda_{1+}r})\hat{\xi}(r)\,\mathrm{d}r \nonumber\\ &\quad +\varPsi_-\exp({\lambda_{1-}t})\int^t \exp({-\lambda_{1-}r})\hat{\xi}(r)\,\mathrm{d}r, \end{align}$$

(3.13)$$\begin{align}\hat{a}_1 &={-}A\exp({\lambda_{1+}t})\int^t \exp({-\lambda_{1+}r})\hat{\xi}(r)\,\mathrm{d}r\nonumber\\ &\quad +A \exp({\lambda_{1-}t})\int^t \exp({-\lambda_{1-}r})\hat{\xi}(r)\,\mathrm{d}r, \end{align}$$

where

(3.14)\begin{equation} \left.\begin{gathered} \varPsi_{{\pm}}=\frac{\sigma}{2k^2 Q}\left(\frac{\mathrm{i}\beta}{k}- \mu-\nu k^2+\eta k^2 \pm Q\right),\quad A= \frac{\mathrm{i}\sigma B_0}{kQ}, \\ \lambda_{1\pm}=\frac{1}{2}\left(\frac{\mathrm{i}\beta}{k}-\mu-\nu k^2-\eta k^2\pm Q\right),\\ Q=\left[\left(\frac{\mathrm{i}\beta}{k}-\mu-\nu k^2+\eta k^2\right)^2-4k^2B_0^2\right]^{1/2}. \end{gathered}\right\} \end{equation}

Note that $\lambda _{1\pm }$ are the two eigenvalues of the matrix in (3.11). In the limit $B_0\rightarrow 0$,

(3.15a,b)\begin{equation} \lambda_{1+}\rightarrow \frac{\mathrm{i}\beta}{k}-\mu-\nu k^2=\lambda_1,\quad \lambda_{1-}\rightarrow -\eta k^2. \end{equation}

Hence, $\lambda _{1+}$ recovers the hydrodynamic eigenvalue of Rossby waves and $\lambda _{1-}$ is the rate of Ohmic damping of a magnetic field mode.

3.3. Instability analysis

Upon the basic state, we add the disturbances of zonostrophic instability, denoted by a subscript ‘2’:

(3.16a–d)\begin{align} \psi'=\psi_1+\psi_2 + \cdots ,\quad a'=a_1+a_2 + \cdots ,\quad U=U_2 + \cdots ,\quad B=B_2 + \cdots . \end{align}

Substituting (3.16a–d) into (3.6)–(3.9) and then linearising the ‘2’ components, we obtain for the fluctuating fields,

(3.17)\begin{align} & \zeta_{2,t}+\beta\psi_{2,x}-B_0 \nabla^2 a_{2,x}+\mu \zeta_2-\nu \nabla^2 \zeta_2 \nonumber\\ &\quad ={-}U_2 \zeta_{1,x}+U_{2,yy}\psi_{1,x}+B_2 \nabla^2 a_{1,x}-B_{2,yy}a_{1,x}, \end{align}

(3.18)$$\begin{align} &a_{2,t}-B_0\psi_{2,x}-\eta \nabla^2 a_2 ={-}U_2 a_{1,x}+B_2 \psi_{1,x}, \end{align}$$

(3.19)$$\begin{align} &\zeta_2 =\nabla^2 \psi_2, \end{align}$$

and for the mean fields,

(3.20)$$\begin{gather} U_{2,t}+\mu U_2-\nu U_{2,yy} =\overline{\psi_{1,x}\psi_{2,yy}+\psi_{2,x}\psi_{1,yy}-a_{1,x}a_{2,yy}-a_{2,x}a_{1,yy}}, \end{gather}$$

(3.21)$$\begin{gather}B_{2,t}-\eta B_{2,yy} =(\overline{a_1\psi_{2,x}+a_2\psi_{1,x}})_{yy}. \end{gather}$$

Note that here we are referring to $B_2$ as the ‘mean’ (magnetic) field for succinctness, and will continue to do so, though really it is only the perturbation component of the full mean field $B_0 + B_2(y,t) + \cdots$, with the mean always taken in $x$.

Similarly to the hydrodynamic case, we seek solutions in the form

(3.22) \begin{equation} \left.\begin{gathered} U_2 =\hat{U}(t) \exp (\textrm{i} my)+\mathrm{c.c.}, \\ \psi_2=\skew{3}{\hat}{\psi}_{21}(t) \exp({\mathrm{i}kx+\mathrm{i}my})+ \skew{3}{\hat}{\psi}_{22}(t) \exp({-\mathrm{i}kx+\mathrm{i}my})+\mathrm{c.c.} , \\ B_2 =\hat{B}(t) \exp (\textrm{i} kx)+\mathrm{c.c.}, \\ a_2=\hat{a}_{21}(t) \exp({\mathrm{i}kx+\mathrm{i}my})+\hat{a}_{22}(t) \exp({-\mathrm{i}kx+\mathrm{i}my})+\mathrm{c.c.} \end{gathered}\right\} \end{equation}

Substituting (3.22) into (3.17)–(3.21), we obtain the evolution equations for the fluctuating amplitudes,

(3.23)\begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix} \skew{3}{\hat}{\psi}_{21} \\ \hat{a}_{21} \end{matrix}\right)&=\left(\begin{matrix} \displaystyle\dfrac{\beta\mathrm{i}k}{m^2+k^2}-\mu-\nu(m^2+k^2) & \mathrm{i}kB_0 \\ \mathrm{i}kB_0 & -\eta(m^2+k^2) \end{matrix}\right) \left(\begin{matrix} \skew{3}{\hat}{\psi}_{21} \\ \hat{a}_{21} \end{matrix}\right) \nonumber\\ &\quad -\mathrm{i}k \left(\begin{matrix} \varLambda (\skew{3}{\hat}{\psi}_1\hat{U}-\hat{a}_1\hat{B}) \\ \hat{a}_{1}\hat{U}-\skew{3}{\hat}{\psi}_1\hat{B} \end{matrix}\right), \end{align}

(3.24)\begin{align} \frac{\mathrm{d}}{\mathrm{d}t}\left(\begin{matrix} \skew{3}{\hat}{\psi}_{22} \\ \hat{a}_{22} \end{matrix}\right)&=\left(\begin{matrix} \displaystyle-\dfrac{\beta\mathrm{i}k}{m^2+k^2}-\mu-\nu(m^2+k^2) & -\mathrm{i}kB_0 \\ -\mathrm{i}kB_0 & -\eta(m^2+k^2) \end{matrix}\right) \left(\begin{matrix} \skew{3}{\hat}{\psi}_{22} \\ \hat{a}_{22} \end{matrix}\right) \nonumber\\ &\quad +\mathrm{i}k \left(\begin{matrix} \varLambda (\skew{3}{\hat}{\psi}_1^*\hat{U}-\hat{a}_1^*\hat{B}) \\ \hat{a}_{1}^*\hat{U}-\skew{3}{\hat}{\psi}_1^*\hat{B} \end{matrix}\right), \end{align}

and for the mean flow and field,

(3.25)$$\begin{gather} \frac{\mathrm{d}\hat{U}}{\mathrm{d}t}+\mu \hat{U}+\nu m^2\hat{U} =\mathrm{i}k m^2 (\skew{3}{\hat}{\psi}_1^*\skew{3}{\hat}{\psi}_{21}-\hat{a}_1^*\hat{a}_{21})-\mathrm{i}km^2 (\skew{3}{\hat}{\psi}_1\skew{3}{\hat}{\psi}_{22}-\hat{a}_1\hat{a}_{22}), \end{gather}$$

(3.26)$$\begin{gather}\frac{\mathrm{d}\hat{B}}{\mathrm{d}t}+\eta m^2\hat{B} =\mathrm{i}km^2(\skew{3}{\hat}{\psi}_1^*\hat{a}_{21}-\hat{a}_1^*\skew{3}{\hat}{\psi}_{21})-\mathrm{i}km^2 (\skew{3}{\hat}{\psi}_1\hat{a}_{22}-\hat{a}_1\skew{3}{\hat}{\psi}_{22}). \end{gather}$$

The solutions of (3.23) are

(3.27) \begin{align} \skew{3}{\hat}{\psi}_{21} &= \int^t \left\{[-\varLambda D_+\skew{3}{\hat}{\psi}_1(\tau) +M\hat{a}_1(\tau)]\hat{U}(\tau)\right.\nonumber\\ &\left.\quad -[M\skew{3}{\hat}{\psi}_1(\tau)-\varLambda D_+\hat{a}_1(\tau)]\hat{B}(\tau)\right\} \exp[{\lambda_{2+}(t-\tau)}]\,\mathrm{d}\tau \nonumber\\ &\quad +\int^t \left\{[\varLambda D_-\skew{3}{\hat}{\psi}_1(\tau) - M\hat{a}_1(\tau)]\hat{U}(\tau)\right.\nonumber\\ &\quad \left. +\,[M\skew{3}{\hat}{\psi}_1(\tau)-\varLambda D_-\hat{a}_1(\tau)]\hat{B}(\tau)\right\} \exp[{\lambda_{2-}(t-\tau)}]\,\mathrm{d}\tau, \end{align}

(3.28) \begin{align} \hat{a}_{21} &= \int^t\left\{[\varLambda M\skew{3}{\hat}{\psi}_1(\tau) +D_-\hat{a}_1(\tau)] \hat{U}(\tau)\right.\nonumber\\ &\left.\quad -[D_-\skew{3}{\hat}{\psi}_1(\tau)+\varLambda M\hat{a}_1(\tau)]\hat{B}(\tau)\right\} \exp[{\lambda_{2+}(t-\tau)}]\,\mathrm{d}\tau \nonumber\\ &\quad +\int^t\left\{-[\varLambda M\skew{3}{\hat}{\psi}_1(\tau) +D_+\hat{a}_1(\tau)] \hat{U}(\tau)\right.\nonumber\\ &\left. \quad +[D_+\skew{3}{\hat}{\psi}_1(\tau)+\varLambda M\hat{a}_1(\tau)]\hat{B}(\tau)\right\} \exp[{\lambda_{2-}(t-\tau)}]\,\mathrm{d}\tau, \end{align}

where

(3.29) \begin{equation} \left.\begin{gathered} \lambda_{2\pm}=\frac{\mathrm{i}\varOmega_2-\mu-(\nu+\eta)k_2^2\pm Q_2}{2},\\ Q_2=\left[(\mathrm{i}\varOmega_2-\mu-\nu k_2^2+\eta k_2^2)^2-4k^2B_0^2\right]^{1/2}, \\ D_{{\pm}}=\frac{\mathrm{i}k(\mathrm{i}\varOmega_2-\mu-\nu k_2^2+\eta k_2^2\pm Q_2)}{2Q_2},\\ M=\frac{k^2B_0}{Q_2},\quad \varOmega_2=\frac{\beta k}{m^2+k^2},\quad k_2^2=m^2+k^2. \end{gathered}\right\} \end{equation}

For the solution of (3.24), we notice that the right-hand side of (3.23) is the complex conjugate of the right-hand side of (3.24), except that $\hat {U}$ and $\hat {B}$ remain the same. Therefore, we may find the solutions for the ‘22’ components by simply taking the complex conjugate of (3.27) and (3.28), and then replace all occurrences of $\hat {U}^*$ and $\hat {B}^*$ by $\hat {U}$ and $\hat {B}$.

We then substitute (3.27) and (3.28) into (3.25) and (3.26) to find a system of equations for the mean velocity and field:

(3.30) \begin{align} &\frac{\mathrm{d}\hat{U}}{\mathrm{d}t}+\mu\hat{U}+\nu m^2\hat{U}\nonumber\\ &\quad = \mathrm{i}m^2k\int^t\left\{[-\varLambda D_+\skew{3}{\hat}{\psi}_1^*(t) \skew{3}{\hat}{\psi}_1(\tau)+M\skew{3}{\hat}{\psi}_1^*(t)\hat{a}_1(\tau)\right. \nonumber\\ &\qquad -\varLambda M \hat a^*_1(t)\skew{3}{\hat}{\psi}_1(\tau)-D_-\hat{a}_1^*(t)\hat{a}_1(\tau)]\hat{U}(\tau) \nonumber\\ &\qquad +[{-}M \skew{3}{\hat}{\psi}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau)+\varLambda D_+\skew{3}{\hat}{\psi}_1^*(t) \hat{a}_1(\tau)+D_-\hat{a}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau) \nonumber\\ &\qquad \left. +\,\varLambda M \hat{a}_1^*(t) \hat{a}_1(\tau)]\hat{B}(\tau)\right\} \exp[{\lambda_{2+}(t-\tau)}]\,\mathrm{d}\tau \nonumber\\ &\qquad +\mathrm{i}m^2k\int^t \left\{[\varLambda D_-\skew{3}{\hat}{\psi}_1^*(t)\skew{3}{\hat}{\psi}_1 (\tau)-M\skew{3}{\hat}{\psi}_1^*(t)\hat{a}_1(\tau) \right. \nonumber\\ &\qquad +\varLambda M\hat a^*_1(t) \skew{3}{\hat}{\psi}_1(\tau)+D_+\hat{a}_1^*(t)\hat{a}_1(\tau)]\hat{U}(\tau) \nonumber\\ &\qquad +[M \skew{3}{\hat}{\psi}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau)-\varLambda D_-\skew{3}{\hat}{\psi}_1^*(t) \hat{a}_1(\tau)-D_+\hat{a}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau) \nonumber\\ &\left.\qquad -\varLambda M \hat{a}_1^*(t) \hat{a}_1(\tau)]\hat{B}(\tau)\right\} \exp[{\lambda_{2-}(t-\tau)}]\,\mathrm{d}\tau \nonumber\\ &\qquad + \mathrm{c.c.e.}\hat{U}.\hat{B}., \end{align}

(3.31) \begin{align} &\frac{\mathrm{d}\hat{B}}{\mathrm{d}t}+\eta m^2\hat{B}\nonumber\\ &\quad =\mathrm{i}m^2k\int^t\left\{[\varLambda M\skew{3}{\hat}{\psi}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau)+D_-\skew{3}{\hat}{\psi}_1^*(t) \hat{a}_1(\tau) \right.\nonumber\\ &\qquad +\varLambda D_+\hat{a}^*_1(t)\skew{3}{\hat}{\psi}_1(\tau)-M\hat{a}_1^*(t)\hat{a}_1(\tau)]\hat{U}(\tau) \nonumber\\ &\qquad -[D_-\skew{3}{\hat}{\psi}^*_1(t)\skew{3}{\hat}{\psi}_1(\tau)+\varLambda M \skew{3}{\hat}{\psi}^*_1(t) \hat{a}_1(\tau)-M\hat{a}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau) \nonumber\\ &\left.\qquad +\varLambda D_+\hat{a}^*_1(t)\hat{a}_1(\tau)] \hat{B}(\tau)\right\}\exp[{\lambda_{2+}(t-\tau)}]\,\mathrm{d}\tau \nonumber\\ &\qquad +\mathrm{i}m^2k\int^t\left\{[-\varLambda M\skew{3}{\hat}{\psi}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau)-D_+ \skew{3}{\hat}{\psi}_1^*(t)\hat{a}_1(\tau) \right.\nonumber\\ &\qquad -\varLambda D_-\hat{a}^*_1(t) \skew{3}{\hat}{\psi}_1(\tau)+M\hat{a}_1^*(t)\hat{a}_1(\tau)]\hat{U}(\tau) \nonumber\\ &\qquad + [D_+\skew{3}{\hat}{\psi}^*_1(t)\skew{3}{\hat}{\psi}_1(\tau)+\varLambda M \skew{3}{\hat}{\psi}^*_1(t) \hat{a}_1(\tau)-M\hat{a}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau) \nonumber\\ &\left.\qquad +\varLambda D_-\hat{a}^*_1(t) \hat{a}_1(\tau)]\hat{B}(\tau)\right\}\exp[{\lambda_{2-}(t-\tau)}]\,\mathrm{d}\tau \nonumber\\ &\qquad +\mathrm{c.c.e.}\hat{U}.\hat{B}. \end{align}

The notation ‘c.c.e.$\hat {U}$.$\hat {B}$’. means the complex conjugate of the previous terms except $\hat {U}$ and $\hat {B}$ remain unchanged (cf. ‘c.c.e.$s$.’ for (2.28) and (2.29)). There is a proliferation of terms in the MHD mean flow and field equations compared with the hydrodynamic case (cf. (2.22)), as we have doubled the number of fields and waves in the basic state ($\hat {\psi }_1$ and $\hat {a}_1$) and in the harmonics ($\psi _{21}$, $a_{21}$, $\psi _{22}$, $a_{22}$), and now have two mean fields ($\hat {U}$ and $\hat {B}$).

The next step is to take the expectation of (3.30) and (3.31), and to seek exponentially growing solutions for the mean flow and mean field expectations,

(3.32a,b) \begin{equation} \mathbb{E}[\hat{U}(t)]=\tilde{U} \exp (st), \quad \mathbb{E}[\hat{B}(t)]=\tilde{B} \exp (st). \end{equation}

As for hydrodynamic instability, we assume that the mean flow and the mean field are both statistically uncorrelated with the fundamental waves, leading to the separation of expectations following our partial ergodic assumption,

(3.33a,b)\begin{equation} \left.\begin{gathered} \mathbb{E}[\hat{U}(\tau)\skew{3}{\hat}{\psi}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau)]= \mathbb{E}[\hat{U}(\tau)]\mathbb{E}[\skew{3}{\hat}{\psi}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau)],\\ \mathbb{E}[\hat{B}(\tau)\skew{3}{\hat}{\psi}_1^*(t)\hat{a}_1(\tau)]= \mathbb{E}[\hat{B}(\tau)] \mathbb{E}[\skew{3}{\hat}{\psi}_1^*(t)\hat{a}_1(\tau)], \end{gathered}\right\} \end{equation}

and similarly for related terms. Equation (3.33a,b) again is weaker than the corresponding full ergodic assumption, namely

(3.34a,b)\begin{equation} \hat{U}= \mathbb{E}[\hat{U}],\quad \hat{B}= \mathbb{E}[\hat{B}], \end{equation}

that both the mean flow and field have no stochasticity (Constantinou & Parker Reference Constantinou and Parker2018).

The expectations of the fundamental waves $\hat {\psi }_1^*(t)\hat {\psi }_1(\tau )$, $\hat {\psi }_1^*(t)\hat {a}_1(\tau )$, etc. can now be computed as in (2.27), and then applying (3.32a,b) to the expectation of (3.30) and (3.31) provides the dispersion relation, with details deferred to Appendix C. The final result is

(3.35)\begin{equation} \left(\begin{matrix} (s+\mu+\nu m^2) \tilde{U} \\ (s+\eta m^2) \tilde{B} \end{matrix}\right)=\left(\begin{matrix} N_{UU} & N_{UB} \\ N_{BU} & N_{BB} \end{matrix}\right) \left(\begin{matrix} \tilde{U} \\ \tilde{B} \end{matrix}\right). \end{equation}

The terms $N_{UU}$, $N_{UB}$, $N_{BU}$ and $N_{BB}$, which are functions of $s$ and the parameters, with their expressions being given in (C9), have a clear physical meaning: $N_{UB}$ represents the feedback to the mean flow $\tilde {U}$ from the mean field $\tilde {B}$, etc. The condition for a non-trivial solution for $\tilde {U}$ and $\tilde {B}$ yields the dispersion relation which determines the eigenvalues $s$:

(3.36)\begin{equation} (s+\mu+\nu m^2-N_{UU})(s+\eta m^2-N_{BB})=N_{UB}N_{BU}. \end{equation}

When the basic state magnetic field is switched off, i.e. $B_0=0$, the coupling terms $N_{UB}$ and $N_{BU}$ are zero, and the two eigenvalues are given by

(3.37a,b)\begin{equation} s+\mu+\nu m^2=N_{UU},\quad s+\eta m^2 =N_{BB}. \end{equation}

Expression (3.37a) is the dispersion relation for the hydrodynamic zonostrophic instability (2.29), and (3.37b) is

(3.38)\begin{equation} s+\eta m^2+\frac{m^2\sigma^2}{2k^2(\mu+\nu k^2)} \left(\frac{1}{s+\mu+\nu k^2+\eta(m^2+k^2)+\mathrm{i}\beta/k}+\mathrm{c.c.e.}s.\right)=0. \end{equation}

It always gives a negative growth rate, representing the damping of the mean magnetic field by the flow. As in the hydrodynamic case, we undertake a rescaling to reduce the number of parameters by applying (2.30a–e) together with

(3.39a,b)\begin{equation} B_0=B_{0\star}\sigma^{{2}/{3}}k^{{-}1},\quad \eta=\eta_\star k^{{-}2}\sigma^{{2}/{3}}, \end{equation}

and then (3.36) can be expressed as

(3.40)\begin{equation} s_\star{=}s_\star(m_\star, \mu_\star, \nu_\star, \beta_\star, B_{0\star}, \eta_\star). \end{equation}

It is difficult to compare our dispersion relation directly with that of Constantinou & Parker (Reference Constantinou and Parker2018) because, as they state, their expression is complicated and uninformative, and so is not published in their paper. Their result is also for ring forcing, unlike ours for point forcing. Nonetheless, we will see that our numerical solution of the dispersion relation bears some similarity to theirs.

3.4. Asymptotic dispersion relations in limiting cases

The expression for the full dispersion relation, (3.36) or (3.40), is complicated, but in certain limits of the parameters, it is possible to obtain simpler expressions, as explained in Appendix C, that can provide deeper insights into the effect of the magnetic field and the mechanisms driving instability. These limits include strong magnetic diffusivity, and weak magnetic diffusivity with strong and weak magnetic fields. When discussing the applicability of these expressions, we assume $\mu _\star \ll \nu _\star \ll 1$, while $\beta _\star$ is of order of unity or larger, as previously.

3.4.1. The limit of large $\eta _\star$

In the limit $\eta _\star \rightarrow \infty$, the dispersion relation approximates as

(3.41)\begin{align} & s_\star{+}\mu_\star{+}\nu_\star m_\star^2 \nonumber\\ &\quad =\frac{m_\star^2 (1-m_\star^2)}{2\left(\mu_\star{+}\nu_\star{+} \dfrac{B_{0\star}^2}{\eta_\star}\right) \left\{[s_\star{+}2\mu_\star{+}\nu_\star(m_\star^2+2)](1 +m_\star^2)+ \dfrac{B_{0\star}^2}{\eta_\star}(2+m_\star^2)+ \mathrm{i}\beta_\star m_\star^2\right\}} \nonumber\\ &\qquad +\mathrm{c.c.e.}s . \end{align}

For finite $\eta _\star$, (3.41) is valid when $\eta _\star \gg \beta _\star$ and $B_{0\star }^2/\eta _\star \ll 1$. Equation (3.41) may be interpreted as the hydrodynamic expression (2.28) with $\lambda _1$ and $\lambda _{2}$ replaced by $\lambda _{1+}$ and $\lambda _{2+}$, which have the asymptotic expressions

(3.42a,b)\begin{align} \lambda_{1+}\sim \frac{\mathrm{i}\beta}{k}-\mu- \nu k^2-\frac{B_0^2}{\eta},\quad \lambda_{2+}\sim \frac{\mathrm{i}\beta k}{k^2+m^2}-\mu-\nu (m^2+k^2)-\frac{k^2 B_0^2}{\eta (m^2+k^2)}, \end{align}

in the large $\eta$ limit. Hence the action of the magnetic field is to enhance the damping of the waves, and thus to have a stabilising effect.

The solutions of $s_\star$ predicted by (3.41) are real, as in the hydrodynamic case. As $B_{0\star }$ increases, $s_\star$ will reduce and finally become negative, so that the instability is suppressed by the field. Equation (3.41) remains a sound approximation in this limit, as we will demonstrate in § 3.5 via comparison with numerical solutions of the full dispersion relation. Note that the right-hand side of (3.41) arises solely from the Reynolds stress, which is strongly damped by the magnetic field; the Maxwell stress itself is negligible in this regime.

3.4.2. The limit of small $\eta _\star$ and small $B_{0\star }$

In the limit of $B_{0\star },\eta _\star \rightarrow 0$ while $B_{0\star }^2/\eta _\star$ remains finite, the dispersion relation reduces to a different expression

(3.43)\begin{align} s_\star{+}\mu_\star{+}\nu_\star m_\star^2 &= \frac{m_\star^2(1-m_\star^2)}{2(\mu_\star{+}\nu_\star)} \frac{1}{(s_\star{+}2\mu_\star{+}2\nu_\star{+}\nu_\star m_\star^2) (1+m_\star^2)+\mathrm{i}\beta_\star m_\star^2}- \frac{m_\star^2B_{0\star}^2}{2\eta_\star \beta_\star^2 s_\star} \nonumber\\ &\quad +\mathrm{c.c.e.}s. \end{align}

This expression corresponds to a simplified mean-flow evolution equation,

(3.44) \begin{align} \frac{\mathrm{d}\hat{U}}{\mathrm{d}t}+\mu\hat{U}+\nu m^2\hat{U} & =mk^2\int^t\left\{\varLambda \skew{3}{\hat}{\psi}_1^*(t)\skew{3}{\hat}{\psi}_1(\tau) \exp[{\lambda_{2+}(t-\tau)}]\right. \nonumber\\ &\left.\quad -\hat{a}_1^*(t)\hat{a}_1(\tau) \exp[{\lambda_{2-}(t-\tau)}]+\mathrm{c.c.}\vphantom{\skew{3}{\hat}{\psi}_1^*}\right\}\hat{U}(\tau)\,\mathrm{d}\tau. \end{align}

The magnetic term in (3.43) therefore comes from the quadratic term $\hat {a}_1^*(t)\hat {a}_1(\tau )$, its negative sign exhibiting a stabilising effect of the magnetic field. From the original mean flow equation (3.8), we observe that $\hat {a}_1^*(t)\hat {a}_1(\tau )$ and $\hat {\psi }_1^*(t)\hat {\psi }_1(\tau )$ arise from the Maxwell and Reynolds stresses, respectively, and so the stabilising effect of the magnetic field can be interpreted as the Maxwell stress counteracting the Reynolds stress. The quadratic term from the magnetic field shares similarities with equation  (39) of Chen & Diamond (Reference Chen and Diamond2020), but our system incorporates more details of the time dependent waves, unlike in their system where the stochastic magnetic field is taken to be static.

For finite $\eta _\star$ and $B_{0\star }$, (3.43) is valid when $\eta _\star, B_{0\star }^2\ll \nu _\star$. The magnetic term raises the order of the equation from cubic in (2.31) to quartic, and with this we will shortly find some branches of complex solutions for $s_\star$. The presence of complex eigenvalues $s_\star$ has important implications for the stochasticity of the flow, as we will see in § 3.6 via numerical simulations. The expression (3.43) can be used to predict the field strength that makes $s_\star$ complex for all wavenumbers, but not to predict the field that makes $\operatorname {\mathrm {Re}} s_\star <0$ and thus suppresses the exponential growth. Such a magnetic field turns out to be very strong, and we will demonstrate this regime subsequently in § 3.4.3.

3.4.3. The limit of small $\eta _\star$ and large $B_{0\star }$

Finally we derive an asymptotic dispersion relation for a relatively strong magnetic field and weak magnetic diffusivity. In the limit of $\eta _\star \rightarrow 0$ and $B_{0\star }\rightarrow \infty$, we find

(3.45)\begin{equation} (s_\star{+}\mu_\star{+}\nu_\star m_\star^2)s_\star{=}N_{UB\star}N_{BU\star}, \end{equation}

with

(3.46)\begin{align} N_{UB\star} &= \frac{\mathrm{i}m_\star^2(m_\star^2+2)[\mathrm{i}\beta_\star{-} \nu_\star(m_\star^4+m_\star^2)]}{2B_{0\star}(\mu_\star{+}\nu_\star)(m_\star^2+1) [\nu_\star m_\star^4+(\mathrm{i}\beta_\star{+}2s_\star{+}2\mu_\star{+}3\nu_\star) m_\star^2+2(s_\star{+}\mu_\star{+}\nu_\star)]} \nonumber\\ &\quad +\mathrm{c.c.e.}s., \end{align}

(3.47)\begin{align} N_{BU\star} &= \frac{-\mathrm{i}m_\star^2[2\nu_\star (m_\star^4+2m_\star^2+1)+ 2\mu_\star(m_\star^2+1)+(2s_\star{+}\mathrm{i}\beta_\star )m_\star^2]}{2B_{0\star} (\mu_\star{+}\nu_\star)[\nu_\star m_\star^4+(\mathrm{i}\beta_\star{+}2s_\star{+} 2\mu_\star{+}3\nu_\star)m_\star^2+2(s_\star{+}\mu_\star{+}\nu_\star)]} \nonumber\\ &\quad +\mathrm{c.c.e.}s. \end{align}

Here $(N_{UB\star }, N_{BU\star })=\sigma ^{-({2}/{3})}(N_{UB}, N_{BU})$ (see (3.36)), while $N_{UU}$ and $N_{BB}$ may be neglected at leading order. The relation (3.45) has some distinctive properties. In the previous two asymptotic dispersion relations (3.41) and (3.43), only $N_{UU}$ plays a role (corresponding also to (3.37a,b)), hence the hydrodynamic instability still has a major contribution, although it is modified by the magnetic field. On the contrary, in (3.45), $N_{UU}$ becomes negligible and $N_{UB}$ and $N_{BU}$ are the leading-order terms. This means the hydrodynamic instability is suppressed by the strong magnetic field, and the interaction between the mean flow and mean field is the dominant effect, which can yield zonostrophic instability. Another distinguishing feature of this dispersion relation is that $\eta _\star$ does not appear in (3.45)–(3.47). Hence when the magnetic diffusivity is weak enough, it no longer affects the instability.

For finite $B_{0\star }$ and $\eta _\star$, (3.45) is valid when $B_{0\star }$ is at order $\nu _\star ^{-1}$ or larger and $\eta _\star$ is much smaller than $\nu _\star$. The eigenvalue $s_\star$ in this regime is generally complex. The relation (3.45) can be used to predict the transition from $\operatorname {\mathrm {Re}} s_\star >0$ to $\operatorname {\mathrm {Re}} s_\star <0$ as $B_{0\star }$ increases, and so provide insights about the instability threshold.

In the next section, we will compare the predictions of the asymptotic dispersion relations with numerical solutions. We will also derive scaling laws for instability thresholds.

3.6. Numerical simulation

We now perform numerical simulations for the MHD flow governed by (3.1) and (3.2). We use a pseudospectral method, as described for the hydrodynamic case, except that for weak diffusivity $\eta =10^{-5}$ corresponding to $\eta _\star =0.019$, we use a higher resolution of $512\times 512$ mesh points. We use the same initial condition $\psi =10^{-4}\cos (my)$ for the flow, and there is no magnetic perturbation, $B=0$ initially; only the uniform background field $B_0$ is present. We will consider various values of $B_0$ and $\eta$, as studied in figures 5–7, and use the same values for the other parameters, namely $k=16$, $m=5$, $\sigma =0.05$, $\beta =5$, $\mu =0$, $\nu =10^{-4}$, corresponding to $\beta _\star =2.3$, $\nu _\star =0.19$, $m_\star =0.31$, as we did in the purely hydrodynamic simulation in § 2.5.

We first study the case of relatively high magnetic diffusivity $\eta =10^{-2}$, corresponding to the dispersion relations shown in figure 5. The evolution of the vorticity $\zeta$ and current density $j$ are shown in figure 10(a,b), and the evolution of the mean flow $U$ and mean field $B$ are shown in figure 10(c–f). The behaviour of $\zeta$ and $U$ is similar to the hydrodynamic case – zonostrophic instability causes exponential growth of the zonal flow, forming two zonal jets in opposite directions. The zonal jets shear the vorticity, causing sinuous winding which becomes stronger over time. The exponential growth of $U_{rms}$ in figure 10(d) agrees very well with that of the expectation predicted by theory. The spatial structure of $j$ in figure 10(b) behaves in a similar way to that of $\zeta$ in figure 10(a), but has a much weaker amplitude. The mean field $B$ grows exponentially during the exponential growth of the zonal flow in figure 10( f), but then falls back to very small values. It thus appears that the field is largely controlled by the flow.

Figure 10. Numerical simulation for the MHD flow for $B_0=0.01$, $\eta =10^{-2}$, $\beta =5$, $k=16$, $\nu =10^{-4}$, $\mu =0$, $\sigma =0.05$, $m=5$, corresponding to $\beta _\star =2.3$, $\nu _\star =0.19$, $B_{0\star }=1.2$, $\eta _\star =19$, $m_\star =0.31$. Panels (a,b) show the evolution of vorticity $\zeta$ and current density $j$; panels (c,e) are the Hovmöller diagram for the mean flow $U(y,t)$ and field $B(y,t)$ and panels (d,f) show the r.m.s. value of $U$ and $B$, respectively; solid lines are the results of the numerical simulation, and the red straight lines are the theoretical prediction, with growth rate $s=4.35\times 10^{-2}$.

To further explore the stochastic behaviour of the MHD system, we consider three different magnetic field strengths $B_0=0.01, 0.03, 0.05$ at the same $\eta =10^{-2}$, and run 10 simulations for each case. The numerical results for $U_{rms}$ are shown in figure 11(a–c). We also plot a thick straight line to indicate the exponential growth of the expectation, as predicted by theory. The value of the growth rate $s$ is indicated in the title of each plot. Although there is significant variability, the theoretical prediction for the expectation generally agrees with the growth or decay in each realisation, thus providing support for the full ergodic assumption (3.34a,b). In figure 11(d), we show the ensemble averages and the theoretical prediction. There is good agreement for all three cases, confirming the theory. The stabilising effect of the magnetic field is also clearly demonstrated.

Figure 11. Temporal evolution of $U_{rms}$ for $\eta =10^{-2}$ ($\eta _\star =19$) and (a) $B_0=0.01$ $(B_{0\star }=1.2)$, (b) $B_0=0.03$, (c) $B_0=0.05$. The other parameters are the same as in figure 10. Solid lines represent 10 realisations of the numerical simulation of the stochastic flow, and the thick red lines give the theoretical prediction. In (d), the ensemble average of the 10 realisations for each $B_0$ is shown together with the theoretical prediction plotted in straight lines.

Next, we consider the zonostrophic instability at a moderate diffusivity $\eta =10^{-4}$; results from the dispersion relation for this case have been set out in figure 6. In figure 12, we present the evolution of $j$, $\eta$, $U$ and $B$ for one realisation with $B_0=0.01$. Zonal flows again emerge as a result of zonostrophic instability, but the weaker magnetic diffusivity results in a stronger influence of the magnetic field. For example, the strengths of $j$ and $B$ are now of the same order as $\zeta$ and $U$, respectively, indicating similar importance of the flow and the field. The spatial pattern of $j$ in figure 12(b) is characterised by elongated thin filaments, somewhat different from the vortices seen for $\zeta$ in figure 12(a). The growth of $U_{rms}$ in figure 12(d) is somewhat faster than that predicted for the expectation by theory, and the agreement is now only qualitative.

Figure 12. Numerical simulation for the MHD flow for $B_0=0.01$, $\eta =10^{-4}$ corresponding to $B_{0\star }=1.2$, $\eta _\star =0.19$, with the theoretical growth rate $s=3.31\times 10^{-2}$. The other parameters and the contents of the panels are the same as in figure 10.

Following the typical realisation in figure 12 for a moderate diffusivity $\eta =10^{-4}$ with $B=0.01$, we have undertaken a series of runs. Figure 13 shows 10 realisations for each of the four magnetic field strengths $B_0=5\times 10^{-3}, 10^{-2}, 3\times 10^{-2}$ and $10^{-1}$ in figure 13(a–d); the theoretical growth rate of the expectation is plotted using thick straight lines. The ensemble average for each field strength is shown in figure 13(e) and also compared with theory. The growth rates $s$ for $B_0=3\times 10^{-2}$ and $10^{-1}$ in figure 13(c,d) are complex, and we only plot the exponential growth or decay from the real part of $s$; we will elaborate more on this issue shortly. Note that the wavenumber is fixed to be $m=5$ for different $B_0$, hence the instability property of a particular $(B_0,\eta )$ simulation set here could be slightly different from that in figure 8 which is based on the most unstable wavenumber. For very weak magnetic field $B_0=5\times 10^{-3}$ in figure 13(a), each realisation has a growth rate reasonably close to that the theory predicts for the expectation, with good agreement for the ensemble average in figure 13(e), similar to the hydrodynamic case. As $B_0$ increases to $10^{-2}$ in figure 13(b), for realisations where exponential growth is prominent, the theoretical result still captures the behaviour fairly well, resulting in good agreement again for the average in figure 13(e). However, there are also many realisations in figure 13(b) where the zonal flow does not grow (up to the length $t=200$ of our simulations), and the system now shows an increased degree of randomness.

Figure 13. Temporal evolution of $U_{rms}$ at $\eta =10^{-4}$ ($\eta _\star =0.19$) and (a) $B_0=5\times 10^{-3}$, (b) $B_0=10^{-2}$ ($B_{0\star }=1.2$), (c) $B_0=3\times 10^{-2}$, (d) $B_0=10^{-1}$. The other parameters are the same as in figure 10. Solid lines represent 10 realisations of the numerical simulation of the stochastic flow, and the thick red lines represent the prediction of the zonostrophic instability. In (e), the ensemble average of the 10 realisations for each $B_0$ is shown by solid lines, compared with the red straight lines giving theoretical growth rates.

At a stronger background field $B_0=3\times 10^{-2}$, where $s=3.2\times 10^{-3}+6.5\times 10^{-3}\mathrm {i}$ becomes complex, individual realisations become more chaotic, in figure 13(c): the growing realisations often have much larger growth rates than the theory for the expectation, while the decaying ones may reach very small amplitudes. Their ensemble average in figure 13(e) also shows a poorer agreement with the expectation. Finally, at the largest $B_0=0.1$, the theory predicts that the expectation should decay; individual realisations may again behave differently, but these appear less chaotic in nature. At earlier times, the ensemble average has good agreement with the expectation, whilst at later times, occasional growth of some realisations make the ensemble average diverge from the expectation. In summary, at this lower value $\eta =10^{-4}$ of the magnetic diffusivity, zonal flows of individual realisations have a higher degree of randomness and their growth may diverge from that of the expectation as predicted by theory. The full ergodic assumption for the mean flow and field stated in (3.34a,b) is therefore very questionable here, and in fact does not appear to operate in any meaningful, qualitative, way. However, our instability analysis does not use the full ergodic assumption in the form (3.34a,b). Rather, we used the partial ergodic assumption (3.33a,b) which allows variation of individual realisations from the expectation, and is thus a better approximation in this situation. The use of the partial ergodic assumption is supported by the fairly good agreement between the expectation from the theory and the ensemble average of the simulations.

In figure 13(c,d), we also observe prominent high-frequency oscillations of the zonal flow. We find their frequencies are very well described by

(3.56)\begin{equation} \omega=2k B_0, \end{equation}

i.e. the zonal flow oscillates at twice the characteristic Alfvén wave frequency. We have confirmed that this frequency is not the frequency of the magneto-Rossby waves of either the basic state or the disturbances, represented by the imaginary parts of $\lambda _{1\pm }$ and $\lambda _{2\pm }$ (cf. (3.14) and (3.29)), respectively. Nor is it related to the imaginary part of $s$ whose physical meaning as yet remains obscure. Our understanding of this oscillation remains limited. From the appearance of (3.56), the oscillation results from the coupling between the forcing with wavenumber $k$ and the basic magnetic field $B_0$, but we do not have a quantitative argument to describe this. We have confirmed that this oscillation also appears in the linearised mean flows described by (3.30) and (3.31), but not in their expectations (C3) and (C4), thus it seems that the flow stochasticity is a key element to excite this oscillation.

Finally, we explore the case of very weak magnetic diffusivity $\eta =10^{-5}$, corresponding to the zonostrophic instability studied in figure 7. In figure 14, we show a realisation for $B_0=0.01$. The weak magnetic diffusivity renders very fine filaments in the spatial pattern of the current $j$, which also influences the pattern of $\zeta$. The exponential growth of $U_{rms}$ and $B_{rms}$ is much faster than that predicted by the theory for the expectation, indicating a high degree of stochasticity. The oscillations of the mean flow and the mean field with frequency twice the Alfvén frequency become more prominent. A key point we emphasise is that as the theoretical growth rate of zonostrophic instability is $s=0.016+0.039\mathrm {i}$, this case falls into the regime of ‘unstable with complex growth rates’ in figure 8. Hence, in contrast to Tobias et al. (Reference Tobias, Diamond and Hughes2007) who found this region to be stable, we have a concrete example of zonostrophic instability taking place and generating zonal flows. The data of the white noise for this realisation has been documented and is available online (see the data access statement at the end of the paper).

Figure 14. Numerical simulation for the MHD flow for $B_0=0.01$, $\eta =10^{-5}$, corresponding to $B_{0\star }=1.2$, $\eta _\star =0.019$, with theoretical growth rate $s=1.60\times 10^{-2}+3.94\times 10^{-2}\mathrm {i}$. The other parameters and panels are the same as in figure 10. In (d,f), only the real part of $s$ has been used to plot the theoretical growth.

We then show 10 realisations at various magnetic field strengths $B_0=2.5\times 10^{-3}, 10^{-2}, 2\times 10^{-2}$ and $5\times 10^{-2}$ with $\eta =10^{-5}$ fixed in figure 15. The behaviour shows many similarities to figure 13, but our main objective is to confirm the statistical behaviour of modes which are unstable with complex growth rates. At the low diffusivity of $\eta =10^{-5}$, a wide range of field strengths $B_0$ fall into this category, and as shown in figure 7, unstable complex modes are dominant for stronger magnetic fields. In figure 15, the cases of $B_0=10^{-2}$ and $2\times 10^{-2}$ in figure 15(b,c) are unstable with complex growth rates. As we see, the realisations for $B_0=10^{-2}$ and $2\times 10^{-2}$ are highly chaotic. The theory developed predicts an exponential growth of the expectation, but individual realisations either have much faster growth or decay significantly; very few evolve as theory predicts. The high frequency oscillations with frequency twice the Alfvén frequency also become very strong, making the evolution even more disorganised. There appear to be more decaying realisations for $B_0=2\times 10^{-2}$ than for $B_0=10^{-2}$, an indication of the stabilising effect of the field. The agreement between the ensemble average over the 10 realisations simulated and the growth rate of the expectation shown in figure 15(e) for $B_0=10^{-2}$ and $2\times 10^{-2}$ is adequate, but not as good as the case of $B_0=2.5\times 10^{-3}$ where the growth rate is real and positive. At an even stronger field $B_0=5\times 10^{-2}$ when $s_{r}$ becomes negative, individual realisations remain highly chaotic, but there is good agreement between their ensemble average and the decay predicted by the theory at the early stage. Hence, we have confirmed that the modes that are unstable with complex growth rates are highly stochastic; for an individual realisation, zonal flow may or may not emerge, and the growth rate predicted for the expectation of the mean fields has little relevance. While the full ergodic assumption (3.34a,b) clearly does not work, the partial ergodic assumption (3.33a,b) is also put into doubt given the theoretical prediction does not agree well with the numerical ensemble average. We could also question the validity of the quasilinear approximation (3.6)–(3.9), which is a fundamental assumption employed by the theory.

Figure 15. Temporal evolution of $U_{rms}$ for $\eta =10^{-5}$ ($\eta _\star =0.019$) and (a) $B_0=2.5\times 10^{-3}$, (b) $B_0=10^{-2}$ ($B_{0\star }=1.2$), (c) $B_0=2\times 10^{-2}$, (d) $B_0=5\times 10^{-2}$. Other information is the same as in figure 13.

Our conclusion seems different to that of Tobias et al. (Reference Tobias, Diamond and Hughes2007) who suggest that the region of ‘unstable with complex growth rates’ in figure 8 should have no zonal flow formation. We think the reason for this disparity may lie in the forcing: our simulations indicate that different realisations yield wildly different evolution of the mean flow, so that any single simulation is unrepresentative. It is also possible that the spatial structure of the forcing could make a difference: our forcing only has one wavenumber in the $x$-direction, while Tobias et al. (Reference Tobias, Diamond and Hughes2007) used a range of wavenumbers in both the $x$ and $y$-directions. We leave the issue of more realistic forcing for further investigation.

Despite our investigations of the modes with complex growth rates, there are still many aspects that we do not fully understand, in particular the physical meaning of the imaginary part of $s$. A straightforward interpretation is that it represents oscillation of the expectation, but oscillation is more a property of waves than mean flows. Indeed, the ensemble mean from numerical simulations does not pass through zero as the theory would otherwise predict. Mainly for this reason, we have not included the imaginary part of $s$ in comparisons with the numerical simulations. On the other hand, when $s$ is complex with a positive real part, the agreement between the theory and the simulations is not so good, suggesting that the imaginary part of $s$ has a role that we do not yet understand. Linked with this is the observation that when $s$ is complex, runs show highly disorganised evolution of mean flow and field; we do not know if these are related, nor can we currently explore this in depth given the expense of computing a sufficiently large ensemble. Since complex growth rates can also arise for hydrodynamic zonostrophic instability with a more complicated forcing structure (Ruiz et al. Reference Ruiz, Parker, Shi and Dodin2016), we suspect similar phenomena could take place in these non-magnetic systems. Studying this further, perhaps via a Fokker–Planck equation, is a topic for future research.