On statistical zonostrophic instability and the effect of magnetic fields

Abstract

Zonal flows are mean flows in the east–west direction, which are ubiquitous on planets, and can be formed through ‘zonostrophic instability’: within turbulence or random waves, a weak large-scale zonal flow can grow exponentially to become prominent. In this paper, we study the statistical behaviour of the zonostrophic instability and the effect of magnetic fields. We use a stochastic white noise forcing to drive random waves, and study the growth of a mean flow in this random system. The dispersion relation for the growth rate of the expectation of the mean flow is derived, and properties of the instability are discussed. In the limits of weak and strong magnetic diffusivity, the dispersion relation reduces to manageable expressions, which provide clear insights into the effect of the magnetic field and scaling laws for the threshold of instability. The magnetic field mainly plays a stabilising role and thus impedes the formation of the zonal flow, but under certain conditions it can also have destabilising effects. Numerical simulation of the stochastic flow is performed to confirm the theory. Results indicate that the magnetic field can significantly increase the randomness of the zonal flow. It is found that the zonal flow of an individual realisation may behave very differently from the expectation. For weak magnetic diffusivity and moderate magnetic field strengths, this leads to considerable variation of the outcome, that is whether zonostrophic instability takes place or not in individual realisations.

1. Introduction

Zonal flows are mean flows in the east–west direction, and are commonly found on Earth and other planets. They are perhaps most prominent on Jupiter, where belts of strong zonal jets are the main visible feature of its surface. Numerous studies have been undertaken to understand the properties of zonal flows. A review of this topic can be found in the recent book by Galperin & Read (2019); here we summarise the representative literature most relevant to our study.

Rhines (1975) identified that zonal jets have a scale associated with the wavenumber $k_\beta = \sqrt {\beta /2U}$, where $U$ is the root-mean square (r.m.s.) zonal velocity and $\beta$ is the gradient of the Coriolis parameter. This scale is now known as the Rhines scale. Above this length scale, the inverse cascade of turbulence is suppressed by the $\beta$-effect, and turbulence transfers energy to the zonal jets. Williams (1978) performed numerical simulations that reproduce the zonal flows on both Earth and Jupiter. From the conditions under which zonal flows emerge, he concluded that the $\beta$-effect and forcing are the two elements essential for the formation of zonal jets. Vallis & Maltrud (1993) identified the important scales of turbulence in the zonal jets, and also showed that in addition to the $\beta$-effect, topography can also generate zonal flows. Smith (2004) and Scott & Polvani (2007) considered quasigeostrophic flows and found that the effect of introducing a finite deformation radius is to increase the level of $\beta$ required for zonal jets to form. When the deformation radius is small enough, zonal jets are confined near the equator. Galperin et al. (2006) studied the energy spectrum of turbulent flows with strong zonal jets, and found that they exhibit a $-5$ power law, in what they termed a ‘zonostrophic regime’, and in contrast to the celebrated Kolmogorov $-\frac {5}{3}$ power law of isotropic turbulence. Dritschel & McIntyre (2008) proposed a mechanism for the formation of zonal jets: the mixing of potential vorticity results in staircases in the potential vorticity profile.

To further understand mechanisms for jet formation, Farrell & Ioannou (2003,  2007) established a compact system that combines the evolution of turbulence and mean flow. They coined the term ‘stochastic structural stability’, and showed that the state with zero mean flow can be unstable, leading to the formation of zonal jets. Srinivasan & Young (2012) advanced the theory by undertaking an analytical study of the instability problem and deriving an explicit dispersion relation, also introducing the term ‘zonostrophic instability’. Parker & Krommes (2013,  2014) incorporated weak nonlinearity and derived a Ginzburg–Landau equation for the zonal-flow amplitude, which they used to model the generation of zonal jets in terms of pattern formation.

In astrophysical contexts, it is necessary to consider the effect of magnetic fields on the flow and study the magnetohydrodynamic (MHD) problem. Diamond et al. (2005) have given a comprehensive review of zonal flows in plasmas, where the magnetic field plays a key role. Recent studies have also revealed the important link between zonal flows and cross-helicity (Heinonen et al. 2023), which is a central topic of MHD. Zonal flows may also be important for the structure of the solar tachocline, which is generally believed to be the source of the Sun's magnetic field (for a comprehensive review of the solar tachocline, see the monograph edited by Hughes, Rosner & Weiss (2007)). In this vein, Tobias, Diamond & Hughes (2007); Tobias, Dagon & Marston (2011) have undertaken numerical simulation of MHD $\beta$-plane flow in Cartesian and spherical geometry. Their results indicate that even a weak magnetic field can suppress the formation of zonal jets, and this may explain the lack of observations so far of zonal jets in the Sun, where strong magnetic fields interact with plasma flows. In terms of theory, Constantinou & Parker (2018) presented an MHD zonostrophic stability analysis using the method of Srinivasan & Young (2012). Their results confirmed that the magnetic field indeed reduces the growth rate of the zonostrophic instability. Durston & Gilbert (2016) and Algatheem, Gilbert & Hillier (2023) considered MHD zonostrophic instability for large-scale shear flows and for deterministic Kolmogorov flow, respectively. Analytical approximations can be made by exploiting scale separation between large-scale flow and small-scale waves, and these indicate that the magnetic field can inhibit the hydrodynamic zonostrophic instability, but it can also generate a new branch of unstable modes for such shear flows.

While it appears that a magnetic field can suppress zonostrophic instability and zonal flows, the underlying physical mechanism requires further exploration. To this end, Constantinou & Parker (2018) and Parker & Constantinou (2019) considered magnetic fields subject to the tilting present in shear flows. They showed that the Maxwell stress imposed by the magnetic field opposes the Reynolds stress, and therefore inhibits the formation of zonal flows. Chen & Diamond (2020) and Chen et al. (2021) considered turbulent magnetic fields that are strong and highly disordered. Via averaging the system over an intermediate scale, they showed that the resulting Maxwell stress again inhibits the growth of any zonal flow.

In the present paper, we will study zonostrophic instability with an emphasis on its statistical properties and the effect of magnetic fields. In particular, we consider the growth of weak mean flows in a system of random waves driven by white noise. We derive a dispersion relation for exponential growth of the expectation of the mean flow, and analyse its properties in detail. We also compare the theory with numerical simulations of the stochastic system, to assess its validity and the assumptions made. Our method for the instability analysis differs from previous studies which are mainly based on spatiotemporal correlation functions, e.g. Farrell & Ioannou (2003,  2007), Srinivasan & Young (2012) and Constantinou & Parker (2018). The use of spatiotemporal correlation functions works well for hydrodynamic zonostrophic instability, but becomes complicated when a magnetic field is added. For example, Constantinou & Parker (2018) describe the calculations leading to the dispersion relation for the MHD zonostrophic instability as ‘complicated and unilluminating’, and they do not document its full expression. This is mainly because the presence of the magnetic field introduces several new terms that are quadratic in the fields: for example, one has to include the correlation between the velocity and the magnetic field to close the system. Our study, on the other hand, directly analyses random waves and minimises these complications. Our derivation only involves the temporal correlations of the random waves, which significantly simplifies the mathematics. We are able derive a compact dispersion relation in certain limiting cases of the parameters, which provides simple scaling laws for the threshold of instability.

An outcome of our simplified dispersion relation is a clear explanation for the mechanism by which the magnetic field affects zonostrophic instability. The physical explanation of Constantinou & Parker (2018), Parker & Constantinou (2019), Chen & Diamond (2020) and Chen et al. (2021) is based on the form of the magnetic terms in the mean-flow equation. Our interpretation, on the other hand, is based on the dispersion relation for zonostrophic instability, in which the effect of the magnetic field on the growth of the mean flow is evident. Also, we do not assume a priori the time and length scales of the flow but solve the full mathematical problem under the quasilinear approximation, allowing a comprehensive survey of parameter space. With this, we observe further effects that can arise in the presence of a magnetic field in addition to the Maxwell stress discussed by previous authors: under certain conditions, the magnetic field can change the effective viscosity of the flow, or affect the instability through the interaction between the mean flow and the mean field.

In many previous studies of zonostrophic instability, an ergodic assumption has been employed: namely that the zonal mean velocity in each realisation is the same as in the ensemble average, in other words the zonal flow itself has no stochasticity. This ergodic assumption has been widely used for analyses of mean flows in various problems, for example, zonostrophic instability (Srinivasan & Young 2012; Constantinou & Parker 2018), the structural stability of turbulent jets (Farrell & Ioannou 2003,  2007) and pattern formation of zonal flows (Parker & Krommes 2013,  2014; Ruiz et al. 2016). Its validity has not been established in general, though Bouchet, Nardini & Tangarife (2013) provide a systematic argument based on separation of time scales between fluctuating and mean fields, making use of a Fokker–Planck formulation. More recently, Allawala, Tobias & Marston (2020) have shown that direct statistical simulation using an ensemble average and a zonal average may not produce exactly the same results, making it reasonable to question the ergodic assumption. In the present study, we will use numerical simulations to examine the ergodic assumption in the context of zonostrophic instability with and without a magnetic field. We run multiple simulations to see if individual realisations of the mean flow behave similarly, or not.

Another issue that arises for MHD zonostrophic instability concerns modes with complex growth rates. Constantinou & Parker (2018) showed that in the presence of a magnetic field, zonostrophic modes can possess complex growth rates, a phenomenon uncommon for pure hydrodynamic zonostrophic instability (e.g. Srinivasan & Young (2012), but see Ruiz et al. (2016)). The nature of these modes remains poorly understood: the real parts of the growth rates can be positive, but these modes lie in the parameter regime where Tobias et al. (2007) reported no zonal flow formation, thus whether zonostrophic instability takes place for these modes remains a question. We will undertake numerical simulations to explore whether these modes are indeed unstable and give rise to zonal flows.

The organisation of the paper is as follows. In § 2, we analyse purely hydrodynamic zonostrophic instability. Although the hydrodynamic problem is relatively well understood, we use it as a basic example to establish our method of analysis, and to lay the foundation for the MHD problem. We present the governing equations and the quasilinear approximation in § 2.1. Then we establish the ‘basic state’ in § 2.2, which is a forced–dissipative system without mean flow. The stability of this state is studied in § 2.3, and a dispersion relation is derived for the growth rate of the zonal mean flow. Its properties are briefly discussed in § 2.4. Finally in § 2.5, numerical simulations are presented to verify the theory. The analysis for the MHD zonostrophic instability is undertaken in § 3 following the same methodology, but it is a more complicated and interesting problem. In particular, the full dispersion relation derived in § 3.3 is a rather unwieldy expression; however, we derive simplified dispersion relations in § 3.4 in certain asymptotic limits of the parameters. In the discussion of the dispersion relation in § 3.5, we give a comprehensive investigation of the effect of the magnetic field strength and magnetic diffusivity, and compare this with previous studies. In numerical simulations in § 3.6, we demonstrate the strong impact of the magnetic field on the stochasticity of the evolving zonal flows. Concluding remarks are given in § 4.

2. Hydrodynamic zonostrophic instability

In this section, we present the analysis for the hydrodynamic zonostrophic instability. We give a straightforward derivation of a dispersion relation for the unstable growth rate of the zonal flow; this takes a relatively simple form and is equivalent to that obtained by Srinivasan & Young (2012). We do this to set out the methodology for the magnetic zonostrophic instability that will be studied in the next section.

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 (2003,  2007). It also has similarities with the approach of Srinivasan & Young (2012), 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 2003,  2007; Srinivasan & Young 2012). 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 2012; Farrell & Ioannou 2003; Marston, Conover & Schneider 2008). 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 (2012) 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 2001; Durston & Gilbert 2016). 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.4. Results and discussion

We now briefly discuss the properties of the dispersion relation (2.29) or (2.31), which governs hydrodynamic zonostrophic instability. First, we observe that our dispersion relation would exactly result from the analysis in Srinivasan & Young (2012), if we apply their derivation to our forcing. Although their discussion focused on a ring forcing, their equation (C16) is a dispersion relation for general forms of forcing, and in our case results in (2.29). We give the details in Appendix B, and note that the agreement is not trivial, since we followed a very different derivation.

Focusing on geophysical applications, Srinivasan & Young (2012) studied the effects of the drag coefficient $\mu$ in detail, and mostly set the viscosity $\nu$ to be zero. Our present study focuses on astrophysical applications, and we will mainly consider the situation where the drag $\mu$ is zero and study the effect of the viscosity $\nu$. We will choose $k=16$ and $\sigma =0.05$ as the scale and amplitude of the forcing, respectively. Tobias et al. (2007) considered wavenumbers with $14< k<20$, which includes our $k=16$. The small amplitude $\sigma =0.05$ is intended to generate weak turbulence where the quasilinear approximation is expected to be valid (cf. Srinivasan & Young 2012), and it is also of the same order as that used by Constantinou & Parker (2018).

In figure 1, we show the (real) growth rate $s$ versus vertical wavenumber $m$ at different $\nu$ and $\beta$ in figure 1(a,b) respectively, calculated from the dispersion relation (2.29). One feature of the dispersion relation is that as the drag coefficient $\mu$ is zero the growth rate $s$ approaches zero as $m\rightarrow 0$. This contrasts with the case $\mu >0$ studied by Srinivasan & Young (2012), where $s$ is negative at $m=0$. As seen in figure 1(a), increasing the viscosity reduces the growth rate and finally suppresses instability, as one would expect. In figure 1(b), we see that increasing $\beta$ also reduces instability (Srinivasan & Young 2012), but the growth rates at smaller wavenumbers are unaffected, as is evident from the dispersion relation (2.29). Note that zonostrophic instability continues to exist even for $\beta =0$. Even though there is now no background vorticity gradient present to give a preferred direction in the system, the basic state of the fluid system remains anisotropic because of the forcing we employ in (2.6). As is known, zonostrophic instability cannot occur in fluid systems that are isotropically forced, without some further mechanism to break symmetry, such as a $\beta$-effect or magnetic field (Srinivasan & Young 2012; Bakas & Ioannou 2013).

Figure 1. Dispersion relation for the hydrodynamic instability with $k=16$, $\mu =0$, $\sigma =0.05$. (a) Plots of the growth rates at different $\nu$ with $\beta =5$, and (b) plots of the growth rates at different $\beta$ with $\nu =10^{-4}$. Here $\nu =10^{-3}$ corresponds to $\nu _\star =1.9$, and $\beta =10$ corresponds to $\beta _\star =4.6$.

From figure 1, we see that neutral stability $s\to 0$ occurs in the limit $m\rightarrow 0$. Applying $m\rightarrow 0$ to the dispersion relation (2.29) or (2.31) with $\mu =0$, we find that the condition for zonostrophic instability is

(2.32)\begin{equation} \nu_\star{<}2^{-{1}/{3}}\quad \mathrm{or}\quad \nu<2^{-{1}/{3}}k\sigma^{{2}/{3}}; \end{equation}

the forcing must be strong enough to overcome the effect of viscous dissipation. As we just discussed, $\beta$ is absent from this condition, and we stress that (2.32) only applies to the case $\mu _\star =0$, otherwise the instability threshold takes place at a finite $m$, which involves $\beta$ in the stability condition (cf. Srinivasan & Young 2012). In the rest of the paper, we will mainly pay attention to the situation where $\mu _\star \ll \nu _\star \ll 1$ so that (2.32) is relevant, while $\beta _\star$ is of order of unity or larger.

We comment on another interesting relation between our dispersion relation (2.31) and that of Srinivasan & Young (2012): both these dispersion relations have the properties that all unstable modes have real $s_\star$ and instability only exists for $m_\star <1$, despite the difference in the spatial structure of the forcing. Srinivasan & Young (2012) indicated that they could not show these two important properties analytically in their case of ring forcing. However, we can make some progress with our simpler dispersion relation: for real $s_\star$, we need the right-hand side of (2.31) to be positive so that $s_\star >0$ on the left-hand side is possible, and this requires $m_\star <1$. We have not yet been able to explain why $s_\star$ is always real when $\operatorname {\mathrm {Re}} s_\star >0$. In any case, this appears not to be a generic property: Ruiz et al. (2016) have shown that for forcing with other spatial structures, unstable modes with complex growth rates indeed exist for hydrodynamic zonostrophic instability.

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 2012; Parker & Krommes 2013), 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. 2008), when the drag $\mu$ is sufficiently strong (Bouchet et al. 2013), or when there is time scale or space scale separation between the mean flow and the waves (Bouchet et al. 2013; Durston & Gilbert 2016). 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. The MHD zonostrophic instability

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 2018).

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 (2018) 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 (2020), 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.5. Results and discussion

We now discuss the solution for the growth rate $s$ determined by the full dispersion relation (3.36), focusing on the effect of the magnetic field $B_0$ and the magnetic diffusivity $\eta$. We fix the other parameters at $\beta =5$, $\sigma =0.05$, $\nu =10^{-4}$, $\mu =0$, $k=16$ corresponding to $\beta _\star =2.3$, $\nu _\star =0.19$, as used previously for the hydrodynamic case.

Figure 5 shows $s$ calculated from the dispersion relation for various $B_0$ at a relatively high magnetic diffusivity $\eta =10^{-2}$, corresponding to a large $\eta _\star =19$. The solid lines are the numerical solutions of (3.36), and the dashed lines are the results of the asymptotic solution (3.41) derived under the condition of large $\eta _\star$. All the solutions for the growth rate $s$ are real. The agreement between the asymptotic and numerical solutions is good for all of the cases presented. Increasing the magnetic field reduces the growth rate, through damping the flow as discussed above, and finally suppresses the instability.

Figure 5. The dispersion relation giving $s$ for different magnetic field strengths $B_0$ at $\eta =10^{-2}$, corresponding to $\eta _\star =19$. The other parameters are $\beta =5$, $\sigma =0.05$, $\nu =10^{-4}$, $\mu =0$, $k=16$ corresponding to $\beta _\star =2.3$, $\nu _\star =0.19$; $B_0=0.01$ corresponds to $B_{0\star }=1.2$. Solid lines represent numerical solutions, and dashed lines represent the results of the asymptotic solution (3.41).

Figure 6 shows $s$ calculated from the dispersion relation at a moderate magnetic diffusivity $\eta =10^{-4}$, equal to the viscosity $\nu =10^{-4}$. The red lines correspond to real solutions for $s$, and the blue lines give complex solutions $s=s_{r}+\mathrm {i}s_{i}$. Unfortunately, we do not have analytical results for a comparison in this regime and only show numerical solutions. The magnetic field again has a stabilising effect, but the behaviour of the dispersion relation is significantly different from figure 5. In particular, there are two branches of real modes which can be unstable; the lower branch originates from the stable modes at $B_0=0$, corresponding to the dissipation of mean field governed by (3.38). When two branches of real modes merge at a certain wavenumber (e.g. $m\approx 7$ for $B_0=0.01$), a complex mode branches out. As the magnetic field increases, the waveband of real modes shrinks and that of complex modes broadens. Complex modes can also have $s_{r}>0$ and thus be unstable, but the real modes generally have larger growth rates. A magnetic field reduces the real part of the growth rates (except for weak instabilities at large wavenumbers), and suppresses instability when strong enough. The imaginary part $s_{i}$ first increases with the field, and then decreases when the field becomes strong. Our dispersion relation in this case is similar to figure 2 of Constantinou & Parker (2018).

Figure 6. The dispersion relation giving $s$ for different magnetic field strengths $B_0$ at $\eta =10^{-4}$, corresponding to $\eta _\star =0.19$; $B_0=0.01$ corresponds to $B_{0\star }=1.2$. The other parameters are the same as figure 5. The real part $s_{r}$ of the growth rate is shown in (a) and the imaginary part $s_{i}$ in (b). Red lines represent real solutions, and blue lines represent complex solutions.

Figure 7 shows the dispersion relation at a weak diffusivity $\eta =10^{-5}$, corresponding to $\eta _\star =0.019$. In this situation, our approximate solutions give good predictions for different regimes; the dashed lines are the results of the asymptotic dispersion relation (3.43) for weak magnetic field, and the dash–dot lines are the result of (3.45) for a strong field. The dispersion relation behaves very differently from the previous two cases (figures 5 and 6). For weak field $B_0\le 0.0032$, the curves for real modes form a family of closed loops, and complex modes branch out at the left-hand and right-hand edges of each loop. The loops shrink significantly as the field is increased from low levels, a behaviour that can be predicted from (3.43), indicating that the Maxwell stress is responsible for the stabilising effect. The loops of real modes disappear at $B_0\approx 0.0033$, and for fields stronger than this all unstable modes are complex. The growth rate of complex modes then reduces as the field becomes stronger, but less dramatically; instability is suppressed only when $B_0$ reaches the order of $0.1$, corresponding to $B_{0\star }=12$. The behaviour at such strong fields is well captured by the asymptotic dispersion relation (3.45), indicating that the interaction between the mean flow and mean field is the main source for the instability.

Figure 7. The dispersion relation giving $s$ for different magnetic field strengths $B_0$ at $\eta =10^{-5}$, corresponding to $\eta _\star =0.019$; $B_0=10^{-3}$ corresponds to $B_{0\star }=0.12$. The other parameters are the same as figure 5. The solid lines are numerical solutions; the dash and dash–dot lines are the solutions of the asymptotic dispersion relations (3.43) for weak field and (3.45) for strong field, respectively.

In figure 8, we show the stability diagram in the $( B_0, \eta )$ plane for zonostrophic instability in the MHD system, summarising the stable and unstable regions given in figures 5–7. The boundaries between the stable and unstable regions are indicated by the solid lines, with data points represented by the stars. Scaling laws (to be derived shortly) are plotted in dashed lines. If the magnetic field $B_0$ is strong enough, it can suppress the instability for any prescribed $\eta$ and the flow is stable. On the other hand, for any fixed $B_0$, if the diffusivity $\eta$ is high enough, it can counteract the effect of the field and we recover the hydrodynamic instability. There is a region labelled ‘unstable with complex growth rates’, which corresponds to the situation in figure 7(a): for $0.0032< B_0<0.5$, there are no growth rates that are real and positive, but complex growth rates may have positive real parts. For typical flow instability problems, this region would be simply regarded as ‘unstable’, but there are uncertainties around its interpretation for our problem which concern the statistical behaviour of zonostrophic instability. Thus, we regard this region with caution and will discuss its subtle properties later on.

Figure 8. Stability diagram of the mean flow/field zonostrophic instability in the $(B_0,\eta )$-plane. The classification is based on the maximum $s_{r}$ over all wavenumbers $m$. The parameters are $\beta =5$, $\sigma =0.05$, $\nu =10^{-4}$, $\mu =0$, $k=16$ corresponding to $\beta _\star =2.3$, $\nu _\star =0.19$. The region ‘unstable with complex growth rates’ means all real growth rates are negative, but complex growth rates have positive real parts. Dashed lines are the indicated scaling laws.

We now derive scaling laws for the stability boundaries based on our reduced dispersion relations, which will provide new insights and general conclusions. Again, we assume $\mu _\star \ll \nu _\star \ll 1$ while $\beta _\star$ is of order of unity or larger. At high magnetic diffusivity, the instability threshold can be predicted by the asymptotic dispersion relation (3.41). For neutral stability $s=0$, if we balance the viscosity term on the left-hand side with the right-hand side, assuming small $m_\star$ (as we see in figure 5), we have

(3.48)\begin{equation} \nu_\star m_\star^2\sim \frac{m_\star^2}{2\times ({B_{0\star}^2}/{\eta_\star})\times ({B_{0\star}^2}/{\eta_\star})\times 2}\times 2, \end{equation}

which yields a scaling law

(3.49)\begin{equation} B_{0\star}^2\sim \frac{\eta_\star}{\sqrt{2\nu_\star}}\quad \mathrm{or} \quad B_0^2\sim \frac{\eta \sigma}{k\sqrt{2\nu}}. \end{equation}

This scaling law provides the very precise prediction seen in the top right of figure 8.

At low magnetic diffusivity, the boundary between the region of ‘unstable’ and ‘unstable with complex growth rates’ corresponds to when the loops of real modes disappear in figure 7(a), and may be estimated using the asymptotic dispersion relation (3.43). When the loops disappear, $s_\star$ approximately reduces by half, so we may very roughly estimate that the magnetic term in (3.43) reaches half of the hydrodynamic term. Noting $s_\star \gg \nu _\star$ and $m_\star$ is relatively small when the loop disappears, we have the balance

(3.50)\begin{equation} \frac{m_\star^2 B_{0\star}^2}{2\eta_\star \beta_\star^2 s_\star}\sim \frac{1}{2}\times \frac{m_\star^2}{2\nu_\star s_\star}, \end{equation}

and this provides the scaling law

(3.51)\begin{equation} B_{0\star}^2\sim \frac{\eta_\star \beta_\star^2}{2\nu_\star}\quad \mathrm{or}\quad B_{0}^2\sim \frac{\eta \beta^2}{2\nu k^4}. \end{equation}

As seen in figure 8, this scaling law also gives fairly good predictions, and agrees qualitatively with that found by Tobias et al. (2007) and Parker & Constantinou (2019), as we will elaborate subsequently.

Finally, the boundary between the region of ‘unstable for complex growth rates’ and ‘stable’ corresponds to when $s_{r}$ of the complex mode in figure 7(a) changes sign, which can be described by the reduced dispersion relation (3.45). Due to the complicated expressions for $N_{UB\star }$ and $N_{BU\star }$, we are not able to factor out the real part of $s_\star$, and we will only roughly estimate the threshold based on the orders of terms. Numerical solutions suggest that at the stability threshold, the purely imaginary $s_\star$ is of the order of $\nu _\star$ which is small, hence we can deduce

(3.52)\begin{equation} N_{UB\star}, N_{BU\star}=O\left(\frac{1}{\beta_\star B_{0\star}}\right). \end{equation}

Note that the $O(\nu ^{-1}_\star )$ terms are cancelled by the c.c.e.$s$. Then balancing the left-hand and right-hand sides of (3.45) indicates that $B_{0\star }\propto (\beta _\star \nu _\star )^{-1}$. At this point, we can only use data fitting to find the constant factor in this relation: we find the factor $4$ fits the numerical solution well. Hence, we have the scaling law,

(3.53a,b)\begin{equation} B_{0\star}\sim \frac{4}{\beta_\star\nu_\star},\quad B_0\sim \frac{4\sigma^2}{\beta \nu k^2}, \end{equation}

which appears in figure 8 as the vertical dashed line. Although we have fitted the constant to obtain this law, the analysis reveals key underlying physics: we have $B_{0\star }\propto \nu _\star ^{-1}\gg 1$, showing that complex unstable modes survive strong magnetic fields, and that as $\eta _\star$ becomes small it no longer affects the stability boundary.

For moderate magnetic diffusivity $\eta \sim 10^{-4} = \nu$, we do not have an asymptotic dispersion relation, unfortunately. We see that the stability behaviour in figure 8 is rather complicated: the three regions are all present, and the boundaries wobble as $\eta$ increases. The explanation may be that since the magnetic diffusivity is similar in magnitude to the viscosity, the interaction between the flow and the field is very active. All of the mechanisms that we have identified (e.g. increased viscosity, Maxwell stress and mean flow–mean field interaction) are all present. In other words, none of the terms or effects in (3.30) and (3.31) may be neglected.

Tobias et al. (2007) and Constantinou & Parker (2018) also considered the stability diagram in the $(B_0,\eta )$ parameter plane. We show their results in figure 9 and compare them with our figure 8. Tobias et al. (2007) performed numerical simulations and examined the conditions for which the large-scale zonal flows emerge. Their results are shown in figure 9(a), where a plus or diamond sign represents conditions for which a large-scale zonal flow did or did not emerge, respectively. From these data points, they found that the boundary between these two situations obeys a scaling law $B_0^2/\eta =\mathrm {constant}$. Compared with our study, this scaling law corresponds to our result (3.51) at low magnetic diffusivity, which also obeys $B_0^2\propto \eta$.

Figure 9. Stability diagram in the $(B_0, \eta )$ plane taken from (a) Tobias et al. (2007) and (b) Constantinou & Parker (2018) (figures reproduced with permission from the authors). In (a), the plus signs represent simulations in which large-scale zonal flows emerge, and the diamond symbols simulations in which they do not. The solid line is defined by $B_0^2/\eta =\mathrm {constant}$. In (b), the plus signs represent unstable modes with real growth rate, the stars represent complex growth rates, and the circles represent stable modes. The solid line is defined by (3.54).

Constantinou & Parker's (2018) study was based on a zonostrophic instability analysis. One of their results is shown in figure 9(b), where a plus sign represents a real and positive growth rate, a star represents a complex growth rate and a circle represents a real negative growth rate. We see the structure of the three regions is similar to our figure 8 around $\eta \sim 10^{-4}$ (however, we note a difference in definition: their region of complex growth rates counts those with both positive and negative real part, while ours only includes those with positive real part). Their empirical boundary is defined by

(3.54)\begin{equation} \frac{\omega_A^2}{\omega_R^2}\frac{1+Pr_m^2}{Pr_m}=1, \end{equation}

where $\omega _A$ and $\omega _R$ are the frequencies of shear Alfvén and undamped Rossby waves, respectively, and $Pr_m=\nu /\eta$ is the magnetic Prandtl number. Translated to our notation, this scaling law becomes

(3.55)\begin{equation} B_0^2=\frac{\nu\eta \beta^2}{(\eta+\nu)^2k^4}. \end{equation}

For small magnetic field strengths, it agrees with our scaling law (3.51) up to a constant factor. Constantinou & Parker (2018) inferred this scaling law from the form of the spatiotemporal correlation function of the magnetic field. Our derivation further clarifies the physics.

We will return to the discussion of the region of ‘unstable with complex growth rates’ in our figure 8 shortly. Constantinou & Parker (2018) also reported complex growth rates with a positive real part. However, compared with figure 9(a), this region seems to correspond to the conditions where Tobias et al. (2007) found no zonal flow forms, i.e. which should be regarded as zonostrophically stable. The properties of these modes therefore remain curious. We will investigate the behaviour of these complex modes using numerical simulations in the next section, considering in particular the stochastic behaviour of the flow and field.

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. (2007) 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. (2007) 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. (2007) 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. 2016), 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.

4. Conclusions and remarks

In this paper, we have studied zonostrophic instability focusing on its statistical properties and the effect of a magnetic field. We apply a stochastic forcing with its amplitude varying in the form of a white noise, which generates random waves. Weak zonal flows can then grow exponentially to generate strong and stable zonal jets. We study the expectation of the zonal flow and we have derived the dispersion relation for its exponential growth. We have also undertaken numerical simulations of many realisations of stochastic flows to compare with theory, and examined the validity of the widely used mean-flow ergodic assumption, which assumes the zonal mean flow remains the same in different realisations.

In the zonostrophic instability analysis, we have developed a method that does not depend on the spatiotemporal correlation functions employed in many previous studies. Instead, we analyse stochastic waves directly taking advantage of the temporal delta-correlation property of white noise. This allows us to derive simplified dispersion relations in the limits of weak and strong magnetic field and magnetic diffusivity. Our analysis has revealed the key role played by the temporal correlation of the stochastic waves in the zonostrophic instability. Regarding the mean-flow ergodicity, our derivation depends on a weaker assumption compared with previous studies. We assume that the mean flow can still be stochastic, but its stochasticity is uncorrelated to the stochasticity of the waves. We refer to our assumption as a ‘partial ergodic assumption’, in contrast to the ‘full ergodic assumption’ where the mean flow has no stochasticity.

The dispersion relations that we derive provide straightforward insights into the effect of the magnetic field and scaling laws for the instability thresholds. We find that when the magnetic diffusivity is strong, the effect of the magnetic field is equivalent to increasing the damping of the waves. When the magnetic diffusivity is very weak, for weak uniform magnetic fields, it acts through the Maxwell stress to counteract the Reynolds stress, whereas at strong magnetic field, the interaction between the mean flow and mean field is the dominant dynamics. The magnetic field mainly plays a stabilising role and inhibits zonal flows, but in the regime of weak diffusivity and strong field, the interaction between the mean flow and mean field can have a destabilising effect.

In the numerical simulations, we have seen the unstable formation of zonal flow. In order to take into account the stochastic nature of the system, we have run multiple simulations in each parameter regime. Comparing the ensemble average of the mean flow with the growth of expectation predicted by the zonostrophic theory, we find good agreement in general. When comparing the theoretical expectation with individual realisations, we find that there is good agreement in the purely hydrodynamic case, or if the magnetic field is weak or if the magnetic diffusivity is strong. Otherwise, significant differences between individual realisations and the theoretical growth of the expectation typically take place. Specifically, when the growth rate of the zonostrophic instability is complex with positive real part, the individual realisations can have very chaotic behaviours, bearing rather weak relation to the expectation, and whether the zonal flow will emerge for any one realisation is highly unpredictable. Use of the full ergodic assumption, i.e. assuming the mean flow to be the same in different realisations, thus appears to be very questionable in these situations. The partial ergodic assumption, on which our analysis is based, is found to be a better approximation as it allows stochasticity of the mean flow.

The increased randomness of the zonal flow caused by the magnetic field for weak magnetic diffusion deserves further research. For example, what is the mechanism by which the magnetic field makes the zonal flow more disorganised? Is there a way to quantify and analyse this, for example through the covariance of the stochastic process? We may also explore the effect of more realistic forcing, for example, an isotropic ring forcing with wavevectors in all directions. The stochasticity of the zonal mean flow was often neglected in previous studies (e.g. Farrell & Ioannou 2003; Srinivasan & Young 2012). In those cases, the equations of the mean flow are deterministic, controlled by cumulants or correlation functions. Our finding that the magnetic field can make the mean flow much more stochastic needs further study.