On the structure of the drifton phase space and its relation to the Rayleigh–Kuo criterion of the zonal-flow stability

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).²³ 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 predicts^12–14 nonlinear structures à la Bernstein–Greene–Kruskal (BGK) waves,²⁶ 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.²⁷ 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.²⁹ Consider a collisionless plasma in a uniform magnetic field $B0$ in the z direction, with the equilibrium gradient of the background electron density n₀ in the y direction (Fig. 1). Ions are assumed cold, while electrons are assumed to have a finite temperature T_e. Suppose that perturbations to the electric field E are electrostatic, $E=−∇δφ$⁠, where $δφ(t,x)$ is the corresponding electrostatic potential on the two-dimensional plane $x≐(x,y)$⁠. The electron response to E is adiabatic (yet see below), while the ion response can be described by the $E×B0$ drift and the polarization drift. Then, assuming the quasi-neutrality condition, the evolution of δφ is described by

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

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

where $β≐V*/cs$ is treated as a (positive) constant. Equation (2) represents the original Hasegawa–Mima equation.

Let us introduce the zonal average as $⟨f⟩≐∫0Lxfdx/Lx$⁠, where L_x is the system length in the x direction. Then, perturbations governed by Eq. (1) include ZFs and 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:

where $â$ is an operator such that $â=1$ for DWs and $â=0$ for ZFs.^30,31 Equations (3) and (4) constitute the so-called gHME,²³ 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

which in the GO limit can be understood as the drifton phase-space distribution.²² (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

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 work³⁴ has also gone beyond the quasilinear approximation.) Then, the evolution equations for W and U are^21,22

where $pD2≐1+px2+py2$ and {⋅, ⋅} is the canonical Poisson bracket, namely,

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

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,

where the subscript “t” stands for “traditional.” The iWKE conserves the total enstrophy (per unit length in x) Z_total = Z_DW + Z_ZF and the total energy (per unit length in x) E_total = E_DW + E_ZF of the ZF–DW system, where

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.²² Due to the presence of nonzero Γ [Eq. (11)], Z_DW 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) becomes^21,35

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

Equations (19) and (20) describe the drifton dynamics within the iWKE. Since $H$ does not depend on x, p_x is also conserved along the trajectory. Therefore, the drifton dynamics can be studied on the (y, p_y) plane with p_x 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,

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

and Eqs. (19) and (20) become

Due to the assumed GO approximation, we limit our consideration to the regime where q² ≪ 1 (in dimensional form, $q2≪ρs−2$⁠). We also assume that p_x 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 u₀ compares with the two critical values (the derivations are given in  Appendix)

The GO approximation implies 0 < u_c,1 ≪ u_c,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.

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

and the unstable stationary point (saddle) at

(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 tWKE^12–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,²⁶ 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 $ẏ∝py−3→0$⁠, 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 u₀. 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 u₀. In the second regime, when u_c,1 ≤ u₀ ≤ u_c,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 q² is effectively set to zero. This leads to u_c,1 = β/2, which remains finite. Therefore, passing trajectories also disappear entirely if u₀ ≥ u_c,1. On the other hand, setting q² to zero leads to u_c,2 = +∞; hence, the regime 3 below is impossible in the tWKE.

Regime 3—The third regime corresponds to u₀ ≥ u_c,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

Here, ±y_* are the locations where $U″=β$⁠; namely,

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,²⁷ 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

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 iWKE^21,25

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.²⁰ We launch the simulation with an initial Gaussian DW distribution

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

with small U_q. It is found that depending on the strength of the DW amplitude W₀ [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.

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,

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

where

and

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

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

[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 u₀ > u_c,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

A similar argument leads to

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, E_DW → 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 p_x ≪ 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 q² > 1, while a GO model relies on the assumption that q² ≪ 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

where W₁ 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 Γ₀ are stationary, Eq. (43) is equivalent to

(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 W₁ 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

where p_x is a constant parameter, and we assume p_x > 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 $E≐H(y0,py0)$ is determined by the initial location in the phase space. Then, we obtain p_y as a function of y

where we introduced

and

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

where $E=max H∞$ is achieved at $(y0=0,py0=∞), E=min H∞$ is achieved at $(y0=±π/q,py0=∞)$⁠, and $E=min H0$ 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 C⁰ nor C^∞, then p_y(y) is always finite, and hence, the trajectory is passing. If L intersects C⁰ at some location, then p_y(y) = 0. But dp_y/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 p_y space while approaching a particular spatial location. (However, the drifton will never reach such locations, because dp_y/dt is finite everywhere.) This indicates a runaway trajectory.

For illustration purpose, let us plot C⁰ and C^∞ in Fig. 6 for different values of u₀ and study the corresponding trajectories. Figure 6(a) is for u₀ = 0.1, which corresponds to regime 1; L₁, L₂, and L₃ drawn there correspond to runaway, trapped, and passing trajectories. Figure 6(b) is for u₀ = 2, which corresponds to regime 2; L₁ and L₂ drawn there correspond to runaway and trapped trajectories, while no L can simultaneously avoid intersecting both C⁰ and C^∞, and hence no passing trajectory exists. Figure 6(c) is for u₀ = 10, which corresponds to regime 3; in this case, C⁰ 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 C⁰ nor C^∞; in other words, if the criterion

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

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

Consequently, the conditions for each type of trajectory to exist are given as follows. First, passing trajectories exist when $max H0<min H∞$ (otherwise, the criterion P is never satisfied for any $E$⁠); this gives

where the equality sign applies at p_x = 0. Next, trapped trajectories exist when $min H0<min H∞$ (otherwise, the criterion T is never satisfied for any $E$⁠); this gives

Note that under the GO assumption q² ≪ 1, we have u_c,1 ≪ u_c,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 u₀ < u_c,1, trapped and runaway trajectories exist; passing trajectories also exist because u₀ satisfies Eq. (A9) for small enough p_x; hence, we have regime 1. If u_c,1 ≤ u₀ < u_c,2, only trapped and runaway trajectories exist; hence, we have regime 2. If u₀ ≥ u_c,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 q² to zero. In this case, the two critical ZF magnitudes become u_c,1 = β/2 and u_c,2 = +∞. Then, u₀ < u_c,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 q² to zero, so runaway trajectories are still possible. Moreover, following the same argument as above, we find that, when u₀ ≥ u_c,1, passing trajectories also disappear entirely.