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.

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 model15–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.

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 n0 in the y direction (Fig. 1). Ions are assumed cold, while electrons are assumed to have a finite temperature Te. 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

t[(ρs221)δφ]+uE·[(ρs221)δφ]+V*δφx=0.
(1)

Here, ρscs/Ωi is the ion sound radius (we use to denote definitions), csZTe/mi is the ion sound speed, Z is the ion charge number, ΩiZ|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(lnn0/y)1 is the characteristic length scale of n0. Also, 22/x2+2/y2 is the Laplacian.

FIG. 1.

The assumed coordinate system. Here, B0 is the magnetic field, n0 is the background electron density, and v is the ZF velocity.

FIG. 1.

The assumed coordinate system. Here, B0 is the magnetic field, n0 is the background electron density, and v is the ZF velocity.

Close modal

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(21)φ. Then, Eq. (1) can be written in the following dimensionless form:

wt+(ẑ×φ)·w+βφx=0,
(2)

where βV*/cs is treated as a (positive) constant. Let us also introduce the zonal average as f0Lxfdx/Lx, where Lx 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:

wt+(ẑ×φ)·w+βφx=0,
(3)
w=(2â)φ,
(4)

where â is an operator such that â=1 for DWs and â=0 for ZFs.39,40 Equations (3) and (4) constitute the so-called gHME.41 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,

w̃t+Uw̃x+(βU)φ̃x=0,
(5)
w̃=(21)φ̃.
(6)

Let us search for a solution in the form φ̃=ϕ(y)exp(ikxxiωt) with constant kx and ω. Then, w̃=(2/y2kx21)φ̃, and

(d2dy2α2UβUC)ϕ(y)=0,
(7)

where α21+kx2 and Cω/kx. This equation is identical42 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

A1(UA+βU)ϕ=Cϕ,
(8)

Ad2/dy2α2, so it is understood as an eigenvalue problem (with certain boundary conditions). Without loss of generality, we assume kx > 0. Hence, if Im C >0, then ωkxC has a positive imaginary part, which signifies an instability.

Also note that the complex conjugate of Eq. (7) is

(d2dy2α2UβUC*)ϕ*(y)=0.
(9)

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

U=u0cosqy.
(10)

(This periodic profile, although different from those typically used in studies of neutral fluids28–30,43,44 with fixed boundaries, is often used in numerical or theoretical modeling in plasma physics.23,45,46) For clarity, we assume u0 > 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 form47 

ϕ(y)=ψ(y)eiq¯y,
(11)

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/2q¯<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

ϕ*ϕϕϕ*(UβUCUβUC*)|ϕ|2=0.
(12)

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

0L(UβUCUβUC*)|ϕ|2dy=0L(ϕ*ϕϕϕ*)dy=(ϕ*ϕϕϕ*)|0L=0,
(13)

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

0LUβ|UC|2|ϕ|2dy=0.
(14)

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/dyUβ, 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=β means that

β<q2u0.
(15)

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

0Lϕ*ϕdyα20L|ϕ|2dy=0LUβUC|ϕ|2dy=0L(Uβ)(UC*)|UC|2|ϕ|2dy.
(16)

From Eq. (11), we have

0Lϕ*ϕdy=(ϕ*ϕ)|0L0L|ϕ|2dy=(ψ*ψ+iq¯|ψ|2)|0L0L|ϕ|2dy=0L|ϕ|2dy.
(17)

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

0L|ϕ|2dyα20L|ϕ|2dy<0.
(18)

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

0L(Uβ)(UU*)|UC|2|ϕ|2dy<0.
(19)

This means that (Uβ)(UU*) 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 –u0 ≤ C ≤ u0, which requires a special treatment due to possible singularities at U = C. (The case when C < −u0 or C > u0 will be considered in Sec. III C.) Near a singularity y = ys, two linearly independent solutions ϕ1,2 of Eq. (7) can be obtained using the Frobenius method28 and have the following asymptotics:

ϕ1=(yys)+a2(yys)2+a3(yys)3+,
(20)
ϕ2=1+b1(yys)+b2(yys)2++Gsϕ1ln(yys),
(21)

where Gs(Usβ)/Us. (The subscript s means that the corresponding functions are evaluated at y = ys.) The general solution can be written as

ϕ=c1ϕ1+c2ϕ2.
(22)

Near the singularity, one has ϕc2,ϕc2Gsln(yys), and ϕc2Gs/(yys).

Note that the above two solutions are invalid if Us=0, but those modes have infinite enstrophy and are physically irrelevant.50 Modes with nonzero Gs are physically irrelevant too for the same reason. The only exception is when c2 = 0, i.e., ϕ = c1ϕ1 = 0 at the singularity. However, this is impossible, which is seen as follows. Let us compare Eq. (7) with the following equation:

F=UUCF.
(23)

Between two neighboring singularities, we have C > U and α2 + β/(CU) > 0; hence, ϕ oscillates slower than F according to the Sturm comparison theorem.28 However, F = UC is a solution that is zero at both singularities, which means that ϕ cannot be zero at both singularities. This rules out the possibility that c2 = 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 Gs = 0 need be considered. This corresponds to C=Cnβ/q2 for a sinusoidal ZF. (The subscript n stands for “neutral”.) Let us substitute Cn into Eq. (7) along with Eq. (11). This gives

ψ+2iq¯ψq¯2ψ=(α2+UβUCn)ψ=(α2+q2u0cosqyβu0cosqy+β/q2)ψ=(α2q2)ψ.
(24)

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

(λ+q¯)2=q2α2.
(25)

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

α2=q2(q¯+mq)2.
(26)

For clarity, suppose q¯>0. (The case q¯<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

αn12=q2(qq¯)2,ψn1=eiqy,
(27)
αn22=q2q¯2,ψn2=1.
(28)

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

Now, let us consider α2=αn2+Δα2, where Δα2 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 = Cn + ΔC and ψ = ψn + Δψ. The neutral eigenmode ψn satisfies Eq. (24), while ψ is governed by

ψ+2iq¯ψq¯2ψ=(αn2+Δα2+UβUCnΔC)ψ.
(29)

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

ψψn*ψψn*+2iq¯(ψψn*+ψψn*)=(Δα2+UβUCnΔCUβUCn)ψψn*.
(30)

Let us integrate this equation over y from 0 to L. Using the periodicity of ψn and ψ, we obtain

0L(Δα2+UβUCnΔCUβUCn)ψψn*dy=(ψψn*ψψn*+2iq¯ψψn*)|0L=0.
(31)

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

Δα2J(ΔC)N,
(32)

where N0L|ψn|2dy and

J(ΔC)0L(UβUCnΔCUβUCn)|ψn|2dy.
(33)

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 = Cn; hence, the Taylor expansion in ΔC is invalid there. The correct way to calculate J is as follows. Since JC = 0) = 0, to the first order we have

J(ΔC)=J(0)ΔC.
(34)

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

ΔC=NΔα2J(0).
(35)

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 Cn is not equal to –β/q2 but lies out of the range of U (i.e., Cn < –u0 or Cn > u0), then J(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 Cn = –β/q2. 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

J(ΔC)=0LUβ(UCnΔC)2dy
(36)

at ΔC → 0. Since the integrand has a first-order pole at U = Cn when ΔC = 0, it is convenient to use the Plemelj formula

limϵ0+1xiϵ=P.V.(1x)+iπδ(x),
(37)

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

limΔC0Uβ(UCnΔC)2=P.V.Uβ(UCn)2iπδ(yy1)·(UU2)|y1sgn[Im(ΔC)]+iπδ(yy2)·(UU2)|y2sgn[Im(ΔC)],
(38)

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

limΔC0J(ΔC)=E2iDsgn[Im(ΔC)],
(39)

where

E=P.V.0LUβ(UCn)2dy=q2u0P.V.0L1cosqy+β/q2u0dy
(40)

and

D=π(UU2)|y1=π(UU2)|y2.
(41)

A straightforward calculation gives

E=0,
(42)
D=πqu01ϱ2>0,
(43)

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

ΔC=±iu0|Δα2|q21ϱ2
(44)

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 

FIG. 2.

A schematic of the ZF profile U(y) [Eq. (10)] and the locations of y1 and y2, where U = Cn. Note that U(y1)<0 and U(y2)>0.

FIG. 2.

A schematic of the ZF profile U(y) [Eq. (10)] and the locations of y1 and y2, where U = Cn. Note that U(y1)<0 and U(y2)>0.

Close modal

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

γ1,2=|kxu0|q2G(Δα1,22)1ϱ2.
(45)

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

Δα12=α2αn12,Δα22=α2αn22,
(46)

where αn12 and αn22 are given by Eqs. (27) and (28). As mentioned in Sec. III, we have 0αn12<αn22q2; 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, γ2γ1 everywhere too. Correspondingly, the largest TI growth rate is given by

γTI=|kxu0Δα22|q21ϱ2
(47)

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

ϱ2>1,α2<q2q¯2.
(48)

Finally, since |Δα22|=|α2αn22|=q2q¯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 nn2, we obtain the characteristic wavenumber of ψ(y) as

|ky|ψ/ψΔα22=q2α2<q.
(49)

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, α2=kx2; then, α2 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=q2q¯2. Strictly speaking, this conclusion is only valid when α2 is close to αn22. But one can also prove that no unstable eigenmode exists if α2q2. 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

0L|ϕ|2dy+α20L|ϕ|2dy=0L(Uβ)(UC*)|UC|2|ϕ|2dy.
(50)

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

0L|ϕ|2dy+α20L|ϕ|2dy=0L(Uβ)(U+CnCC*)|UC|2|ϕ|2dy=0L(Uβ)(UCn)(UCn)(U+Cn2Cr)|UC|2|ϕ|2dy=0L(Uβ)(UCn)(UCr)2(CnCr)2|UC|2|ϕ|2dy,
(51)

where CrReC. Since U=u0cosqy and Cn = –β/q2, one finds that (Uβ)/(UCn)=q2 is a positive constant. Hence

0L|ϕ|2dy+α20L|ϕ|2dy=q20L(UCr)2(CnCr)2|UC|2|ϕ|2dy<q20L|ϕ|2dy,
(52)

so we obtain

α2<q20L|ϕ|2dy0L|ϕ|2dyq2.
(53)

This shows that no unstable eigenmode is possible if α2q2. In application to the gHME, where α21+kx2, this indicates that not only do we have |ky|<q but the ZF itself must also have q2α2 > 1 for the TI onset. In other words, a ZF is always stable in the GO limit, because in this limit, we have q2 ≪ 1, and thus, α2 > q2 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

ϕ=n=N+Nϕnei(q¯+nq)y,
(54)

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 < q2u0 = 2.56, which indicates that ϱ > 1 is only a necessary condition.

FIG. 3.

(a) and (b) Im(ΔC) as a function of β for (a) α2 = 2.21 and (b) α2 = 1.16. (c) and (d) Im(ΔC) as a function of α2 for (c) β = 0.5 and (d) β = 1.7. (e) and (f) Im(ΔC) as a function of α2 for (e) q¯=0.1 and (f) q¯=0.6. For (a)–(d), the parameters are q =1.6, q¯=0, and u0 = 1. For (e) and (f), the parameters are q =1.6, β = 1, and u0 = 1. The solid curves are calculated using Eq. (44), while the dotted curves are numerical results obtained by solving Eq. (8) with N =50. Only positive Im(ΔC) values are shown, but the complete plots are symmetric with respect to the horizontal axis. This is because if C is a solution, then so is C* (Sec. II B). [Associated dataset available at https://doi.org/10.5281/zenodo.1241546.]52 

FIG. 3.

(a) and (b) Im(ΔC) as a function of β for (a) α2 = 2.21 and (b) α2 = 1.16. (c) and (d) Im(ΔC) as a function of α2 for (c) β = 0.5 and (d) β = 1.7. (e) and (f) Im(ΔC) as a function of α2 for (e) q¯=0.1 and (f) q¯=0.6. For (a)–(d), the parameters are q =1.6, q¯=0, and u0 = 1. For (e) and (f), the parameters are q =1.6, β = 1, and u0 = 1. The solid curves are calculated using Eq. (44), while the dotted curves are numerical results obtained by solving Eq. (8) with N =50. Only positive Im(ΔC) values are shown, but the complete plots are symmetric with respect to the horizontal axis. This is because if C is a solution, then so is C* (Sec. II B). [Associated dataset available at https://doi.org/10.5281/zenodo.1241546.]52 

Close modal

Next, we consider the change of Im(ΔC) with α2 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 Δα2 > 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 α2 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 α2 for nonzero q¯. Then, αn12=q2q¯2<q2 and αn22=q2(qq¯)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 α2 is close to αn1,22. It also captures the qualitative dependence on β. However, the dependence on α2 away from α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.

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 α2=1+kx2, so α2 > 1. Hence, the results are as follows. The two necessary conditions for the TI are (i) q2u0 > β and (ii) q2 > 1. The second condition comes from the fact that if q2 ≤ 1, then α21+kx2q2 for any real kx, 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 kx, such that 1+kx2 is smaller but close to q2. 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

q2>qmin2max{β/u0,1},
(55)

as also illustrated in Fig. 4(a).

FIG. 4.

(a) The analytically calculated TI-stability diagram [Eq. (55)] for a sinusoidal ZF within the gHME model. (b) A similar diagram found numerically for non-sinusoidal ZFs given by Eq. (58) with b = ±0.5, so feff = 2.2 (see Fig. 5). In both cases, β = 1. [Associated dataset available at https://doi.org/10.5281/zenodo.1241546.]52 

FIG. 4.

(a) The analytically calculated TI-stability diagram [Eq. (55)] for a sinusoidal ZF within the gHME model. (b) A similar diagram found numerically for non-sinusoidal ZFs given by Eq. (58) with b = ±0.5, so feff = 2.2 (see Fig. 5). In both cases, β = 1. [Associated dataset available at https://doi.org/10.5281/zenodo.1241546.]52 

Close modal

As discussed in Sec. III C, the TI growth rate γTI|kxIm(ΔC)| is found to be

γTI=|kxu0|ϑH(ϑ)1ϱ2.
(56)

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>ρs2 and (ii) q2u0>ρs2V*, and γTI is still given by Eq. (56), where

ϑ=1(q¯2+ρs2+kx2)/q2,ϱ=q2ρs2u0/V*.
(57)

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ρs2, is not adequate to describe the TI. This conclusion applies to both the traditional WKE1,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>ρs2, 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=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,

U(y)=u0k=0+Ukcos[(2k+1)qy],
(58)

where Uk = bk/(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

qeffk=0+qk2|Uk|k=0+|Uk|,
(59)

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

FIG. 5.

(a) Two examples of non-sinusoidal ZFs [Eq. (58)] with b = ±0.5. (b) The factor feffqeff/q, where the effective wavenumber qeff is given by Eq. (59).

FIG. 5.

(a) Two examples of non-sinusoidal ZFs [Eq. (58)] with b = ±0.5. (b) The factor feffqeff/q, where the effective wavenumber qeff is given by Eq. (59).

Close modal

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 u0 (denoted by u0,min) below which unstable eigenvalues disappear for all kx. For example, Fig. 6(a) shows the results for b = –0.5. It is seen that ϱeff,minqeff2u0,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 ϱeffqeff2u0/β can be considered as an effective RK parameter for non-sinusoidal ZFs; i.e., the first TI criterion holds approximately in the form

ϱeff1.
(60)
FIG. 6.

(a) qeff2u0,min versus q. The effective wavenumber qeff is given by Eq. (59). The ZF profile is given by Eq. (58) with b = –0.5. (b) qminfeff versus b for different u0, where feff is given in Fig. 5(b) and β = 1 is fixed. [Associated dataset available at https://doi.org/10.5281/zenodo.1241546.]52 

FIG. 6.

(a) qeff2u0,min versus q. The effective wavenumber qeff is given by Eq. (59). The ZF profile is given by Eq. (58) with b = –0.5. (b) qminfeff versus b for different u0, where feff is given in Fig. 5(b) and β = 1 is fixed. [Associated dataset available at https://doi.org/10.5281/zenodo.1241546.]52 

Close modal

Next, we consider a fixed b and numerically search for the minimum q (denoted by qmin) below which unstable eigenvalues disappear for all kx. We find that for nonzero b, qmin is smaller than one, and qmin approaches zero as |b| approaches one. However, qeff,minfeffqmin remains of order one. An example with β = 1 is given in Fig. 6(b), where we show different results depending on the choice of u0. In particular, consider b =0, which corresponds to a sinusoidal ZF (hence feff = 1). Then, for u0 = 0.5 < 1, Eq. (55) gives qmin=21.4, whereas for u0 ≥ 1, one obtains qmin = 1. Also, the u0-dependence becomes weak at u0 ≥ 3. The results show that qeff,min decreases slowly at nonzero b, and qeff,min0.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

qeff1.
(61)

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

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.

1.
P. H.
Diamond
,
S.-I.
Itoh
,
K.
Itoh
, and
T. S.
Hahm
,
Plasma Phys. Controlled Fusion
47
,
R35
(
2005
).
2.
A.
Fujisawa
,
Nucl. Fusion
49
,
013001
(
2009
).
3.
Z.
Lin
,
T. S.
Hahm
,
W. W.
Lee
,
W. M.
Tang
, and
R. B.
White
,
Science
281
,
1835
(
1998
).
4.
H.
Biglari
,
P. H.
Diamond
, and
P. W.
Terry
,
Phys. Fluids B: Plasma Phys.
2
,
1
(
1990
).
5.
W.
Dorland
,
F.
Jenko
,
M.
Kotschenreuther
, and
B. N.
Rogers
,
Phys. Rev. Lett.
85
,
5579
(
2000
).
6.
F.
Jenko
,
W.
Dorland
,
M.
Kotschenreuther
, and
B. N.
Rogers
,
Phys. Plasmas
7
,
1904
(
2000
).
7.
A. I.
Smolyakov
,
P. H.
Diamond
, and
V. I.
Shevchenko
,
Phys. Plasmas
7
,
1349
(
2000
).
8.
A. I.
Smolyakov
,
P. H.
Diamond
, and
M.
Malkov
,
Phys. Rev. Lett.
84
,
491
(
2000
).
9.
K.
Srinivasan
and
W. R.
Young
,
J. Atmos. Sci.
69
,
1633
(
2012
).
10.
J. B.
Parker
and
J. A.
Krommes
,
Phys. Plasmas
20
,
100703
(
2013
).
11.
J. B.
Parker
and
J. A.
Krommes
,
New J. Phys.
16
,
035006
(
2014
).
12.
B. F.
Farrell
and
P. J.
Ioannou
,
J. Atmos. Sci.
64
,
3652
(
2007
).
13.
B. F.
Farrell
and
P. J.
Ioannou
,
Phys. Plasmas
16
,
112903
(
2009
).
14.
N. C.
Constantinou
,
B. F.
Farrell
, and
P. J.
Ioannou
,
J. Atmos. Sci.
73
,
2229
(
2016
).
15.
M. A.
Malkov
,
P. H.
Diamond
, and
A.
Smolyakov
,
Phys. Plasmas
8
,
1553
(
2001
).
16.
K.
Miki
,
P. H.
Diamond
,
Ö. D.
Gürcan
,
G. R.
Tynan
,
T.
Estrada
,
L.
Schmitz
, and
G. S.
Xu
,
Phys. Plasmas
19
,
092306
(
2012
).
17.
P. H.
Diamond
,
Y.-M.
Liang
,
B. A.
Carreras
, and
P. W.
Terry
,
Phys. Rev. Lett.
72
,
2565
(
1994
).
18.
E.-J.
Kim
and
P. H.
Diamond
,
Phys. Rev. Lett.
90
,
185006
(
2003
).
19.
H.
Zhu
,
Y.
Zhou
,
D. E.
Ruiz
, and
I. Y.
Dodin
,
Phys. Rev. E
97
,
053210
(
2018
).
20.
B. N.
Rogers
,
W.
Dorland
, and
M.
Kotschenreuther
,
Phys. Rev. Lett.
85
,
5336
(
2000
).
21.
B. N.
Rogers
and
W.
Dorland
,
Phys. Plasmas
12
,
062511
(
2005
).
22.
E.-J.
Kim
and
P. H.
Diamond
,
Phys. Plasmas
9
,
4530
(
2002
).
23.
R.
Numata
,
R.
Ball
, and
R. L.
Dewar
,
Phys. Plasmas
14
,
102312
(
2007
).
24.
D. A.
St-Onge
,
J. Plasma Phys.
83
,
905830504
(
2017
).
25.
R.
Singh
,
H.
Jhang
, and
H. K.
Kaang
,
Phys. Plasmas
23
,
074505
(
2016
).
26.
F.
Rath
,
A. G.
Peeters
,
R.
Buchholz
,
S. R.
Grosshauser
,
F.
Seiferling
, and
A.
Weikl
,
Phys. Plasmas
25
,
052102
(
2018
).
27.
J. B.
Parker
, Ph.D. thesis, Princeton University,
2014
.
29.
F. B.
Lipps
,
J. Fluid Mech.
12
,
397
(
1962
).
30.
P. G.
Drazin
and
L. N.
Howard
,
Adv. Appl. Mechanics
9
,
1
(
1966
).
31.
P. G.
Drazin
and
W. H.
Reid
,
Hydrodynamic Stability
(
Cambridge University Press
,
1981
).
32.
J. B.
Parker
,
J. Plasma Phys.
82
,
595820602
(
2016
).
33.
A. I.
Smolyakov
and
P. H.
Diamond
,
Phys. Plasmas
6
,
4410
(
1999
).
34.
J. C.
Li
and
P. H.
Diamond
,
Phys. Plasmas
25
,
042113
(
2018
).
35.
M. J.
Pueschel
,
B. J.
Faber
,
J.
Citrin
,
C. C.
Hegna
,
P. W.
Terry
, and
D. R.
Hatch
,
Phys. Rev. Lett.
116
,
085001
(
2016
).
36.
D. R.
Hatch
,
F.
Jenko
,
A.
Bañon Navarro
,
V.
Bratanov
,
P. W.
Terry
, and
M. J.
Pueschel
,
New J. Phys.
18
,
075018
(
2016
).
37.
A. E.
Fraser
,
P. W.
Terry
,
E. G.
Zweibel
, and
M. J.
Pueschel
,
Phys. Plasmas
24
,
062304
(
2017
).
38.
A.
Hasegawa
and
K.
Mima
,
Phys. Rev. Lett.
39
,
205
(
1977
).
39.
W. D.
Dorland
, Ph.D. thesis, Princeton University,
1993
.
40.
G. W.
Hammett
,
M. A.
Beer
,
W.
Dorland
,
S. C.
Cowley
, and
S. A.
Smith
,
Plasma Phys. Controlled Fusion
35
,
973
(
1993
).
41.
J. A.
Krommes
and
C.-B.
Kim
,
Phys. Rev. E
62
,
8508
(
2000
).
42.

Equation (7) is also identical to Eq. (18) in Ref. 23 except that x and y are swapped because the notations are different.

43.
P. G.
Drazin
and
L. N.
Howard
,
J. Fluid Mech.
14
,
257
(
1962
).
44.
L. N.
Howard
and
P. G.
Drazin
,
Stud. Appl. Math.
43
,
83
(
1964
).
45.
S.
Kobayashi
and
B. N.
Rogers
,
Phys. Plasmas
19
,
012315
(
2012
).
46.
D. E.
Ruiz
,
J. B.
Parker
,
E. L.
Shi
, and
I. Y.
Dodin
,
Phys. Plasmas
23
,
122304
(
2016
).
47.
V. A.
Yakubovich
and
V. M.
Starzhinskii
,
Linear Differential Equations with Periodic Coefficients
(
Keter Publishing House Jerusalem Ltd.
,
1975
), Vol. 1.
48.
C. C.
Lin
,
The Theory of Hydrodynamic Stability
(
Cambridge University Press
,
1955
).
49.
J.
Pedlosky
,
Geophysical Fluid Dynamics
(
Springer-Verlag
,
Berlin
,
1982
), Chap. 7.
50.

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

51.

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.

52.
H.
Zhu
,
Y.
Zhou
, and
I. Y.
Dodin
, “
On the Rayleigh–Kuo criterion for the tertiary instability of zonal flows
,”
Zenodo
. .
53.

The second part of our instability criterion, q2>ρs2, 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, where the non-adiabatic electron response is assumed. Therefore, our α2=1+kx2 must be replaced with α2=kx2. Then, α2<q2 can also be satisfied even within the GO assumption.

54.
H.
Zhu
,
Y.
Zhou
, and
I. Y.
Dodin
,
Phys. Plasmas
25
,
072121
(
2018
).