The phase space of driftons (drift-wave quanta) is studied within the generalized Hasegawa–Mima collisionless-plasma model in the presence of zonal flows. This phase space is made intricate by the corrections to the drifton ray equations that were recently proposed by Parker [J. Plasma Phys. 82, 595820602 (2016)] and Ruiz et al. [Phys. Plasmas 23, 122304 (2016)]. Contrary to the traditional geometrical-optics (GO) model of the drifton dynamics, it is found that driftons can not only be trapped or passing but also accumulate spatially while experiencing indefinite growth of their momenta. In particular, it is found that the Rayleigh–Kuo threshold known from geophysics corresponds to the regime when such “runaway” trajectories are the only ones possible. On one hand, this analysis helps to visualize the development of the zonostrophic instability, particularly its nonlinear stage, which is studied here both analytically and through wave-kinetic simulations. On the other hand, the GO theory predicts that zonal flows above the Rayleigh–Kuo threshold can only grow; hence, the deterioration of intense zonal flows cannot be captured within a GO model. In particular, this means that the so-called tertiary instability of intense zonal flows cannot be adequately described within the quasilinear wave kinetic equation, contrary to some previous studies.

The interaction between zonal flows (ZFs) and drift-wave (DW) turbulence has a substantial effect on turbulent transport in fusion devices, and hence, has been actively studied in plasma physics for few decades.1–7 One of the popular reduced models of DW dynamics in ZFs is the wave kinetic equation (WKE),8–17 which assumes the geometrical-optics (GO) approximation, i.e., loosely speaking, that the DW wavelengths are vanishingly small compared to ZF scales. Within this approximation, DWs can be understood as a gas of “driftons,” which are quasi-particles described by coordinates x, momenta p (DW wave vectors), and energies H (DW frequencies). In particular, H=H(t,x,p) serves as the quasi-particle Hamiltonian that determines the drifton trajectory for given x and p at time t. This model facilitates understanding of many important effects, including the zonostrophic instability (ZI), i.e., the formation of ZF out of DW turbulence.16–20 

It was shown recently that the drifton Hamiltonian used in previous studies is oversimplified, and an improved Hamiltonian has been proposed in Refs. 21 and 22 based on the generalized Hasegawa–Mima equation (gHME).23 The corresponding improved WKE (iWKE) accounts for the loss of drifton enstrophy to ZFs, and hence, is a more adequate GO model. The advantages of the iWKE are demonstrated in Refs. 21 and 22 by numerical simulations. A numerical comparison between the iWKE and the quasilinear gHME was reported in Ref. 24. However, it is also insightful to explore the single-particle drifton dynamics, i.e., the drifton phase-space trajectories in a prescribed ZF. Such a study can help elucidate the importance of individual terms in the drifton Hamiltonian. It can also help understand the nonlinear dynamics of ZFs, including the nonlinear stage of the ZI, and identify factors that are important for its saturation. Here, we report such a study, which explores in depth the drifton dynamics within the iWKE proposed in Refs. 21 and 22. Some of our results were also highlighted in Ref. 25. The purpose of this paper is to expand the discussion and to elaborate on details.

Our main findings are as follows: (i) Contrary to the traditional WKE (tWKE) of the drifton dynamics, which predicts12–14 nonlinear structures à la Bernstein–Greene–Kruskal (BGK) waves,26 the iWKE predicts that driftons do not have to be just passing or trapped. Instead, they can accumulate in certain spatial locations while experiencing indefinite growth of their momenta. We call such trajectories “runaway.” (ii) Depending on the ZF parameters, the drifton phase space can have three different regimes. In regime 1, passing, trapped, and runaway trajectories coexist. In regime 2, passing trajectories disappear entirely, but both trapped and runaway trajectories can coexist. In regime 3, only runaway trajectories are left. (iii) Remarkably, regime 3 is precisely the regime when the ZF amplitude exceeds the Rayleigh–Kuo threshold known from geophysics.27 Also notably, this regime is not captured by the tWKE. (iv) We apply our phase-space analysis to visualize the development of the ZI, particularly its nonlinear stage, using both theoretical arguments and iWKE simulations. Moreover, we find that the GO theory predicts that ZFs above the Rayleigh–Kuo threshold can only grow; hence, the deterioration of intense ZFs cannot be captured within a GO model.25,28 In particular, this means that the so-called tertiary instability (TI) of intense ZFs cannot be adequately described within a quasilinear WKE (including the tWKE and the iWKE, which both assume the GO limit), contrary to some previous studies. Our results serve as a stepping stone toward revising basic physics of DW–ZF interactions from the new perspective of drifton phase-space dynamics beyond the traditional (tWKE-based) approach.

The rest of this paper is organized as follows. In Sec. II, the gHME and the iWKE are introduced. In Sec. III, the three different regimes of the drifton phase-space structure are described. The two critical ZF magnitudes that separate these three regimes are also given. In Secs. IV and V, the nonlinear ZI and the TI are discussed, respectively. Our main conclusions are summarized in Sec. VI. Auxiliary calculations are given in  Appendix.

First, let us introduce the original Hasegawa–Mima equation.29 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

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 of 1/Ωi and length in units of ρ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. Equation (2) represents the original Hasegawa–Mima equation.

Let us 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 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 DWs. To account for this and thus make the plasma model more realistic, the governing equations can be rewritten as follows:

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

where â is an operator such that â=1 for DWs and â=0 for ZFs.30,31 Equations (3) and (4) constitute the so-called gHME,23 which is the model that we assume below.

The ZF can be described using the average velocity U(y,t)yφ, and DWs can be described using the zonal average of their Wigner function

W(y,p,t)d2seip·sw̃(x+s2,t)w̃(xs2,t),
(5)

which in the GO limit can be understood as the drifton phase-space distribution.22 (The discussion on the positive definiteness of the Wigner function in the concept of quantum mechanics can be found in Ref. 32.) The GO limit itself is defined as the regime where

ϵmax(λDWλZF,ρsλZF)1,
(6)

where λZF and λDW are the wavelengths of ZFs and DWs, respectively. To proceed, the quasilinear approximation is used,18,19,33 which assumes that the DW self-interactions can be ignored. (Recent work34 has also gone beyond the quasilinear approximation.) Then, the evolution equations for W and U are21,22

Wt={H,W}+2ΓW,
(7)
Ut=yd2p(2π)2pxpyWpD4,
(8)

where pD21+px2+py2 and {⋅, ⋅} is the canonical Poisson bracket, namely,

{A,B}Ax·BpAp·Bx.
(9)

The Hermitian and anti-Hermitian parts of the Hamiltonian are given by

H=βpx/pD2+pxU+pxU/pD2,
(10)
Γ=Upxpy/pD4,
(11)

where primes denote derivatives with respect to y.

Equation (7) is the iWKE as described in Refs. 21 and 22. In comparison, the tWKE used in previous studies is given by the same Eq. (7) but with different H and Γ, namely,

Ht=βpx/pD2+pxU,
(12)
Γt=0,
(13)

where the subscript “t” stands for “traditional.” The iWKE conserves the total enstrophy (per unit length in x) Ztotal = ZDW + ZZF and the total energy (per unit length in x) Etotal = EDW + EZF of the ZF–DW system, where

ZDW12d2pdy(2π)2W,
(14)
ZZF12dy(U)2,
(15)
EDW12d2pdy(2π)2WpD2,
(16)
EZF12dyU2.
(17)

In contrast, the tWKE conserves only the DW enstrophy but not the total enstrophy, and as a result, can be unsatisfactory in many respects.21,22

As an integral of the Wigner function over the whole phase space, the DW enstrophy [Eq. (15)] can be considered as the total number of driftons.22 Due to the presence of nonzero Γ [Eq. (11)], ZDW is not conserved, and hence, the iWKE [Eq. (7)] does not conserve the total number of driftons. However, it can be made conservative in the case of stationary U by introducing F(y,p,t)W(y,p,t)/[βU(y)]; then, Eq. (7) becomes21,35

Ft={H,F},
(18)

where H is still given by Eq. (10). This shows that F is conserved along drifton trajectories, which are given by Hamilton's equations

dydt=Hpy=2pxpypD4(βU),
(19)
dpydt=Hy=pxpD2(U+pD2U).
(20)

Equations (19) and (20) describe the drifton dynamics within the iWKE. Since H does not depend on x, px is also conserved along the trajectory. Therefore, the drifton dynamics can be studied on the (y, py) plane with px serving as a parameter.

In more general situations, when U is not stationary, Eq. (18) does not apply. But even in those situations, one can still view Eqs. (19) and (20) as equations of the drifton motion, while Γ only affects the evolution of the drifton density along such trajectories, not the trajectories themselves.

For simplicity, we assume a sinusoidal ZF, namely,

U(y)=u0cosqy,
(21)

where u0 and q are constant (for clarity, we assume u0 > 0 and q >0). Then, Eq. (10) leads to the following expression for the Hamiltonian:

H(y,py)=βpxpD2+pxu0cosqy(1q2pD2),
(22)

and Eqs. (19) and (20) become

dydt=2pxpypD4(β+q2u0cosqy),
(23)
dpydt=pxqu0(1q2pD2)sinqy.
(24)

Due to the assumed GO approximation, we limit our consideration to the regime where q2 ≪ 1 (in dimensional form, q2ρs2). We also assume that px is nonzero. Then, by studying the drifton phase-space trajectories governed by Eqs. (23) and (24), one can identify three distinct regimes depending on how the ZF magnitude u0 compares with the two critical values (the derivations are given in  Appendix)

uc,1β2q2,uc,2βq2.
(25)

The GO approximation implies 0 < uc,1uc,2. The phase-space structures are illustrated in Fig. 2 that shows typical contour plots of H corresponding to three distinct regimes. (In a stationary ZF considered here, driftons travel along constant-energy surfaces.) Specifically, these three regimes are as follows.

FIG. 2.

Contour plots of H given by Eq. (22). The arrows show the phase-space velocity fields (ẏ,pẏ) given by Eqs. (23) and (24). Three different regimes are shown: (a) regime 1, which corresponds to a weak ZF (u0 = 0.1); (b) regime 2, which corresponds to a moderate ZF (u0 = 2); and (c) regime 3, which corresponds to a strong ZF (u0 = 10). The labels “T,” “P,” and “R” denote trapped, passing, and runaway trajectories, respectively. The white dashed lines in (c) are |y|=y*, where y* is given by Eq. (29). In all cases, the parameters are β = 1, q =0.5, and px = 0.5.

FIG. 2.

Contour plots of H given by Eq. (22). The arrows show the phase-space velocity fields (ẏ,pẏ) given by Eqs. (23) and (24). Three different regimes are shown: (a) regime 1, which corresponds to a weak ZF (u0 = 0.1); (b) regime 2, which corresponds to a moderate ZF (u0 = 2); and (c) regime 3, which corresponds to a strong ZF (u0 = 10). The labels “T,” “P,” and “R” denote trapped, passing, and runaway trajectories, respectively. The white dashed lines in (c) are |y|=y*, where y* is given by Eq. (29). In all cases, the parameters are β = 1, q =0.5, and px = 0.5.

Close modal

Regime 1—The first regime corresponds to u0 < uc,1 (weak ZF). In this regime, there are two types of phase-space stationary points, namely, the stable stationary points (centers) at

y=±πq,py=0
(26)

and the unstable stationary point (saddle) at

y=0,py=0.
(27)

(Due to the periodicity of the system, we limit our consideration to a single period, y ∈ [–π/q, π/q].) The trajectories in this regimes are of three different types [Fig. 2(a)]: passing (labeled by “P”), trapped (labeled by “T”), and runaway (labeled by “R”). Passing trajectories reside near the saddle, while trapped trajectories reside near the centers. Passing and trapped trajectories are qualitatively similar to those predicted by the tWKE12–14 and are also reminiscent of the corresponding trajectories of charged particles interacting with plasma waves. If these were the only trajectories, a nonlinear ZF–DW system in this regime would have been similar to a BGK wave,26 but the runaway trajectories make the picture qualitatively different. These trajectories are localized spatially around y =0 but extend to infinity along the momentum axis; therefore, these driftons tend to accumulate in certain spatial locations. (This is understood from the fact that, while |py| remains growing, the DW group velocity at large |py| is ẏpy30, so a drifton eventually stops moving along y.) The existence of runaway driftons indicates that a ZF–DW system cannot be in an exact steady state if runaway trajectories are populated; otherwise F is a function of H only [Eq. (18)] and the runaway trajectories will make W non-integrable.

Note that runaway trajectories are possible even for arbitrarily small u0. Also note that runaway trajectories can also be obtained from the tWKE [Eqs. (12) and (13)]. This explains that some plots in Refs. 13 and 14 obtained from the tWKE appear similar to our Fig. 2(a), even though the underlying models are different.

Regime 2—The fraction of passing trajectories shrinks with the increase in u0. In the second regime, when uc,1u0uc,2 (moderate ZF), passing trajectories disappear entirely. This is illustrated in Fig. 2(b). The centers and the saddle are also given by Eqs. (26) and (27). The trapped trajectories reside near the center, while the remaining phase space corresponds to runaways.

Note that regime 2 is also possible in the tWKE, where the derivatives of U are omitted in H, so q2 is effectively set to zero. This leads to uc,1 = β/2, which remains finite. Therefore, passing trajectories also disappear entirely if u0uc,1. On the other hand, setting q2 to zero leads to uc,2 = +; hence, the regime 3 below is impossible in the tWKE.

Regime 3—The third regime corresponds to u0uc,2 (strong ZF). In this regime, all the stationary points [Eqs. (26) and (27)] are unstable, so trapped trajectories also disappear, and only runaway trajectories are left. As illustrated in Fig. 2(c), all the driftons move towards

|y|=y*,|py|=.
(28)

Here, ±y* are the locations where U=β; namely,

y*=1q[πarccos(βq2u0)].
(29)

In particular, note that no trajectory can cross the vertical lines |y|=y* [shown as vertical white dashed lines in Fig. 2(c)], since the drifton velocity is always zero at |y|=y* [Eq. (23)].

Remarkably, the condition under which regime 3 is realized (i.e., that u0>uc,2β/q2) is precisely the Rayleigh–Kuo criterion,27 which states that a necessary condition for the ZF instability is the existence of spatial locations where U=β. The connection between regime 3 and the Rayleigh–Kuo criterion is only captured by the improved Hamiltonian [Eq. (10)]. In the tWKE, where the U term in H is neglected, the Rayleigh–Kuo criterion is never reached, and hence regime 3 cannot be realized. However, as will be argued below (Sec. VI), regime 3 in the iWKE does not play quite the same role as that in the Rayleigh–Kuo criterion. In contrast with the full-wave theory, the GO model predicts that driftons in regime 3 amplify a ZF rather than destroy it, as we will now discuss.

Here, we study the nonlinear structures of DW turbulence in ZFs based on the drifton phase-space trajectories presented in Sec. III. Specifically, we consider the ZI, which describes the formation of ZFs out of DW turbulence with a given equilibrium drifton Wigner function W(p).16–20 Assuming perturbations of the form U=Re(Uqeiqy+γZIt) and δW=Re(Wqeiqy+γZIt), the linearized Eqs. (10) and (11) are

Wq=1γZI+2iβqpxpy/pD4×[iqpx(1q2pD2)Wpy+2iq3pxpypD4W],
(30)
Uq=iqγZId2p(2π)2pxpypD4Wq.
(31)

By plugging (30) into (31) and integrating by parts the term that contains W/py, we obtain the dispersion relation for the linear ZI of the iWKE21,25

1=d2p(2π)2q2px2pD4(14py2/pD2)(1q2/pD2)(γZIpD4+2iβqpxpy)2W(p).
(32)

As the ZF amplitude becomes finite, the ZI enters its nonlinear stage and eventually saturates. Previous studies discussed how the saturated state is determined by the interplay of ZFs and passing and trapped orbits.12–14,36 However, this picture is qualitatively altered by runaway trajectories. To show this, we numerically simulate the iWKE [Eqs. (10) and (11)] using the pseudo-spectral method described in Ref. 24 (where the iWKE is termed as “CE2-GO”). A weak 8th-order hyperviscosity is added for numerical stability.20 We launch the simulation with an initial Gaussian DW distribution

W(p)W(t=0,y,p)=4πW0r2exp(|p|22r2)
(33)

(where r is some constant serving as a characteristic DW wavenumber) and an initial ZF perturbation

U(t=0)=Uqcosqy
(34)

with small Uq. It is found that depending on the strength of the DW amplitude W0 [Eq. (33)], the ZF can saturate in one of the three different regimes described in Sec. III. The simulation results are shown in Figs. 3–5, where the structure of the DW Wigner functions W reflects the structure of the underlying drifton trajectories.

FIG. 3.

Results of iWKE simulations with initial conditions given by Eqs. (33) and (34). Here, the parameters are β = 1, q =0.5, W0 = 0.3, r =0.5, Uq = 0.01, and the hyperviscocity coefficient is ν = 1 × 10−7.20 (a) The drifton density n [Eq. (36)] and the ZF velocity U at t =60, 120, and 180. The relation between the change in n and U agrees with Eq. (40). (b) The time evolution of the energy E and enstrophy Z integrated over one spatial period y ∈ [–2π, 2π]. The energy and enstrophy exchange between DWs and ZF is small, since the ZF is weak. (c) The drifton phase-space Wigner function W(y, py) at px = 0.5 at the three different instants. Passing and trapped trajectories are clearly seen. [Associated dataset available at https://doi.org/10.5281/zenodo.1244318.]37 

FIG. 3.

Results of iWKE simulations with initial conditions given by Eqs. (33) and (34). Here, the parameters are β = 1, q =0.5, W0 = 0.3, r =0.5, Uq = 0.01, and the hyperviscocity coefficient is ν = 1 × 10−7.20 (a) The drifton density n [Eq. (36)] and the ZF velocity U at t =60, 120, and 180. The relation between the change in n and U agrees with Eq. (40). (b) The time evolution of the energy E and enstrophy Z integrated over one spatial period y ∈ [–2π, 2π]. The energy and enstrophy exchange between DWs and ZF is small, since the ZF is weak. (c) The drifton phase-space Wigner function W(y, py) at px = 0.5 at the three different instants. Passing and trapped trajectories are clearly seen. [Associated dataset available at https://doi.org/10.5281/zenodo.1244318.]37 

Close modal
FIG. 4.

The same as Fig. 3, except with W0 = 1 and ν = 2 × 10−7. In (b), the decrease in the total enstrophy at t40 is due to the fact that runaway driftons at large |py| are heavily damped by hyperviscosity. It is seen that the energy exchange between DWs and ZF is large due to runaways. In (c), trapped and runaway trajectories are clearly seen. [Associated dataset available at https://doi.org/10.5281/zenodo.1244318.]37 

FIG. 4.

The same as Fig. 3, except with W0 = 1 and ν = 2 × 10−7. In (b), the decrease in the total enstrophy at t40 is due to the fact that runaway driftons at large |py| are heavily damped by hyperviscosity. It is seen that the energy exchange between DWs and ZF is large due to runaways. In (c), trapped and runaway trajectories are clearly seen. [Associated dataset available at https://doi.org/10.5281/zenodo.1244318.]37 

Close modal
FIG. 5.

The same as Figs. 3 and 4, except with W0 = 50 and ν = 10−5. In (b), the decrease in the total enstrophy at t7 is due to the fact that runaway driftons at large |py| are heavily damped by hyperviscosity. It is seen that the energy exchange between DWs and ZF is large due to the runaways. Note that the ZF fails to absorb the whole DW energy, which is explained in the text. In (c), runaway trajectories are clearly seen. However, more intricate structures emerge too, because the ZF is far from sinusoidal. [Associated dataset available at https://doi.org/10.5281/zenodo.1244318.]37 

FIG. 5.

The same as Figs. 3 and 4, except with W0 = 50 and ν = 10−5. In (b), the decrease in the total enstrophy at t7 is due to the fact that runaway driftons at large |py| are heavily damped by hyperviscosity. It is seen that the energy exchange between DWs and ZF is large due to the runaways. Note that the ZF fails to absorb the whole DW energy, which is explained in the text. In (c), runaway trajectories are clearly seen. However, more intricate structures emerge too, because the ZF is far from sinusoidal. [Associated dataset available at https://doi.org/10.5281/zenodo.1244318.]37 

Close modal

At the nonlinear stage of ZI, the time-evolution of ZFs can be qualitatively estimated from the drifton trajectories. In order to demonstrate this, we derive the drifton hydrodynamic equations as follows. Let us write down Eq. (7) explicitly,

Wt=(UpxpD2+Upx)Wpy+2(Uβ)pxpypD4Wy2pxpyUpD4W.
(35)

Integrating Eq. (35) over p leads to the following equation for the drifton density n:

nt+(βU)Vy=0,
(36)

where

n(y,t)d2p(2π)2W(y,p,t)
(37)

and

V(y,t)d2p(2π)22pxpypD4W(y,p,t).
(38)

Then, one can rewrite Eq. (8) as follows:

Ut=12Vy,
(39)

and the combination of Eqs. (37) and (38) gives

Ut=12(Uβ)nt.
(40)

[Even though the denominator is zero at U=β, the right-hand side of Eq. (40) remains finite because, according to Eqs. (36), ∂n/∂t is also zero at such locations.]

By using the knowledge of the phase-space trajectories, which determine the drifton flows, one can predict whether the drifton density grows or decreases at a given location. This gives the sign of ∂n/∂t; then, the sign of ∂U/∂t can be inferred from Eq. (40), so one can tell whether the ZF is peaking or flattening. In regimes 1 and 2, Uβ is always negative; hence, ∂U/∂t and ∂n/∂t have opposite signs, which is consistent with Figs. 3 and 4. Regime 3 is more interesting due to its connection with the Rayleigh–Kuo criterion. Let us consider a ZF of the assumed sinusoidal form [Eq. (21)] that satisfies u0 > uc,2. Then, from Fig. 2(c), it is seen that within |y|<y*, driftons move away from y =0. Therefore, n decreases at y =0. Since Uβ<0 at y =0, from Eq. (40) we have

U(y=0)t>0.
(41)

A similar argument leads to

U(y=±π/q)t<0.
(42)

Hence, the ZF profile gets more peaked, i.e., the ZF is globally amplified. The amplification of the ZF, in turn, will reinforce the drifton runaway. From the expression of the DW energy [Eq. (16)], the bulk motion of driftons to |py|= causes pD2; hence, EDW → 0. Therefore, in regime 3, the ZF tends to absorb all the energy from DWs; the reasons why this does not happen in Fig. 5 are twofold: (i) in order to better visualize drifton trajectories, we have chosen a Gaussian initial distribution centered at p=0 [Eq. (33)]; hence, a large fraction of driftons with px ≪ 1, which move in phase space slowly [Eqs. (19) and (20)], does not significantly participate in the energy exchange; also, (ii) some driftons do not runaway as can be seen in Fig. 5(c), since the ZF is far from sinusoidal in this highly nonlinear stage.

In addition to the ZI of a ZF–DW system, which leads to the ZF amplification, it is also of interest to examine whether the iWKE is applicable to describe the so-called tertiary instability (TI).10,38–43 Specifically, we define the TI as the instability of a DW on top of a prescribed non-turbulent ZF equilibrium, i.e., an instability of a Kelvin–Helmholtz type. Note that this definition is different from that in Refs. 39 and 40, where the TI was attributed to the ion-temperature gradient (absent in our model), but similar to that in the majority of relevant papers.10,41–43

As speculated in Refs. 21 and 41 and later elaborated in Refs. 25 and 28, the Rayleigh–Kuo criterion is a necessary condition for this instability, so one might expect the TI to develop in regime 3. However, as shown above, ZF can only grow in regime 3 rather than deteriorate, so in principle, there is no TI in the WKE under the GO assumption. The reason is that, as shown in Ref. 28, within the gHME, the TI requires q2 > 1, while a GO model relies on the assumption that q2 ≪ 1.

Another explanation for the absence of the TI in our model is as follows. Let us consider a small DW perturbation around a stationary ZF. Linearizing Eq. (7) gives

W1t={H0,W1}+2Γ0W1,
(43)

where W1 is the Wigner function of driftons. (The zeroth-order DW Wigner function is zero because, as mentioned earlier, we assume a non-turbulent background for the TI.) Since H0 and Γ0 are stationary, Eq. (43) is equivalent to

ddt(W1Uβ)=0
(44)

(as we mentioned in Sec. II C for a similar equation), where d/dt is taken along the drifton trajectories determined by Eqs. (19) and (20). As shown in  Appendix, drifton runaway trajectories do not reach locations where U=β. Hence, Uβ remains finite along trajectories, and W1 cannot grow exponentially with time. This rules out the TI and, accordingly, this also means that the TI cannot be described by the quasilinear WKE, because the WKE relies on the GO approximation. (However, full-wave quasilinear models are perfectly capable of capturing the TI; for example, see Refs. 25, 28, and 44.)

In summary, this paper presents the first study of the drifton phase-space dynamics within the iWKE proposed in Refs. 21 and 22. Contrary to the traditional GO model of the drifton dynamics, it is found that driftons can not only be trapped or passing, but also accumulate spatially while experiencing indefinite growth of their momenta. In particular, it is found that the Rayleigh–Kuo threshold known from geophysics corresponds to the regime when such “runaway” trajectories are the only ones possible. On one hand, this analysis helps to visualize the development of the ZI, particularly its nonlinear stage, which we study both analytically and through iWKE simulations. On the other hand, the GO theory predicts that ZFs above the Rayleigh–Kuo threshold can only grow; hence, the deterioration of intense ZFs cannot be captured within a GO model. In particular, this means that the so-called tertiary instability of intense zonal flows cannot be adequately described within the quasilinear WKE, contrary to some previous studies.

The authors thank J. B. Parker for providing a copy of his code for our wave-kinetic simulations. 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.

Here, we provide a detailed description of the drifton phase-space trajectories governed by Eqs. (23) and (24). We also derive the two critical ZF magnitudes given by Eq. (25). Assuming a sinusoidal ZF [Eq. (21)], the drifton Hamiltonian can be written as

H(y,py)=pxu0cosqypx1+px2+py2(β+q2u0cosqy),
(A1)

where px is a constant parameter, and we assume px > 0 without loss of generality. Since driftons move along constant-H surfaces, the drifton trajectory can be determined by equating H(y,py) to a constant E, where EH(y0,py0) is determined by the initial location in the phase space. Then, we obtain py as a function of y

py2(y)=(1+px2)EH0(y)H(y)E,
(A2)

where we introduced

H0(y)H(y,py=0)=pxu0cosqypx1+px2(β+q2u0cosqy)
(A3)

and

H(y)H(y,py=)=pxu0cosqy.
(A4)

It is straightforward to show that maxH0=H0(y=0) and minH0=H0(y=±π/q) (and the same for H); one can also quickly show that the range of E is given by

EH(y0,py0)[min(minH,minH0),maxH],
(A5)

where E=maxH is achieved at (y0=0,py0=),E=minH is achieved at (y0=±π/q,py0=), and E=minH0 is achieved at (y0=±π/q,py0=0).

The drifton trajectories can be conveniently studied on the (y,H) plane using the following method. We plot two curves, namely, C0:H=H0(y) and C:H=H(y) (Fig. 6) and draw a horizontal line L that represents H=E [note that E should be within the range (A5)]. If L intersects neither C0 nor C, then py(y) is always finite, and hence, the trajectory is passing. If L intersects C0 at some location, then py(y) = 0. But dpy/dt is nonzero according to Eq. (24); hence, the drifton will bounce back at that location, which indicates a trapped trajectory (provided that L does not intersect C). If L intersects C at some location, then |py(y)|=, which indicates that the drifton is running away in py space while approaching a particular spatial location. (However, the drifton will never reach such locations, because dpy/dt is finite everywhere.) This indicates a runaway trajectory.

FIG. 6.

H=H0(y) (solid curves), H=H(y) (dotted-dashed curves), and H=E (horizontal lines) for (a) u0 = 0.1 (regime 1), (b) u0 = 2 (regime 2), and (c) u0 = 10 (regime 3). In all figures, the parameters are β = 1, q =0.5, and px = 0.5. The definitions of H0(y) and H(y) are given by Eqs. (A3) and (A4). Each horizontal line represents a drifton trajectory. Note that drifton trajectories are confined within the spatial regions indicated by the solid-line parts of these horizontal lines.

FIG. 6.

H=H0(y) (solid curves), H=H(y) (dotted-dashed curves), and H=E (horizontal lines) for (a) u0 = 0.1 (regime 1), (b) u0 = 2 (regime 2), and (c) u0 = 10 (regime 3). In all figures, the parameters are β = 1, q =0.5, and px = 0.5. The definitions of H0(y) and H(y) are given by Eqs. (A3) and (A4). Each horizontal line represents a drifton trajectory. Note that drifton trajectories are confined within the spatial regions indicated by the solid-line parts of these horizontal lines.

Close modal

For illustration purpose, let us plot C0 and C in Fig. 6 for different values of u0 and study the corresponding trajectories. Figure 6(a) is for u0 = 0.1, which corresponds to regime 1; L1, L2, and L3 drawn there correspond to runaway, trapped, and passing trajectories. Figure 6(b) is for u0 = 2, which corresponds to regime 2; L1 and L2 drawn there correspond to runaway and trapped trajectories, while no L can simultaneously avoid intersecting both C0 and C, and hence no passing trajectory exists. Figure 6(c) is for u0 = 10, which corresponds to regime 3; in this case, C0 and C intersect at |y|=y* [Eq. (29)], and hence every L intersects C, giving a runaway trajectory. Note that the horizontal lines in Fig. 6 are divided into solid-line parts and dashed-line parts; the actual drifton trajectories are confined in the spatial regions indicated by the solid-line parts.

More formally, this method of classifying trajectories can also be presented as follows. For a passing trajectory, L intersects neither C0 nor C; in other words, if the criterion

P:E>maxH0andE<minH
(A6)

is satisfied, then we have a passing trajectory. For a trapped trajectory, L intersects C0 but not C, and hence the corresponding criterion is

T:minH0EmaxH0andE<minH.
(A7)

For a runaway trajectory, L intersects C, and hence the corresponding criterion is

R:minHEmaxH.
(A8)

Consequently, the conditions for each type of trajectory to exist are given as follows. First, passing trajectories exist when maxH0<minH (otherwise, the criterion P is never satisfied for any E); this gives

u0<β2(1+px2)q2uc,1β2q2,
(A9)

where the equality sign applies at px = 0. Next, trapped trajectories exist when minH0<minH (otherwise, the criterion T is never satisfied for any E); this gives

u0<uc,2βq2.
(A10)

Note that under the GO assumption q2 ≪ 1, we have uc,1uc,2. Finally, runaway trajectories always exist, since one can always find E within the range (A5) that satisfy the criterion R.

The above criteria help us to quickly identify the three regimes as follows. If u0 < uc,1, trapped and runaway trajectories exist; passing trajectories also exist because u0 satisfies Eq. (A9) for small enough px; hence, we have regime 1. If uc,1u0 < uc,2, only trapped and runaway trajectories exist; hence, we have regime 2. If u0uc,2, only runaway trajectories exist; hence, we have regime 3.

We also emphasize that runaway trajectories are possible even in the tWKE [Eq. (12)], which is equivalent to setting q2 to zero. In this case, the two critical ZF magnitudes become uc,1 = β/2 and uc,2 = +. Then, u0 < uc,2 is satisfied automatically, and the system is always either in regime 1 or in regime 2. However, the criterion R can still be satisfied even after setting q2 to zero, so runaway trajectories are still possible. Moreover, following the same argument as above, we find that, when u0 ≥ uc,1, passing trajectories also disappear entirely.

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
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.
C.
Connaughton
,
S.
Nazarenko
, and
B.
Quinn
,
Phys. Rep.
604
,
1
(
2015
).
8.
R.
Trines
,
R.
Bingham
,
L. O.
Silva
,
J. T.
Mendonça
,
P. K.
Shukla
, and
W. B.
Mori
,
Phys. Rev. Lett.
94
,
165002
(
2005
).
9.
M. A.
Malkov
,
P. H.
Diamond
, and
A.
Smolyakov
,
Phys. Plasmas
8
,
1553
(
2001
).
10.
E.-J.
Kim
and
P. H.
Diamond
,
Phys. Plasmas
9
,
4530
(
2002
).
11.
A. A.
Vedenov
,
A. V.
Gordeev
, and
L. I.
Rudakov
,
Plasma Phys.
9
,
719
(
1967
).
12.
P.
Kaw
,
R.
Singh
, and
P. H.
Diamond
,
Plasma Phys. Controlled Fusion
44
,
51
(
2002
).
13.
R.
Singh
,
R.
Singh
,
P.
Kaw
,
Ö. D.
Gürcan
, and
P. H.
Diamond
,
Phys. Plasmas
21
,
102306
(
2014
).
14.
M.
Sasaki
,
T.
Kobayashi
,
K.
Itoh
,
N.
Kasuya
,
Y.
Kosuga
,
A.
Fujisawa
, and
S.-I.
Itoh
,
Phys. Plasmas
25
,
012316
(
2018
).
15.
A. I.
Smolyakov
and
P. H.
Diamond
,
Phys. Plasmas
6
,
4410
(
1999
).
16.
A. I.
Smolyakov
,
P. H.
Diamond
, and
V. I.
Shevchenko
,
Phys. Plasmas
7
,
1349
(
2000
).
17.
A. I.
Smolyakov
,
P. H.
Diamond
, and
M.
Malkov
,
Phys. Rev. Lett.
84
,
491
(
2000
).
18.
K.
Srinivasan
and
W. R.
Young
,
J. Atmos. Sci.
69
,
1633
(
2012
).
19.
J. B.
Parker
and
J. A.
Krommes
,
Phys. Plasmas
20
,
100703
(
2013
).
20.
J. B.
Parker
and
J. A.
Krommes
,
New J. Phys.
16
,
035006
(
2014
).
21.
J. B.
Parker
,
J. Plasma Phys.
82
,
595820602
(
2016
).
22.
D. E.
Ruiz
,
J. B.
Parker
,
E. L.
Shi
, and
I. Y.
Dodin
,
Phys. Plasmas
23
,
122304
(
2016
).
23.
J. A.
Krommes
and
C.-B.
Kim
,
Phys. Rev. E
62
,
8508
(
2000
).
24.
J. B.
Parker
,
Phys. Plasmas
25
,
055708
(
2018
).
25.
H.
Zhu
,
Y.
Zhou
,
D. E.
Ruiz
, and
I. Y.
Dodin
,
Phys. Rev. E
97
,
053210
(
2018
).
26.
I. B.
Bernstein
,
J. M.
Greene
, and
M. D.
Kruskal
,
Phys. Rev.
108
,
546
(
1957
).
28.
H.
Zhu
,
Y.
Zhou
, and
I. Y.
Dodin
, e-print arXiv:1805.02233.
29.
A.
Hasegawa
and
K.
Mima
,
Phys. Rev. Lett.
39
,
205
(
1977
).
30.
W. D.
Dorland
, Ph.D. thesis,
Princeton University
,
1993
.
31.
G. W.
Hammett
,
M. A.
Beer
,
W.
Dorland
,
S. C.
Cowley
, and
S. A.
Smith
,
Plasma Phys. Controlled Fusion
35
,
973
(
1993
).
32.
N. D.
Cartwright
,
Physica A
83
,
210
(
1976
).
33.
J. R.
Herring
,
J. Atmos. Sci.
20
,
325
(
1963
).
34.
D. E.
Ruiz
,
M. E.
Glinsky
, and
I. Y.
Dodin
, e-print arXiv:1803.10817.
35.
R. D.
Wordsworth
,
Phys. Fluids
21
,
056602
(
2009
).
36.
A recently study [
J. C.
Li
and
P. H.
Diamond
,
Phys. Plasmas
25
,
042113
(
2018
)] indicates “resonant vorticity mixing” as another mechanism for ZF saturation. However, the effect arises from the non-adiabatic part of the electron response to DWs, which is absent in our model.
37.
H.
Zhu
,
Y.
Zhou
, and
I. Y.
Dodin
(
2018
). “On the structure of the drifton phase space and its relation to the Rayleigh--Kuo criterion of the zonal-flow stability,”
Zenodo
. .
38.
F.
Rath
,
A. G.
Peeters
,
R.
Buchholz
,
S. R.
Grosshauser
,
F.
Seiferling
, and
A.
Weikl
,
Phys. Plasmas
25
,
052102
(
2018
).
39.
B. N.
Rogers
,
W.
Dorland
, and
M.
Kotschenreuther
,
Phys. Rev. Lett.
85
,
5336
(
2000
).
40.
B. N.
Rogers
and
W.
Dorland
,
Phys. Plasmas
12
,
062511
(
2005
).
41.
R.
Numata
,
R.
Ball
, and
R. L.
Dewar
,
Phys. Plasmas
14
,
102312
(
2007
).
42.
D. A.
St-Onge
,
J. Plasma Phys.
83
,
905830504
(
2017
).
43.
R.
Singh
,
H.
Jhang
, and
H. K.
Kaang
,
Phys. Plasmas
23
,
074505
(
2016
).
44.
J. B.
Marston
,
W.
Qi
, and
S. M.
Tobias
, e-print arXiv:1412.0381.