The wave kinetic equation (WKE) describing drift-wave (DW) turbulence is widely used in the studies of zonal flows (ZFs) emerging from DW turbulence. However, this formulation neglects the exchange of enstrophy between DWs and ZFs and also ignores effects beyond the geometrical-optics limit. We derive a modified theory that takes both of these effects into account, while still treating DW quanta (“driftons”) as particles in phase space. The drifton dynamics is described by an equation of the Wigner–Moyal type, which is commonly known in the phase-space formulation of quantum mechanics. In the geometrical-optics limit, this formulation features additional terms missing in the traditional WKE that ensure exact conservation of the total enstrophy of the system, in addition to the total energy, which is the only conserved invariant in previous theories based on the WKE. Numerical simulations are presented to illustrate the importance of these additional terms. The proposed formulation can be considered as a phase-space representation of the second-order cumulant expansion, or CE2.
The formation of zonal flows (ZFs) is a problem of fundamental interest in many contexts, including the physics of planetary atmospheres, astrophysics, and fusion science.1–7 In particular, the interaction of ZFs and drift-wave (DW) turbulence in laboratory plasmas significantly affects the transport of energy, momentum, and particles, so understanding it is critical to improving the plasma confinement. But modeling the underlying physics remains a difficult problem. The workhorse approach to describing the DW-ZF coupling is the wave kinetic equation (WKE),5,8 but it is limited to the ray approximation9 and, in fact, is oversimplified even as a geometrical-optics10 (GO) model. That leads to missing essential physics, as was recently pointed out in Ref. 11 and will be elaborated below. These issues can be fixed by using the more accurate quasilinear approach known as the second-order cumulant expansion, or CE2,12–16 whose applications to DW-ZF physics were pursued in Refs. 17–20. However, the CE2 is less intuitive than the WKE, and its robustness with respect to further approximations remains obscure. Having an approach as accurate as the CE2 and as intuitive as the WKE would be more advantageous.
Here, we propose such an approach for a DW turbulence model based on the generalized Hasegawa–Mima equation (gHME).21,22 The idea is as follows. We start by splitting the gHME into two coupled equations that describe ZFs and fluctuations, respectively, and then linearize the equation for fluctuations, like in the CE2 approach. We notice then that this linearized equation is similar to that describing a quantum particle that is governed by a generalized (non-Hermitian) Hamiltonian. By drawing on this analogy, we then formulate an exact (modulo quasilinear approximation) kinetic equation for such particle, which is akin to the so-called Wigner–Moyal equation in quantum mechanics.23–25,42,43
Compared to the CE2, the Wigner–Moyal formulation is arguably more intuitive, namely, for two reasons: (i) like the traditional WKE (hereafter denoted by tWKE), it permits viewing DW quanta (“driftons”) as particles, except that now driftons are quantumlike particles, i.e., have nonzero wavelengths; and (ii) the separation between Hamiltonian effects and dissipation remains transparent and unambiguous even beyond the GO approximation. Compared to the tWKE, the new approach is also more accurate because (i) it captures effects beyond the GO limit, and (ii) even in the GO limit, it predicts corrections to the tWKE that emerge from the newly found corrections to the drifton dispersion. (In this aspect, our paper can be understood as an expansion of the GO approximation introduced in Ref. 11.) These corrections are essential as they allow DW-ZF enstrophy exchange, which is not included in the tWKE. By deriving the GO limit from first principles, we eliminate this discrepancy and obtain a formulation that exactly conserves the total enstrophy (as opposed to the DW enstrophy conservation predicted by the tWKE) and the total energy, in precise agreement with the underlying gHME. We also illustrate the substantial difference between the GO limit of our formulation and the tWKE using numerical simulations.
The paper is organized as follows. In Sec. II, we introduce the gHME and its quasilinear approximation. In Sec. III, we derive the Wigner–Moyal formulation. In Sec. IV, we rederive the dispersion relation for the linear growth rate of ZFs. In Sec. V, we derive a corrected WKE that, in contrast to the tWKE, conserves both the total enstrophy and energy. Numerical simulations are presented to compare the new WKE with the tWKE. In Sec. VI, we summarize our main results. Auxiliary calculations are presented in Appendixes. This includes a brief introduction to the Weyl calculus that we extensively use in our paper ( Appendix A), a spectral representation of our formulation ( Appendix B), and proofs of the conservation properties of our models ( Appendix C).
II. BASIC MODEL
Our formulation is based on the gHME21,22
which is widely used to describe the electrostatic two-dimensional (2-D) turbulent flows both in a magnetized plasma with a density gradient and in an atmospheric fluid on a rotating planet, where the role of DWs is played by Rossby waves.1,20 Both contexts will be described on the same footing, so our results are applicable to DWs and Rossby waves equally. We assume the usual geophysical coordinate system, where x = (x, y) is a 2-D coordinate, the x-axis is the ZF direction, and the y-axis is the direction of the local gradient of the plasma density or of the Coriolis parameter. (In the context of fusion plasmas, a different choice of coordinates is usually preferred in literature, where x and y are swapped.) The constant β is a measure of this gradient. The function ψ(x, t) is the electric potential or the stream function, is the fluid velocity on the x plane, and ez is a unit vector normal to this plane. The function w(x, t) is the generalized vorticity given by , where is an operator such that in parts of the spectrum corresponding to DWs and in those corresponding to ZFs. (The symbol denotes definitions.) Also, LD is the plasma sound radius or the deformation radius. (For plasmas, one can take LD = 1 in normalized units.21 Also, the barotropic model used in geophysics is recovered in the limit LD → ∞.12–15) The term Q(x, t) describes the external forces and dissipation. Systems with Q = 0 will be called isolated.
Let us decompose the fields into their zonal-averaged and fluctuating components, denoted with bars and tildes, respectively. (For any g, its zonal average is , where Lx, henceforth assumed equal to one, is the system length along x axis.) In particular, , where the two components of the generalized vorticity are18
where . Equations for and are obtained by taking the zonal-average and fluctuating parts of Eq. (1). This gives
It is convenient to rewrite Eqs. (4) in terms of the ZF velocity , whose only component U(y, t) is . Specifically, one has , and . We will also assume and . Here, is some external force with zero zonal average (eventually, we will assume it to be a white noise), and the constant coefficients μdw and μzf are intended to emulate the dissipation of DWs and ZFs caused by the external environment. Then, Eqs. (4) become
Equations (6) are the same model as the one that underlies the CE2. Although not exact, this model is useful because it captures the key aspects of ZF dynamics, such as formation and merging of zonal jets.15,20,26 Below, we use it to derive a formulation of DW-ZF interactions alternative to the CE2.
III. WIGNER–MOYAL FORMULATION
A. State vector
Consider a family of all reversible linear transformations of of the form . These transformations map into some family of image functions. Since these functions are mutually equivalent up to an isomorphism, the resulting family can be viewed as a single object, a time-dependent “state vector” . (Analogous definitions will be assumed also for and .) The original function is then understood as a projection of , namely, as its “coordinate representation” given by . Here, are the eigenstates of the position operator normalized such that . This definition of a field is similar to that used in quantum mechanics for describing probability amplitudes.27 Hence, it is convenient to describe the dynamics of using a quantumlike formalism. This is done as follows.
In addition to the coordinate operator , we introduce a momentum (wave-vector) operator such that, in the coordinate representation, . Accordingly, , where
Hence, Eq. (6a) can be represented in the following form:
The operator is given by
Also, , and the prime above U henceforth denotes ∂y; in particular, .
B. Generalized von Neumann equation
Let us express Eq. (9) as , where and are the Hermitian and anti-Hermitian parts of , correspondingly. Explicitly, these operators can be written as
where . In particular, taking the trace of this equation also gives an equation for the “total number of DW quanta,” ; namely,
This indicates that determines the loss of quanta, or dissipation of DWs. [In particular, the term in Eq. (10b) is responsible for DW dissipation to the external environment, whereas the term destroys DW quanta while conserving the total enstrophy, as will be discussed in Sec. III E.] Also, determines the conservative dynamics of DWs and thus can be understood as the drifton Hamiltonian. (The non-Hermitian operator will be attributed as the generalized Hamiltonian.) Notice that the distinction between dissipation and Hamiltonian effects remains unambiguous even beyond the GO approximation.
Equation (11) can be understood as a generalized von Neumann equation akin to the one that commonly emerges in quantum mechanics. A standard approach to such equation is to project it on the phase space using the Weyl transform. Hence, we proceed as follows. (Readers who are not familiar with the Weyl calculus are encouraged to read Appendix A before continuing further.)
C. Wigner–Moyal equation
Let us introduce W as the Weyl symbol of , i.e.,
which is real because is Hermitian. In quantum mechanics, a similar construct is known as the Wigner function,28 so one can readily identify the physical meaning of W. Specifically, in the regime, when the ray approximation applies and dissipation is negligible, W/(2π)2 represents the phase-space probability density of driftons [the numerical coefficient (2π)2 comes from Eq. (A4)], while beyond the GO limit, it can be considered as a generalization of this probability density.29 Using the fact that is real, one can also cast W as
which also implies
One can interpret the right-hand side of Eq. (14) as the local spatial spectrum of the correlation function of w. Hence, W will be called the DW spectral function.
Here and are Moyal's “sine bracket” [Eq. (A10)] and “cosine bracket” [Eq. (A12)]. The functions HH(y, p, t), HA(y, p, t), and F(x, p, t) are the Weyl symbols of , and , respectively. In particular, using Eq. (A5) and the fact that U is independent of x, one obtains
where . By analogy with quantum mechanics, we call Eq. (16) a Wigner–Moyal equation.
Next, let us consider the zonal average of this equation
where . We adopt the ergodic assumption, namely, that the zonal average is equivalent to the ensemble average [denoted ] over realizations of the random force (e.g., as done in Ref. 18). To calculate , consider integrating Eq. (8) on a time interval (t0, t). The result can be written as where the indexes denote the times at which functions are evaluated, and . We assume
where Ξ is a correlation function that is homogeneous in x but not necessarily in y.15,20 Since can be affected by only if , one has . Hence,
where ‘h.c.’ denotes the “Hermitian conjugate.” In other words, once the correlation function Ξ of is specified, can be readily calculated as the Fourier image of Ξ.
This concludes the calculation of the functions that determine the evolution of through Eq. (19). However, these functions generally depend on U, so an additional equation for U is needed to make the theory self-consistent. This equation is derived as follows.
D. Equation for the zonal-flow velocity
Returning to Eq. (6b), we rewrite the nonlinear term as
E. Main equations and conservation laws
Let us slightly change the notation and summarize the above equations in the following form:
The function can be understood as the Weyl symbol of the drifton Hamiltonian, whereas Γ determines the dissipation of DW quanta that is caused specifically by DW interaction with ZFs. This is explained as follows. Since Eqs. (25) are exact within the quasilinear approximation (modulo the ergodic assumption), they inherit the same conservation laws as the original quasilinear model given by Eqs. (6). Specifically, for isolated systems ( and μdw, zf = 0), Eqs. (25) and (26) exactly conserve the total enstrophy and energy [Eq. (5)]
rather than their DW and ZF components. (A direct proof is given in Appendix C 1.) For completeness, we present expressions for these components:
The conservative equations (25) and (26), which we attribute as the Wigner–Moyal formulation of DW-ZF interactions, constitute the main result of our work. This formulation can be understood as an alternative phase-space representation of the CE2 since it is derived from the same quasilinear model. However, the Wigner–Moyal formulation is arguably more intuitive than the CE2, namely, for two reasons: (i) Like in the tWKE, driftons are treated as particles, except now they are quantumlike particles, i.e., have nonzero wavelengths; hence, one is not constrained to the GO limit. (ii) Also, the separation between Hamiltonian effects and dissipation remains transparent and unambiguous even beyond the GO approximation. The Wigner–Moyal formulation elucidates the link between the WKE formalism and the CE2 and also helps make approximations rigorous by making them systematic. Below, these and other applications are discussed in further detail.
IV. GROWTH RATE OF ZONAL FLOWS
A. Basic equations
To demonstrate the convenience of the Wigner-Moyal formulation, let us apply it to rederive the rate of the linear zonostrophic instability, i.e., the growth rate of weak ZFs. Suppose a homogeneous equilibrium with zero ZF velocity and some DW spectral function . [As pointed out in Sec. III C, the corresponding represents the phase-space probability distribution of driftons.] Consider small perturbations to this equilibrium; namely,
Here, the constant q serves as the modulation wave number, and the constant γ is the instability rate to be found. The linearization of Eq. (25a) leads to
where we assume the notation for any A. Solving for in terms of Uq leads to
Then, Eq. (25b) yields
Due to Eq. (A16), this can be simplified as follows:
After substituting the expression for , one gets
B. Zonal flows with nonzero group velocity
As a side note, it is commonly thought that ZFs only grow in situ; i.e., Re γ > 0 with Im γ = 0. There have been questions over whether it is possible to have unstable zonal modes at nonzero Im γ, which implies nonzero group velocity.32 Here we show, by presenting an example, that the answer is yes. Let us consider the following steady state:
After integrating, Eq. (32) can be cast as follows:
where is the DW group velocity in the absence of ZFs, , and . Numerical solutions of Eq. (34) are presented in Fig. 1. As shown, solutions for γ are complex over some interval of ky. This example shows the existence of “traveling” unstable ZF modes.
Additional insights on this phenomenon can be obtained by considering the limit, where and q ≪ ky. In this case, Eq. (34) simplifies to
One may consider this as the GO limit of Eq. (34). Equation (35) predicts four nontrivial solutions for X, which are given by . Different regimes for the solutions can be deduced. When σ ≥ 1, the solutions γ are purely real. For the parameters in Fig. 1, this regime corresponds to . In the interval 0 < σ < 1 corresponding to , γ is complex-valued. In the interval −1/8 ≤ σ ≤ 0, which corresponds to , the solutions are purely imaginary. Finally, in the interval σ < −1/8 corresponding to , two solutions γ are purely imaginary, and two other solutions are purely real. The different regimes identified by solving Eq. (35) are consistent with the observed numerical solutions of the exact dispersion relation (34). In Sec. V, we will explore the GO limit of the DW-ZF interactions in more detail.
V. GEOMETRICAL-OPTICS LIMIT AND THE WAVE KINETIC EQUATION
Let us assume that the characteristic wavelengths for ZFs and DWs are λzf and λdw, respectively, and
Hence, the following estimates will be adopted:
where H denotes both and Γ. (The latter estimate is given for the maximum of , which is realized at .37) This gives
Notice that our WKE differs from the tWKE, which assumes a simpler dispersion of DWs; namely,
The numerical results presented in Figs. 2–4 illustrate the importance of the difference between our WKE and the tWKE [subfigures (a) and (b), respectively]. As seen in Fig. 2, while our WKE model predicts ZFs with a particular λzf, the scale of ZFs predicted by tWKE is determined by nothing but the grid size that is used in simulations. This is because the tWKE predicts that the rate of the zonostrophic instability γ (Sec. IV) scales linearly with the ZF wave number q, so ZFs are produced at the largest q that is allowed (cf. Ref. 11).
Consider also the enstrophy plots in Fig. 3. To aid our discussion, we added the plots of the enstrophy that the external forcing injects into the DW-ZF system; namely, . Within our model, the total enstrophy remains always smaller than , which is natural, since the simulation is done for μdw,zf > 0. In contrast, the tWKE model predicts that can surpass , which is unphysical. In addition, the values of the ZF and total enstrophies predicted by the tWKE are several times larger than those predicted by our model.
For the sake of completeness, Fig. 4 also presents the corresponding energies and the energy introduced by the external force, . In both cases, , which is in agreement with the fact that both models conserve the total energy of an isolated system. Still, the tWKE predicts very different results quantitatively even though the tWKE model [Eqs. (43)] is seemingly close to ours [Eqs. (42)].
The goal of this paper is to propose a new formulation of DW-ZF interactions that is more accurate than the tWKE and, simultaneously, more intuitive than the CE2. We adopt the same model [Eq. (6)] that was previously applied to derive the CE2. Then, we manipulate it using the Weyl calculus to produce a phase-space formulation of DW-ZF interactions. The resulting formulation [Eqs. (25) and (26)] is akin to a quantum kinetic theory and involves a pseudodifferential Wigner–Moyal equation. To facilitate its numerical implementation in the future, we also present an integral representation of our main equations ( Appendix B).
On one hand, this Wigner–Moyal formulation can be understood as an alternative representation to the CE2 since both models use the same assumptions. For example, we show that it leads to the same linear growth rate of weak ZFs as that obtained from the CE2 (Sec. IV). On the other hand, the Wigner–Moyal formulation is arguably more intuitive than the CE2 for two reasons: (i) it permits treating driftons as particles (i.e., as objects traveling in phase space), except now they are quantumlike particles with nonzero wavelengths; and (ii) the separation between Hamiltonian effects and dissipation remains unambiguous even beyond the GO limit.
Compared to the tWKE, the new approach is also more precise because (i) it captures effects beyond the GO limit and (ii) even in the GO limit, it predicts corrections to the tWKE that emerge from the newly found corrections to the drifton dispersion (Sec. V). These corrections are essential as they allow DW-ZF enstrophy exchange, which is not included in the tWKE. By deriving the GO limit from first principles, we eliminate this discrepancy and arrive at a model that exactly conserves the total enstrophy (as opposed to the DW enstrophy conservation predicted by the tWKE) and the total energy, in agreement with the underlying gHME. We also illustrate the substantial difference between the GO limit of our WKE model and the tWKE using numerical simulations.
This work can be expanded at least in two directions. First, the difference between the Wigner–Moyal formulation and the newly proposed WKE can be assessed quantitatively using numerical simulations. Second, the analytic methods we proposed here can be extended to other turbulence models, such as those in Refs. 39 and 40. The anticipated benefit is that more accurate equations would be derived that would respect the fundamental conservation laws that existing theories may be missing otherwise.
The authors thank J. A. Krommes for valuable discussions. This work was supported by the U.S. DOE through Contract Nos. DE-AC02-09CH11466 and DE-AC52-07NA27344, by the NNSA SSAA Program through DOE Research Grant No. DE-NA0002948, and by the U.S. DOD NDSEG Fellowship through Contract No. 32-CFR-168 a.
APPENDIX A: WEYL CALCULUS
We shall refer to this description of the operators as a phase-space representation since Weyl symbols are functions of the 2 n-dimensional ray phase space (x, p). Conversely, the inverse Weyl transformation is
In particular, notice that, for any operator , its matrix elements in the coordinate representation, , can be expressed as
so A(x, p) can be understood as a spectrum of . In particular,
Other properties of the Weyl transform that we use in this paper are as follows:
As seen from Eq. (A3), for any operator , its trace can be expressed as
If A(x, p) is the Weyl symbol of , then A*(x, p) is the Weyl symbol of . As a corollary, the Weyl symbol of a Hermitian operator is real.
- For any , the corresponding Weyl symbols satisfy23,24(A5)Here, ⋆ is the Moyal product, which is given by(A6)and is the Janus operator, which is given by(A7)The arrows indicate the directions in which the derivatives act, and is the canonical Poisson bracket; namely,(A8)
The Moyal product is associative; i.e.,
- The anti-symmetrized Moyal product defines the so-called Moyal bracket(A10)Because of the latter equality, this bracket is also called the sine bracket. To lowest order,(A11)
- The symmetrized Moyal product is defined as(A12)Because of the latter equality, this bracket is also called the cosine bracket. To lowest order,(A13)
In a wave equation of the Wigner–Moyal type, such as Eq. (25a), neglecting higher-order phase-space derivatives in the sine and cosine brackets leads to a WKE. Higher-order wave effects, such as diffraction and tunneling, are lost in this limit. For this reason, it is called the GO approximation, or the ray approximation.
- Assuming that fields vanish at infinity rapidly enough, the phase-space integral of the Moyal product of two symbols equals the integral of the regular product of these symbols:(A14)
As a corollary(A15)(A15a)(A15b)
For any constant q, one has
As a corollary, one has
For any constants k and q, one can also show that
- Now we tabulate some Weyl transforms of various operators. (We use a two-sided arrow to show the correspondence with the Weyl transform.) First of all, the Weyl transforms of the identity, position, and momentum operators are given by(A19)If f and g are any two functions, then(A20)Similarly, using Eq. (A6), we have(A21)(A22)
One may also notice the connection between these relations and the commutation relation between operators; i.e., .
APPENDIX B: SPECTRAL REPRESENTATION OF THE WIGNER–MOYAL FORMULATION
To facilitate numerical implementations of the Wigner–Moyal formulation in the future, we propose an integral form of Eqs. (25) using a spectral representation. (Numerical simulations of the CE2 theory were reported in Refs. 12–14.) The assumed notation for the Fourier representation of any A(y, p, t) will be
We start by rewriting Eqs. (A17) as
where . This leads to
so one obtains
Also, using Eq. (A18), we obtain
where for any Ar(p, t). Also
By inserting these expressions into Eq. (25a), we obtain
Also, using that
one gets the following representation of Eq. (25b):
APPENDIX C: CONSERVATION OF THE TOTAL ENSTROPHY AND ENERGY
1. Wigner–Moyal model
To show conservation of energy, we obtain
Here, we used the fact that the Taylor expansion of Eq. (A10) for the Moyal bracket consists of total derivatives on y, so its integral over y is zero. The other terms can be expressed as follows. First of all,
Also, using Eq. (A18), we obtain
2. WKE model
where we used Eq. (42b) and the fact that the integral of the Poisson bracket over all phase space is zero. In contrast, the tWKE does not conserve the total enstrophy because Γ = 0 [see Eqs. (43)]. This is also understood as follows: the tWKE manifestly conserves , whereas is obviously not conserved; thus, cannot be conserved either.
To show conservation of energy, we obtain
where the integral of over y is zero because it can be written as a total derivative on y. Then, we have
Note that an analysis using the tWKE model leads to an expression similar to Eq. (C9) with only the first two terms in the integrand. These terms cancel out, thus showing that the tWKE model also conserves the total energy.
Specifically, this means that the traditional WKE has a unique set of characteristics at each phase-space location. Using the terminology that is adopted in the present paper, the ray approximation can be defined formally as a model based on the approximations (A11) and (A13). This model is applicable under the assumption that ZFs have scales larger compared to the characteristic wavelength of DWs, which is also known as the condition of geometrical optics. See below for details.
The only peculiarity of our problem in this respect is that is real rather than complex. However, that is just a matter of initial conditions. The equation governing still has a quantumlike form.
Beyond the ray approximation, the probability density of driftons cannot be defined, because wave quanta cannot be assigned specific locations in phase space due to the uncertainty principle. Unlike the true probability density, the Weyl symbol of the density matrix is not positive-definite in this case; however, it does remain real.
Decreasing the hyperviscosity parameter by an order of magnitude to ν = 0.0001 leads to qualitatively similar results. Quantitatively, the steady-state ZF enstrophy increases by about 25% when using the tWKE. This can be understood by how U in the tWKE model is dominated by short wavelength components, which are more strongly affected by hyperviscosity.
In the barotropic limit (LD → ∞), the ordering (37) can no longer be satisfied. This means that a GO approximation of Eqs. (25) and (26) is not possible in this limit. This problem originates from the fact that the Hamiltonian is divergent at p = 0. In the unnatural case where for p < p0, then we can estimate . Then, a GO approximation is possible if .
The fundamental reason for this discrepancy is that the tWKE assumes to be the true distribution probability of driftons, but is not. The true probability distribution would satisfy a strictly conservative equation. In the GO limit, it can always be found, at least numerically, as will be discussed elsewhere. The special case when is negligible is discussed in Ref. 11.