On the Rayleigh–Kuo criterion for the tertiary instability of zonal flows

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,¹⁹ 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,²⁸ 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,²⁸ 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.

First, let us introduce the original Hasegawa–Mima equation.³⁸ Consider a collisionless plasma in a uniform magnetic field $B0$ in the z direction, with the equilibrium gradient of the background electron density n₀ 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 E are electrostatic, $E=−∇δφ$⁠, where $δφ(t,x)$ is the corresponding electrostatic potential on the two-dimensional plane $x ≐ (x,y)$⁠. The electron response to E is adiabatic (yet see below), while the ion response can be described by the $E×B0$ drift and the polarization drift. Then, assuming the quasi-neutrality condition, the evolution of δφ is described by the (original) Hasegawa–Mima equation

Here, $ρs≐cs/Ωi$ is the ion sound radius (we use $≐$ to denote definitions), $cs≐ZTe/mi$ is the ion sound speed, Z is the ion charge number, $Ωi≐Z|e|B0/mi$ is the ion gyrofrequency, e is the electron charge, $uE≐ẑ×∇δφ/B0$ is the $E×B0$ velocity, $ẑ$ is the unit vector along the z axis, $V*≐Te/(LnB0|e|)$ is the electron diamagnetic drift velocity, and $Ln≐(−∂ ln n0/∂y)−1$ is the characteristic length scale of n₀. Also, $∇2≐∂2/∂x2+∂2/∂y2$ is the Laplacian.

Let us measure time in units 1/Ω_i and length in units ρ_s. Let us also introduce a normalized potential $φ≐eδφ/Te$ and a normalized “generalized vorticity” $w≐(∇2−1)φ$⁠. Then, Eq. (1) can be written in the following dimensionless form:

where $β≐V*/cs$ is treated as a (positive) constant. Let us also introduce the zonal average as $⟨f⟩≐∫0Lxfdx/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:

where $â$ is an operator such that $â=1$ for DWs and $â=0$ for ZFs.^39,40 Equations (3) and (4) constitute the so-called gHME.⁴¹ Below, we use this model to study the TI.

Consider a stationary ZF with $φ=φ¯(y)$ and $w=φ¯″$⁠. (Hereafter, the prime denotes the derivative with respect to y.) The ZF velocity is $v=U(y)x̂$⁠, where $U(y)≐−φ¯′$ and $x̂$ is the unit vector in the x direction. Consider a DW perturbation $φ̃≐φ−φ¯$ to this ZF. As mentioned in Sec. I, we focus on the case when the initial state is non-turbulent. Then, $φ̃$ is small and can be described by the linearized Eqs. (3) and (4), namely,

Let us search for a solution in the form $φ̃=ϕ(y) exp(ikxx−iωt)$ with constant k_x and ω. Then, $w̃=(∂2/∂y2−kx2−1)φ̃$⁠, and

where $α2≐1+kx2$ and $C≐ω/kx$⁠. This equation is identical⁴² to the barotropic vorticity equation studied by Kuo in Ref. 28. The only difference is that in Kuo's case, α² is not necessarily larger than one. However, we will not make use of our specific expression for α² 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≐d2/dy2−α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_xC has a positive imaginary part, which signifies an instability.

Also 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.]

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 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⁴⁷ 

where ψ(y) is periodic such that $ψ(y+L)=ψ(y)$ for any y, and $q¯$ is a constant. Assuming that ϕ is bounded (as it would be the case, for instance, at periodic boundary conditions), $q¯$ must be real. Without loss of generality, we limit the value of $q¯$ to the first Brillouin zone, i.e., $−q/2≤q¯<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.

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″−β$ must change its sign somewhere in the integration domain; i.e., there must be a location where $U″=β$⁠. (This location can be understood as the point where the “vorticity” Z, defined via $dZ/dy≐U″−β$⁠, has an extremum.⁴⁸) This is the RK necessary condition for a ZF to be unstable.²⁸ For the sinusoidal profile [Eq. (10)], the existence of $U″=β$ 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,³¹ has also been generalized to barotropic and baroclinic flows,⁴⁹ 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″=β$⁠, and add the resulting equation to Eq. (16). Using Eq. (18), this leads to

This means that $(U″−β)(U−U*)$ 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.

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₀ ≤ C ≤ u₀, which requires a special treatment due to possible singularities at U = C. (The case when C < −u₀ or C > u₀ 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²⁸ and have the following asymptotics:

where $Gs≐(Us″−β)/Us′$⁠. (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 $ϕ≈c2, ϕ′≈c2Gs ln (y−ys)$⁠, and $ϕ″≈c2Gs/(y−ys)$⁠.

Note that the above two solutions are invalid if $Us′=0$⁠, but those modes have infinite enstrophy and are physically irrelevant.⁵⁰ Modes with nonzero G_s are physically irrelevant too for the same reason. The only exception is when c₂ = 0, i.e., ϕ = c₁ϕ₁ = 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 α² + β/(C – U) > 0; hence, ϕ oscillates slower than F according to the Sturm comparison theorem.²⁸ 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₂ = 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≐−β/q2$ for a sinusoidal ZF. (The subscript n stands for “neutral”.) Let us substitute C_n into Eq. (7) along with Eq. (11). This gives

This is a linear ordinary differential equation with constant coefficients, so its solutions can be searched in the form $ψ=exp (iλ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¯>0$⁠. (The case $q¯<0$ can be analyzed similarly.) Then, α² ≥ 0 is possible only if m = –1 or m = 0. Therefore, we have two choices of α², and the corresponding solutions are

Since $q¯<q/2$⁠, we have $0≤αn12<αn22≤q2$⁠. This shows that, indeed, C = C_n corresponds to neutral modes. More precisely, two neutral modes correspond to the two choices of α² given by Eqs. (27) and (28).

Now, let us consider $α2=αn2+Δα2$⁠, where Δα² is a small perturbation to a neutral-mode solution (and the subscript n stands for n1 or n2). The corresponding perturbed eigenvalues and eigenmodes are some C = C_n + ΔC and ψ = ψ_n + Δψ. The neutral eigenmode ψ_n satisfies Eq. (24), while ψ is governed by

Let us multiply Eq. (29) by $ψ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 Δα² and ΔC. Hence, we can approximate $ψψn*$ with $|ψn|2$⁠. Therefore

where $N≐∫0L|ψ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 Δα² < 0 and Im(ΔC) < 0 if Δα² > 0; then, unstable eigenmodes exist only when Δα² < 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² but lies out of the range of U (i.e., C_n < –u₀ or C_n > u₀), then $J′(0)$ is real. Accordingly, Eq. (35) gives that ΔC is real, because Δα² 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². Using Eq. (27) or (28), we have $|ψn|2=1$⁠, so $N=L$⁠. Also, the expression for $J′(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′|y1<0$ and $U′|y2>0$ (Fig. 2). Hence, we obtain

where

and

A straightforward calculation gives

where $ϱ≐u0q2/β$ is what we call the Rayleigh–Kuo parameter.¹⁹ Then, using Eq. (35) and L = 2π/q, we obtain

for Δα² < 0. In contrast, if Δα² > 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^*.⁵¹ 

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

Here, $G(Δα2)≐|Δα2|H(−Δα2)$⁠, where H is the Heaviside step function, and

where $αn12$ and $αn22$ are given by Eqs. (27) and (28). As mentioned in Sec. III, we have $0≤αn12<αn22≤q2$⁠; hence, $Δα12>Δα22$⁠. If $Δα12>0$⁠, then $G(Δα12)=|Δα12|H(−Δα12)=0$⁠; hence, $G(Δα22)≥G(Δα12)$⁠; if $Δα12<0$⁠, then $Δα22<Δα12<0$⁠; hence, we still have $G(Δα22)≥G(Δα12)$⁠. This shows that $G(Δα22)≥G(Δα12)$ everywhere, and thus, γ₂ ≥ γ₁ everywhere too. Correspondingly, the largest TI growth rate is given by

at $Δα22<0$⁠, and otherwise γ_TI = 0. Accordingly, the TI develops when

Finally, since $|Δα22|=|α2−αn22|=q2−q¯2−α2$ (assuming $Δα22<0$⁠), the largest γ_TI is realized at $q¯=0$⁠; thus, ϕ(y) = ψ(y) [Eq. (11)]. Then, using Eq. (29) with $q¯=0$ and n → n₂, 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.²⁸ The only difference is that in the latter case, $α2=kx2$⁠; then, α² is allowed to be zero, so Eq. (49) is less stringent, namely, $|ky|≤q$⁠.

Since $0≤αn12<αn22$⁠, the above results indicate that there is no instability at $α2>αn22=q2−q¯2$⁠. Strictly speaking, this conclusion is only valid when α² is close to $αn22$⁠. But one can also prove that no unstable eigenmode exists if $α2≥q2$⁠. 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

Let us multiply Eq. (14) by C – C_n and add the result to Eq. (50). This gives

where $Cr≐Re C$⁠. Since $U=u0 cos qy$ and C_n = –β/q², one finds that $−(U″−β)/(U−Cn)=q2$ is a positive constant. Hence

so we obtain

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

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 $α2≈αn1,22$⁠, which is not the case here. Also, in Fig. 3(b), the instability already vanishes at β = 1.6 < q²u₀ = 2.56, which indicates that ϱ > 1 is only a necessary condition.

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

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

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

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 α²; specifically, in the gHME, one has $α2=1+kx2$⁠, so α² > 1. Hence, the results are as follows. The two necessary conditions for the TI are (i) q²u₀ > β and (ii) q² > 1. The second condition comes from the fact that if q² ≤ 1, then $α2≡1+kx2≥q2$ 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_x, such that $1+kx2$ is smaller but close to q². 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 $γTI≐|kxIm(ΔC)|$ is found to be

Here, H is the Heaviside step function, $ϑ≐1−(q¯2+1+kx2)/q2$⁠, and $ϱ=u0q2/β$⁠. In a dimensional form, the two necessary conditions found above are (i) $q2>ρs−2$ and (ii) $q2u0>ρs−2V*$⁠, and γ_TI is still given by Eq. (56), where

Equations (56) and (57) are among the main results of our paper.

As a corollary, the WKE, which assumes the GO limit that relies on the assumption $q2≪ρs−2$⁠, 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>ρs−2$⁠, 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.⁵³ 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.⁵⁴ 

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′=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″=0$⁠. Also, at arbitrary β, our Eq. (57) shows that when q is large enough, one has $γTI≈|kxu0|$⁠, which too is characteristic of the KHI.

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^k/(2k + 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≐(2k+1)q$⁠. We also define $feff≐qeff/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₀ (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 $ϱeff,min≐qeff2u0,min/β$ 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|≤0.9$⁠. Thus, we conclude that $ϱeff≐qeff2u0/β$ 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≐feffqmin$ 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₀. In particular, consider b = 0, which corresponds to a sinusoidal ZF (hence f_eff = 1). Then, for u₀ = 0.5 < 1, Eq. (55) gives $qmin=2≈1.4$⁠, whereas for u₀ ≥ 1, one obtains q_min = 1. Also, the u₀-dependence becomes weak at u₀ ≥ 3. The results show that q_eff,min decreases slowly at nonzero b, and $qeff,min≳0.5$ within the range $|b|≤0.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.

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;²⁸ 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.

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.