This paper reports the stability conditions for intense zonal flows (ZFs) and the growth rate *γ*_{TI} of the corresponding “tertiary” instability (TI) within the generalized Hasegawa–Mima plasma model. The analytical calculation extends and revises Kuo's analysis of the mathematically similar barotropic vorticity equation for incompressible neutral fluids on a rotating sphere [H.-L. Kuo, J. Meteor. **6**, 105 (1949)]; then, the results are applied to the plasma case. An error in Kuo's original result is pointed out. An explicit analytical formula for *γ*_{TI} is derived and compared with numerical calculations. It is shown that, within the generalized Hasegawa–Mima model, a sinusoidal ZF is TI-unstable if and only if it satisfies the Rayleigh–Kuo criterion (known from geophysics) and that the ZF wave number exceeds the inverse ion sound radius. For non-sinusoidal ZFs, the results are qualitatively similar. As a corollary, there is no TI in the geometrical-optics limit, i.e., when the perturbation wavelength is small compared to the ZF scale. This also means that the traditional wave kinetic equation, which is derived under the geometrical-optics assumption, cannot adequately describe the ZF stability.

## I. INTRODUCTION

Sheared plasma flows driven by turbulence, which are also known as zonal flows (ZFs), significantly affect transport in magnetically confined plasmas.^{1–6} Hence, they have been actively studied in the literature. The linear stage of the zonostrophic instability (ZI) that produces ZFs out of homogeneous turbulence is now largely understood.^{7–14} In contrast, the factors that limit the amplitudes of the nonlinear ZFs have not been identified with certainty. The well-known predator-prey model^{15–18} predicts that the saturation amplitudes depend on the ZF collisional damping rates, which are introduced in an *ad hoc* manner. However, the predator-prey oscillations are also possible in the absence of dissipation,^{19} so the ZF damping rate may not be the only important factor. One may wonder then whether simple parameters can be identified that constrain the ZF amplitude more robustly, i.e., without invoking *ad hoc* energy losses.

Here, we study a specific aspect of this problem, namely, the instability of a prescribed ZF, which is known as the tertiary instability (TI).^{20–26} The plasma is modeled within the generalized Hasegawa–Mima equation (gHME), so our definition of the TI is different from that in Refs. 20 and 21, where this instability was attributed to the ion-temperature gradient (absent in the gHME). However, our definition of the TI is similar to those in the majority of relevant papers.^{22–25}

As pointed out in Ref. 27, the gHME is analogous to the barotropic vorticity equation that describes neutral flows in the atmospheres of rotating planets,^{28,29} where the Coriolis parameter *β* plays the same role as the density gradient in the gHME. Therefore, the TI can be understood as analogous to the Kelvin–Helmholtz instability (KHI)^{30,31} of a neutral flow modified by nonzero *β*. (The analogy between the KHI and the TI is also mentioned in Refs. 22 and 23.) The latter modification was first studied by Kuo,^{28} who generalized the famous Rayleigh's inflection theorem for the KHI to the case of nonzero *β*. The resulting criterion, known as the Rayleigh–Kuo (RK) necessary condition for instability, is widely cited in the geophysics literature but has not been popular in the plasma physics literature with only a few exceptions.^{23,32} In Ref. 23, it was mentioned that the criterion provides a good estimate for the TI threshold, but no rigorous analysis was presented that would address the necessary and sufficient conditions for the TI. This and the fact that Kuo's results are not entirely accurate (see Sec. III) warrant a careful examination of the subject.

In this paper, we restate the RK criterion within the gHME model. Specifically, we identify omissions in Kuo's original paper,^{28} propose an explicit formula for the TI growth rate, and compare it with numerical calculations. It is shown that, within the generalized Hasegawa–Mima model, a sinusoidal ZF is TI-unstable if and only if it satisfies the Rayleigh–Kuo criterion (known from geophysics) and that the ZF wave number exceeds the inverse ion sound radius. For non-sinusoidal ZFs, the results are qualitatively similar. As a corollary, there is no TI in the geometrical-optics (GO) limit, i.e., when the perturbation wavelength is small compared to the ZF scale. This also means that the traditional wave kinetic equation (WKE)^{1,7,8,15,16,22,33,34} cannot adequately describe the ZF stability. In particular, the WKE-based analysis of the TI in Ref. 22 actually addresses a different instability, namely, a branch of the ZI, as will be explained in Sec. V.

Note that, in order to produce a rigorous analytical theory, we simplify the problem by limiting our consideration to the strong-ZF case, namely, the case when the ambient turbulence is negligible and the ZF can be considered laminar. The more general case, when the ambient turbulence is allowed to affect the ZF stability, was studied using the stochastic structural stability theory in Refs. 12–14 and using the equivalent second-order cumulant expansion theory in Refs. 9 and 11. In particular, Refs. 11–14 reported that the ZF-turbulence system undergoes a structural instability even when the laminar ZF with the same amplitude would be RK-stable. Also, Ref. 14 shows that the least-damped eigenmode changes its structure in the presence of the ambient turbulence and hence can become unstable. Correspondingly, the stability criterion that we report here can be considered as determining the *upper bound* of the ZF amplitude in a stable equilibrium. (Notably, the roles of stable eigenmodes are also discussed in Refs. 35–37.)

Our paper is organized as follows. The basic equations are introduced in Sec. II. In Sec. III, we obtain the two necessary conditions and analytically calculate the TI growth rate in the context that includes both geophysical and plasma settings equally. The comparison between the analytical formula for the growth rate and the numerically found eigenvalues is presented in Sec. IV. The ramifications of our theory that are specific to the plasma case (as opposed to Kuo's geophysical problem) are discussed in Sec. V. The generalization to non-sinusoidal ZF is discussed in Sec. VI. Our main conclusions are summarized in Sec. VII.

## II. BASIC EQUATIONS

### A. The generalized Hasegawa–Mima equation

First, let us introduce the original Hasegawa–Mima equation.^{38} Consider a collisionless plasma in a uniform magnetic field $B0$ in the *z* direction, with the equilibrium gradient of the background electron density *n*_{0} in the *y* direction (Fig. 1). Ions are assumed cold, while electrons are assumed to have a finite temperature *T _{e}*. Suppose that perturbations to the electric field

**are electrostatic, $E=\u2212\u2207\delta \phi $, where $\delta \phi (t,x)$ is the corresponding electrostatic potential on the two-dimensional plane $x\u2009\u2250\u2009(x,y)$. The electron response to**

*E***is adiabatic (yet see below), while the ion response can be described by the $E\xd7B0$ drift and the polarization drift. Then, assuming the quasi-neutrality condition, the evolution of**

*E**δφ*is described by the (original) Hasegawa–Mima equation

Here, $\rho s\u2250cs/\Omega i$ is the ion sound radius (we use $\u2250$ to denote definitions), $cs\u2250ZTe/mi$ is the ion sound speed, *Z* is the ion charge number, $\Omega i\u2250Z|e|B0/mi$ is the ion gyrofrequency, *e* is the electron charge, $uE\u2250z\u0302\xd7\u2207\delta \phi /B0$ is the $E\xd7B0$ velocity, $z\u0302$ is the unit vector along the *z* axis, $V*\u2250Te/(LnB0|e|)$ is the electron diamagnetic drift velocity, and $Ln\u2250(\u2212\u2202\u2009ln\u2009n0/\u2202y)\u22121$ is the characteristic length scale of *n*_{0}. Also, $\u22072\u2250\u22022/\u2202x2+\u22022/\u2202y2$ is the Laplacian.

Let us measure time in units 1/Ω_{i} and length in units *ρ _{s}*. Let us also introduce a normalized potential $\phi \u2250e\delta \phi /Te$ and a normalized “generalized vorticity” $w\u2250(\u22072\u22121)\phi $. Then, Eq. (1) can be written in the following dimensionless form:

where $\beta \u2250V*/cs$ is treated as a (positive) constant. Let us also introduce the zonal average as $\u27e8f\u27e9\u2250\u222b0Lxfdx/Lx$, where *L _{x}* is the system length in the

*x*direction. Then, perturbations governed by Eq. (1) include ZFs and drift waves (DWs). The former are identified as zonal-averaged perturbations, and the latter are identified as fluctuations with zero zonal average. Strictly speaking, electrons respond differently to ZFs and drift waves (DWs). Specifically, the above model can be made more realistic if one rewrites the governing equations as follows:

### B. Reduction to an eigenvalue problem

Consider a stationary ZF with $\phi =\phi \xaf(y)$ and $w=\phi \xaf\u2033$. (Hereafter, the prime denotes the derivative with respect to *y*.) The ZF velocity is $v=U(y)x\u0302$, where $U(y)\u2250\u2212\phi \xaf\u2032$ and $x\u0302$ is the unit vector in the *x* direction. Consider a DW perturbation $\phi \u0303\u2250\phi \u2212\phi \xaf$ to this ZF. As mentioned in Sec. I, we focus on the case when the initial state is non-turbulent. Then, $\phi \u0303$ is small and can be described by the linearized Eqs. (3) and (4), namely,

Let us search for a solution in the form $\phi \u0303=\varphi (y)\u2009exp(ikxx\u2212i\omega t)$ with constant *k _{x}* and

*ω*. Then, $w\u0303=(\u22022/\u2202y2\u2212kx2\u22121)\phi \u0303$, and

where $\alpha 2\u22501+kx2$ and $C\u2250\omega /kx$. This equation is identical^{42} to the barotropic vorticity equation studied by Kuo in Ref. 28. The only difference is that in Kuo's case, *α*^{2} is not necessarily larger than one. However, we will not make use of our specific expression for *α*^{2} until Sec. V, where we will consider the application of our results to the gHME explicitly. In this sense, our main results are also applicable in the geophysical context.

Equation (7) can be represented as

$A\u2250d2/dy2\u2212\alpha 2$, so it is understood as an eigenvalue problem (with certain boundary conditions). Without loss of generality, we assume *k _{x}* > 0. Hence, if Im

*C*>

*0, then*

*ω*≡

*k*has a positive imaginary part, which signifies an instability.

_{x}CAlso note that the complex conjugate of Eq. (7) is

This shows that if *ϕ* is a solution corresponding to a certain phase velocity *C*, then *ϕ** is also a solution corresponding to *C**. Hence, an instability is possible whenever Im *C* is nonzero, because a mode with Im *C* of any sign is always accompanied by a mode with Im *C* of the opposite sign. [Notably, this is not the case at nonzero viscosity, because then Eq. (7) acquires complex coefficients.]

### C. Floquet analysis

For simplicity, let us assume (until Sec. VI) an unbounded sinusoidal ZF profile

(This periodic profile, although different from those typically used in studies of neutral fluids^{28–30,43,44} with fixed boundaries, is often used in numerical or theoretical modeling in plasma physics.^{23,45,46}) For clarity, we assume *u*_{0} > 0 and *q *>* *0. Then, Eq. (7) is an ordinary differential equation with periodic coefficients whose period is *L *=* *2*π*/*q*. Then, any solution of Eq. (7) is decomposable into Floquet modes of the form^{47}

where *ψ*(*y*) is periodic such that $\psi (y+L)=\psi (y)$ for any *y*, and $q\xaf$ is a constant. Assuming that *ϕ* is bounded (as it would be the case, for instance, at periodic boundary conditions), $q\xaf$ must be real. Without loss of generality, we limit the value of $q\xaf$ to the first Brillouin zone, i.e., $\u2212q/2\u2264q\xaf<q/2$.

In what follows, we explore conditions under which these restrictions lead to complex *C*, i.e., an instability. Any such instability is by definition considered a TI.

## III. INSTABILITY ONSET

### A. The Rayleigh–Kuo criterion

First, let us repeat the RK argument for completeness. Let us multiply Eq. (7) by *ϕ** and consider the imaginary part of the resulting equation

By integrating this over *y* from 0 to *L*, we obtain

where we used Eq. (11) and the periodicity of *ψ*. If Im *C *≠* *0, then we obtain

Hence, $U\u2033\u2212\beta $ must change its sign somewhere in the integration domain; i.e., there must be a location where $U\u2033=\beta $. (This location can be understood as the point where the “vorticity” *Z*, defined via $dZ/dy\u2250U\u2033\u2212\beta $, has an extremum.^{48}) This is the RK necessary condition for a ZF to be unstable.^{28} For the sinusoidal profile [Eq. (10)], the existence of $U\u2033=\beta $ means that

The RK criterion is a generalization of Rayleigh's famous inflection-point theorem to the case of nonzero *β*. Another famous necessary condition, the Fjørtoft's theorem,^{31} has also been generalized to barotropic and baroclinic flows,^{49} and we briefly mention it here for the sake of completeness. Let us multiply Eq. (7) by *ϕ** and integrate the result over *y* from 0 to *L*

From Eq. (11), we have

Therefore, the left-hand side of Eq. (16) is

Let us multiply Eq. (14) with *C*^{*} – *U*_{*}, where *U*_{*} is the value of *U* where $U\u2033=\beta $, and add the resulting equation to Eq. (16). Using Eq. (18), this leads to

This means that $(U\u2033\u2212\beta )(U\u2212U*)$ must be negative somewhere, which is the generalization of Fjørtoft's theorem in the case of nonzero *β*. For the sinusoidal profile, this criterion is always satisfied if the RK criterion is satisfied.

### B. Neutral eigenmodes

Equation (7) has a special class of solutions called “neutral modes” that correspond to real *C*. (Such perturbations neither grow nor decay in time.) Below, we consider the case –*u*_{0} ≤ *C* ≤ *u*_{0}, which requires a special treatment due to possible singularities at *U *=* C*. (The case when *C* < −*u*_{0} or *C* > *u*_{0} will be considered in Sec. III C.) Near a singularity *y* = *y _{s}*, two linearly independent solutions

*ϕ*

_{1,2}of Eq. (7) can be obtained using the Frobenius method

^{28}and have the following asymptotics:

where $Gs\u2250(Us\u2033\u2212\beta )/Us\u2032$. (The subscript *s* means that the corresponding functions are evaluated at *y* = *y _{s}*.) The general solution can be written as

Near the singularity, one has $\varphi \u2248c2,\u2009\varphi \u2032\u2248c2Gs\u2009ln\u2009(y\u2212ys)$, and $\varphi \u2033\u2248c2Gs/(y\u2212ys)$.

Note that the above two solutions are invalid if $Us\u2032=0$, but those modes have infinite enstrophy and are physically irrelevant.^{50} Modes with nonzero *G _{s}* are physically irrelevant too for the same reason. The only exception is when

*c*

_{2}= 0, i.e.,

*ϕ*=

*c*

_{1}

*ϕ*

_{1}= 0 at the singularity. However, this is impossible, which is seen as follows. Let us compare Eq. (7) with the following equation:

Between two neighboring singularities, we have *C *>* U* and *α*^{2} + *β*/(*C* – *U*) > 0; hence, *ϕ* oscillates slower than *F* according to the Sturm comparison theorem.^{28} However, *F *=* U* – *C* is a solution that is zero at both singularities, which means that *ϕ* cannot be zero at both singularities. This rules out the possibility that *c*_{2} = 0. (Note that this conclusion relies on the specific profile we chose here. It may fail if there is only one singularity, e.g., for a monotonic *U*.)

From the above discussion, we conclude that only those neutral modes that have *G _{s}* = 0 need be considered. This corresponds to $C=Cn\u2250\u2212\beta /q2$ for a sinusoidal ZF. (The subscript

*n*stands for “neutral”.) Let us substitute

*C*into Eq. (7) along with Eq. (11). This gives

_{n}This is a linear ordinary differential equation with constant coefficients, so its solutions can be searched in the form $\psi =exp\u2009(i\lambda y)$. Then, we obtain

Since *ψ* is periodic in *y*, we require that *λ* = *mq*, where *m* can be any integer. Then, Eq. (25) becomes

For clarity, suppose $q\xaf>0$. (The case $q\xaf<0$ can be analyzed similarly.) Then, *α*^{2} ≥ 0 is possible only if *m* = –1 or *m *=* *0. Therefore, we have two choices of *α*^{2}, and the corresponding solutions are

### C. Unstable eigenmodes

Now, let us consider $\alpha 2=\alpha n2+\Delta \alpha 2$, where Δ*α*^{2} is a small perturbation to a neutral-mode solution (and the subscript *n* stands for *n*1 or *n*2). The corresponding perturbed eigenvalues and eigenmodes are some *C* = *C _{n}* + Δ

*C*and

*ψ*=

*ψ*+ Δ

_{n}*ψ*. The neutral eigenmode

*ψ*satisfies Eq. (24), while

_{n}*ψ*is governed by

Let us multiply Eq. (29) by $\psi n*$ and the complex conjugate of Eq. (24) by *ψ*. By subtracting one result from the other results, we obtain

Let us integrate this equation over *y* from 0 to *L*. Using the periodicity of *ψ _{n}* and

*ψ*, we obtain

In the above integral, the term in the bracket is at least of the first order of Δ*α*^{2} and Δ*C*. Hence, we can approximate $\psi \psi n*$ with $|\psi n|2$. Therefore

where $N\u2250\u222b0L|\psi n|2dy$ and

The integrand is analytic in the whole integration domain if Im(Δ*C*) is nonzero. In Ref. 28, Kuo Taylor-expanded this integrand in terms of Δ*C* for all *y* and concluded that Im(Δ*C*) > 0 if Δ*α*^{2} < 0 and Im(Δ*C*) < 0 if Δ*α*^{2} > 0; then, unstable eigenmodes exist only when Δ*α*^{2} < 0. However, Δ*C* cannot be considered small near *U* = *C _{n}*; hence, the Taylor expansion in Δ

*C*is invalid there. The correct way to calculate

*J*is as follows. Since

*J*(Δ

*C*= 0) = 0, to the first order we have

In combination with Eq. (32), this gives that, to the first order in Δ*C*,

Here, the prime denotes the derivative with respect to Δ*C*, unlike in the rest of the text, where the prime denotes the derivative with respect to *y*.

If *C _{n}* is not equal to –

*β*/

*q*

^{2}but lies out of the range of

*U*(i.e.,

*C*< –

_{n}*u*

_{0}or

*C*>

_{n}*u*

_{0}), then $J\u2032(0)$ is real. Accordingly, Eq. (35) gives that Δ

*C*is real, because Δ

*α*

^{2}and $N$ are real by definition. In other words, perturbations to such neutral modes are stable and thus do not need to be considered for our purposes.

Now, let us consider *C _{n}* = –

*β*/

*q*

^{2}. Using Eq. (27) or (28), we have $|\psi n|2=1$, so $N=L$. Also, the expression for $J\u2032(0)$ can be found as the limit of

at Δ*C* → 0. Since the integrand has a first-order pole at *U* = *C _{n}* when Δ

*C*= 0, it is convenient to use the Plemelj formula

where P.V. stands for the principal value. Then, we have

where we used $U\u2032|y1<0$ and $U\u2032|y2>0$ (Fig. 2). Hence, we obtain

where

and

A straightforward calculation gives

where $\u03f1\u2250u0q2/\beta $ is what we call the Rayleigh–Kuo parameter.^{19} Then, using Eq. (35) and *L *=* *2*π*/*q*, we obtain

for Δ*α*^{2} < 0. In contrast, if Δ*α*^{2} > 0, no self-consistent solution exists. This result is at variance with Kuo's result but in agreement with our earlier observation (Sec. II) that, if *C* is an eigenvalue, then so is *C*^{*}.^{51}

Since there are two neutral modes, the above calculation gives two branches of unstable modes. Their growth rates are given by $\gamma =|kxIm(\Delta C)|$, namely,

Here, $G(\Delta \alpha 2)\u2250|\Delta \alpha 2|H(\u2212\Delta \alpha 2)$, where *H* is the Heaviside step function, and

where $\alpha n12$ and $\alpha n22$ are given by Eqs. (27) and (28). As mentioned in Sec. III, we have $0\u2264\alpha n12<\alpha n22\u2264q2$; hence, $\Delta \alpha 12>\Delta \alpha 22$. If $\Delta \alpha 12>0$, then $G(\Delta \alpha 12)=|\Delta \alpha 12|H(\u2212\Delta \alpha 12)=0$; hence, $G(\Delta \alpha 22)\u2265G(\Delta \alpha 12)$; if $\Delta \alpha 12<0$, then $\Delta \alpha 22<\Delta \alpha 12<0$; hence, we still have $G(\Delta \alpha 22)\u2265G(\Delta \alpha 12)$. This shows that $G(\Delta \alpha 22)\u2265G(\Delta \alpha 12)$ everywhere, and thus, *γ*_{2} ≥ *γ*_{1} everywhere too. Correspondingly, the largest TI growth rate is given by

at $\Delta \alpha 22<0$, and otherwise *γ*_{TI} = 0. Accordingly, the TI develops when

Finally, since $|\Delta \alpha 22|=|\alpha 2\u2212\alpha n22|=q2\u2212q\xaf2\u2212\alpha 2$ (assuming $\Delta \alpha 22<0$), the largest *γ*_{TI} is realized at $q\xaf=0$; thus, *ϕ*(*y*) = *ψ*(*y*) [Eq. (11)]. Then, using Eq. (29) with $q\xaf=0$ and *n* → *n*_{2}, we obtain the characteristic wavenumber of *ψ*(*y*) as

Hence, the characteristic spatial scale of the TI mode is actually larger than that of the ZF, so the TI cannot be described within the GO approximation in principle. Also note that the same conclusion holds also for the barotropic vorticity equation.^{28} The only difference is that in the latter case, $\alpha 2=kx2$; then, *α*^{2} is allowed to be zero, so Eq. (49) is less stringent, namely, $|ky|\u2264q$.

### D. There is no instability at $\alpha 2\u2265q2$

Since $0\u2264\alpha n12<\alpha n22$, the above results indicate that there is no instability at $\alpha 2>\alpha n22=q2\u2212q\xaf2$. Strictly speaking, this conclusion is only valid when *α*^{2} is close to $\alpha n22$. But one can also prove that no unstable eigenmode exists if $\alpha 2\u2265q2$. Our proof follows the argument from Ref. 30. The difference is that we extend that argument by allowing for nonzero *β*, which is done as follows.

Let us multiply Eq. (7) by *ϕ** and integrate the result over *y* from 0 to *L*. After integrating by parts, we obtain

where $Cr\u2250Re\u2009C$. Since $U=u0\u2009cos\u2009qy$ and *C _{n}* = –

*β*/

*q*

^{2}, one finds that $\u2212(U\u2033\u2212\beta )/(U\u2212Cn)=q2$ is a positive constant. Hence

so we obtain

This shows that no unstable eigenmode is possible if *α*^{2} ≥ *q*^{2}. In application to the gHME, where $\alpha 2\u22611+kx2$, this indicates that not only do we have $|ky|<q$ but the ZF itself must also have *q*^{2} ≥ *α*^{2} > 1 for the TI onset. In other words, a ZF is *always* stable in the GO limit, because in this limit, we have *q*^{2} ≪ 1, and thus, *α*^{2} > *q*^{2} automatically (see also Sec. V).

## IV. COMPARISON WITH NUMERICAL CALCULATIONS

Here, we compare our approximate analytical result [Eq. (44)] with numerical solutions of Eq. (7). Specifically, we search for *ϕ* in the form of a Floquet mode

where the series has been truncated at a large enough *n *=* N*. Then, the eigenvalues of Eq. (8) are found numerically.

Figures 3(a) and 3(b) show the dependence of Im(Δ*C*) on *β*. It is seen that the perturbation theory [Eq. (44)] approximates the numerical result with reasonable accuracy in Fig. 3(a). The discrepancy in Fig. 3(b) is due to the fact that the perturbation theory assumes $\alpha 2\u2248\alpha n1,22$, which is not the case here. Also, in Fig. 3(b), the instability already vanishes at *β* = 1.6 < *q*^{2}*u*_{0} = 2.56, which indicates that ϱ > 1 is only a necessary condition.

Next, we consider the change of Im(Δ*C*) with *α*^{2} in Figs. 3(c) and 3(d). Since $q\xaf=0$, the two neutral eigenmodes are at $\alpha n12=0$ and $\alpha n22=q2=2.56$. We also see that there are no unstable eigenvalues at Δ*α*^{2} > 0, which agrees with the perturbation theory. Unstable eigenvalues exist when $\alpha 2\u2272\alpha n22=q2$, whose numerical values agree with perturbation theory when $\alpha 2\u2248\alpha n22$. However, the range of *α*^{2} where Im(Δ*C*) > 0 depends on *β*; namely, the range is smaller when *β* is larger.

Finally, let us consider the dependence on $q\xaf$. In Figs. 3(e) and 3(f), we plot the dependence of Im(Δ*C*) on *α*^{2} for nonzero $q\xaf$. Then, $\alpha n12=q2\u2212q\xaf2<q2$ and $\alpha n22=q2\u2212(q\u2212q\xaf)2>0$. In Fig. 3(e), $q\xaf=0.1$ is small, and the existence of Im(Δ*C*) > 0 agrees with perturbation theory. In Fig. 3(f), $q\xaf=0.6$ is large, so $\alpha n12$ and $\alpha n22$ are close to each other. Thus, the branch of unstable eigenvalues starting from $\alpha n22$ intersects with the branch starting from $\alpha n12$.

In summary, our first-order perturbation theory agrees with numerical results when *α*^{2} is close to $\alpha n1,22$. It also captures the qualitative dependence on *β*. However, the dependence on *α*^{2} away from $\alpha n1,22$ is not well described. An alternative approximation for *γ*_{TI}, which is not asymptotically accurate but applicable at all *α*^{2}, was proposed in Ref. 19 based on a different argument.

## V. TERTIARY INSTABILITY IN PLASMAS

Here, we consider the ramifications of our theory that are specific to the plasma problem governed by the gHME as opposed to Kuo's geophysics problem. The difference is in the definition of *α*^{2}; specifically, in the gHME, one has $\alpha 2=1+kx2$, so *α*^{2} > 1. Hence, the results are as follows. The two necessary conditions for the TI are (i) *q*^{2}*u*_{0} > *β* and (ii) *q*^{2} > 1. The second condition comes from the fact that if *q*^{2} ≤ 1, then $\alpha 2\u22611+kx2\u2265q2$ for any real *k _{x}*, and hence, no instability is possible. We also emphasize that the two necessary conditions combined together are also

*sufficient*for the TI. The reason for the sufficiency is that when the two necessary conditions are satisfied, one can find a nonzero

*k*, such that $1+kx2$ is smaller but close to

_{x}*q*

^{2}. Then, according to the analysis in Sec. III C, there exists an unstable eigenmode. Hence, we conclude that for a sinusoidal ZF, the necessary and sufficient condition for the onset of TI is

as also illustrated in Fig. 4(a).

As discussed in Sec. III C, the TI growth rate $\gamma TI\u2250|kxIm(\Delta C)|$ is found to be

Here, *H* is the Heaviside step function, $\u03d1\u22501\u2212(q\xaf2+1+kx2)/q2$, and $\u03f1=u0q2/\beta $. In a dimensional form, the two necessary conditions found above are (i) $q2>\rho s\u22122$ and (ii) $q2u0>\rho s\u22122V*$, and *γ*_{TI} is still given by Eq. (56), where

As a corollary, the WKE, which assumes the GO limit that relies on the assumption $q2\u226a\rho s\u22122$, is not adequate to describe the TI. This conclusion applies to both the traditional WKE^{1,7,8,15,16,22,33,34} and the “improved” WKE proposed recently in Refs. 32 and 46. (The improved WKE too relies on the assumption that *q* is small compared to the characteristic DW wavelength.) An adequate theory of the TI must not assume the GO approximation. This is also discussed in Ref. 19, where an alternative (but similar) approximation for *γ*_{TI} is derived from different arguments.

Finally, let us compare our findings with those in other studies of the TI. First, our findings support the conjecture in Ref. 23 regarding the relevance of the RK criterion, and our approximate formula for *γ*_{TI} [Eq. (56)] is in general agreement with the numerical results in that paper. The second part of our instability criterion, $q2>\rho s\u22122$, is not mentioned in Ref. 23 explicitly but is satisfied for the simulation parameters presented there. Our approach is more rigorous than the truncated-Floquet approach used in Refs. 19, 22, 24, and 26 in terms of predicting the onset of the instability. (A comparison between the growth rates obtained from the two approaches can be found in Ref. 19.) In particular, the parameter *β* does not enter the final dispersion relation in Ref. 22, so the RK criterion is missed.^{53} The discussion of the “generalized KHI” in Ref. 22, as mentioned before, also seems irrelevant to the TI. This is because the generalized KHI is excited by a homogeneous background, which makes it just another branch of the ZI. In any case, as is shown above, the WKE cannot capture the TI in principle, because the underlying GO approximation will always lead to the ZF amplification rather than deterioration.^{54}

Our results are also different from those in Refs. 20 and 21 in the following sense. In those papers, the TI is considered as a mode driven by the ion temperature gradient, which is absent in our model. Also, the mode structure in Refs. 20 and 21 is found to be localized where $U\u2032=0$. In contrast, the TI considered in our paper (as well as in Refs. 22 and 23) is similar to the KHI. The reason is that when *β* = 0, our basic equation (7) can also be used to describe the KHI, in which case the mode amplitude peaks at $U\u2033=0$. Also, at arbitrary *β*, our Eq. (57) shows that when *q* is large enough, one has $\gamma TI\u2248|kxu0|$, which too is characteristic of the KHI.

## VI. NON-SINUSOIDAL ZONAL FLOW

Finally, let us generalize our conclusions to nonsinusoidal flows by considering a class of ZFs with a single control parameter *b*, namely,

where *U _{k}* =

*b*/(2

^{k}*k*+ 1) and –1 <

*b*<

*1 [Fig. 5(a)]. When*

*b*is zero, we recover Eq. (10). When

*b*is nonzero, the ZF contains multiple harmonics, and in order to describe its characteristic length scale, we define the effective (weighted) wavenumber as

where $qk\u2250(2k+1)q$. We also define $feff\u2250qeff/q$ [Fig. 5(b)], which increases as $|b|$ increases.

Knowing the non-sinusoidal ZF profile [Eq. (58)], we can find numerical eigenvalues using the same procedure as in Sec. IV. First, we consider a fixed *q* and numerically search for minimum *u*_{0} (denoted by *u*_{0,min}) below which unstable eigenvalues disappear for all *k _{x}*. For example, Fig. 6(a) shows the results for

*b*= –0.5. It is seen that $\u03f1eff,min\u2250qeff2u0,min/\beta $ remains approximately constant and close to one, namely, ϱ

_{eff,min}≈ 1.2. Similar results are obtained for other choices of

*b*within the range $|b|\u22640.9$. Thus, we conclude that $\u03f1eff\u2250qeff2u0/\beta $ can be considered as an effective RK parameter for non-sinusoidal ZFs; i.e., the first TI criterion holds approximately in the form

Next, we consider a fixed *b* and numerically search for the minimum *q* (denoted by *q*_{min}) below which unstable eigenvalues disappear for all *k _{x}*. We find that for nonzero

*b*,

*q*

_{min}is smaller than one, and

*q*

_{min}approaches zero as $|b|$ approaches one. However, $qeff,min\u2250feffqmin$ remains of order one. An example with

*β*= 1 is given in Fig. 6(b), where we show different results depending on the choice of

*u*

_{0}. In particular, consider

*b*=

*0, which corresponds to a sinusoidal ZF (hence*

*f*

_{eff}= 1). Then, for

*u*

_{0}= 0.5 < 1, Eq. (55) gives $qmin=2\u22481.4$, whereas for

*u*

_{0}≥ 1, one obtains

*q*

_{min}= 1. Also, the

*u*

_{0}-dependence becomes weak at

*u*

_{0}≥ 3. The results show that

*q*

_{eff,min}decreases slowly at nonzero

*b*, and $qeff,min\u22730.5$ within the range $|b|\u22640.9$. Therefore, we conclude that the second necessary condition for the TI still holds approximately for non-sinusoidal ZFs, namely, in the form

Finally, the numerically found stability diagram for nonzero *b* is shown in Fig. 4(b). It is seen that the diagram is very similar to the sinusoidal case. Therefore, just like for sinusoidal ZFs, the two necessary conditions combined together [i.e., Eqs. (60) and (61)] are also sufficient for the TI.

## VII. CONCLUSIONS

In summary, we explored the tertiary instability (TI) of zonal flows within the gHME model. Our analytical calculation extends and revises Kuo's analysis of the mathematically similar barotropic vorticity equation for incompressible neutral fluids on a rotating sphere;^{28} then, the results are applied to the plasma case. An error in Kuo's original results is pointed out. An explicit analytical formula for the TI growth rate *γ*_{TI} is derived [Eqs. (44) and (57)] and compared with numerical calculations. It is shown that, within the generalized Hasegawa–Mima model, a sinusoidal ZF is TI-unstable if and only if it satisfies the Rayleigh–Kuo criterion (known from geophysics) and that the ZF wave number exceeds the inverse ion sound radius. For non-sinusoidal ZFs, the results are qualitatively similar. As a corollary, there is no TI in the GO limit, i.e., when the perturbation wavelength is small compared to the ZF scale. This also means that the traditional wave kinetic equation, which is derived under the GO assumption, cannot adequately describe the ZF stability.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy (DOE), Office of Science, Office of Basic Energy Sciences, and also by the U.S. DOE through Contract No. DE-AC02-09CH11466.

## References

The gHME conserves the total enstrophy. Hence, provided that the initial conditions are physical (smooth enough), infinite-enstrophy modes cannot be excited.

As noted in Sec. II B, this is not the case in the presence of viscosity. Nevertheless, this does not reinstate Kuo's argument, because the latter was developed within a model of a non-viscous fluid.

The second part of our instability criterion, $q2>\rho s\u22122$, is less relevant within the model adopted in Ref. 22. This is because Ref. 22 considers the TI as a flute-like mode with small $k\u2225$, where the non-adiabatic electron response is assumed. Therefore, our $\alpha 2=1+kx2$ must be replaced with $\alpha 2=kx2$. Then, $\alpha 2<q2$ can also be satisfied even within the GO assumption.