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.
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 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, , where is the corresponding electrostatic potential on the two-dimensional plane . The electron response to E is adiabatic (yet see below), while the ion response can be described by the drift and the polarization drift. Then, assuming the quasi-neutrality condition, the evolution of δφ is described by the (original) Hasegawa–Mima equation
Here, is the ion sound radius (we use to denote definitions), is the ion sound speed, Z is the ion charge number, is the ion gyrofrequency, e is the electron charge, is the velocity, is the unit vector along the z axis, is the electron diamagnetic drift velocity, and is the characteristic length scale of n0. Also, is the Laplacian.
Let us measure time in units 1/Ωi and length in units ρs. Let us also introduce a normalized potential and a normalized “generalized vorticity” . Then, Eq. (1) can be written in the following dimensionless form:
where is treated as a (positive) constant. Let us also introduce the zonal average as , 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:
B. Reduction to an eigenvalue problem
Consider a stationary ZF with and . (Hereafter, the prime denotes the derivative with respect to y.) The ZF velocity is , where and 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 with constant kx and ω. Then, , and
where and . 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
, 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
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 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
where ψ(y) is periodic such that for any y, and is a constant. Assuming that ϕ is bounded (as it would be the case, for instance, at periodic boundary conditions), must be real. Without loss of generality, we limit the value of to the first Brillouin zone, i.e., .
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, must change its sign somewhere in the integration domain; i.e., there must be a location where . (This location can be understood as the point where the “vorticity” Z, defined via , 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 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
This means that 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 –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:
where . (The subscript s means that the corresponding functions are evaluated at y = ys.) The general solution can be written as
Near the singularity, one has , and .
Note that the above two solutions are invalid if , 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:
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 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 for a sinusoidal ZF. (The subscript n stands for “neutral”.) Let us substitute Cn 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 . 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 . (The case 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 , 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
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 with . Therefore
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 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 Cn is not equal to –β/q2 but lies out of the range of U (i.e., Cn < –u0 or Cn > u0), then is real. Accordingly, Eq. (35) gives that ΔC is real, because Δα2 and 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.
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
where P.V. stands for the principal value. Then, we have
where we used and (Fig. 2). Hence, we obtain
A straightforward calculation gives
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 , namely,
Here, , where H is the Heaviside step function, and
where and are given by Eqs. (27) and (28). As mentioned in Sec. III, we have ; hence, . If , then ; hence, ; if , then ; hence, we still have . This shows that everywhere, and thus, γ2 ≥ γ1 everywhere too. Correspondingly, the largest TI growth rate is given by
at , and otherwise γTI = 0. Accordingly, the TI develops when
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, ; then, α2 is allowed to be zero, so Eq. (49) is less stringent, namely, .
D. There is no instability at
Since , the above results indicate that there is no instability at . Strictly speaking, this conclusion is only valid when α2 is close to . But one can also prove that no unstable eigenmode exists if . 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 . Since and Cn = –β/q2, one finds that is a positive constant. Hence
so we obtain
This shows that no unstable eigenmode is possible if α2 ≥ q2. In application to the gHME, where , this indicates that not only do we have 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).
IV. COMPARISON WITH NUMERICAL CALCULATIONS
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 , 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.
Next, we consider the change of Im(ΔC) with α2 in Figs. 3(c) and 3(d). Since , the two neutral eigenmodes are at and . We also see that there are no unstable eigenvalues at Δα2 > 0, which agrees with the perturbation theory. Unstable eigenvalues exist when , whose numerical values agree with perturbation theory when . 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 . In Figs. 3(e) and 3(f), we plot the dependence of Im(ΔC) on α2 for nonzero . Then, and . In Fig. 3(e), is small, and the existence of Im(ΔC) > 0 agrees with perturbation theory. In Fig. 3(f), is large, so and are close to each other. Thus, the branch of unstable eigenvalues starting from intersects with the branch starting from .
In summary, our first-order perturbation theory agrees with numerical results when α2 is close to . It also captures the qualitative dependence on β. However, the dependence on α2 away from 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 , 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 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 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
as also illustrated in Fig. 4(a).
As discussed in Sec. III C, the TI growth rate is found to be
Here, H is the Heaviside step function, , and . In a dimensional form, the two necessary conditions found above are (i) and (ii) , and γTI is still given by Eq. (56), where
As a corollary, the WKE, which assumes the GO limit that relies on the assumption , 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, , 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 . 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 . Also, at arbitrary β, our Eq. (57) shows that when q is large enough, one has , 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 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
where . We also define [Fig. 5(b)], which increases as 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 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 remains approximately constant and close to one, namely, ϱeff,min ≈ 1.2. Similar results are obtained for other choices of b within the range . Thus, we conclude that 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 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 approaches one. However, 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 , 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 within the range . 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;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.
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, , 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 , where the non-adiabatic electron response is assumed. Therefore, our must be replaced with . Then, can also be satisfied even within the GO assumption.