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.

## I. INTRODUCTION

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 approximation^{9} and, in fact, is oversimplified even as a geometrical-optics^{10} (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 gHME^{21,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, $v=ez\xd7\u2207\psi $ is the fluid velocity on the **x** plane, and **e**_{z} is a unit vector normal to this plane. The function *w*(**x**, *t*) is the generalized vorticity given by $w\u2009\u2250(\u22072\u2212LD\u22122\alpha \u0302)\psi $, where $\alpha \u0302$ is an operator such that $\alpha \u0302=1$ in parts of the spectrum corresponding to DWs and $\alpha \u0302=0$ in those corresponding to ZFs. (The symbol $\u2250$ denotes definitions.) Also, *L*_{D} is the plasma sound radius or the deformation radius. (For plasmas, one can take *L*_{D} = 1 in normalized units.^{21} Also, the barotropic model used in geophysics is recovered in the limit *L*_{D} → *∞*.^{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 $g\xaf\u2250\u222bdx\u2009g/Lx$, where *L _{x}*, henceforth assumed equal to one, is the system length along

*x*axis.) In particular, $w=w\xaf+w\u0303$, where the two components of the generalized vorticity are

^{18}

where $\u2207D2\u2009\u2250\u2009\u22072\u2212LD\u22122$. Equations for $w\u0303$ and $w\xaf$ are obtained by taking the zonal-average and fluctuating parts of Eq. (1). This gives

*f*

_{NL}will be neglected. Hence,

^{12}In isolated systems, both sets of equations conserve the enstrophy $Z$ and the energy $E$ (strictly speaking, free energy), which are defined as

It is convenient to rewrite Eqs. (4) in terms of the ZF velocity $v\xaf=exU$, whose only component *U*(*y*, *t*) is $U=\u2212\u2202y\psi \xaf$. Specifically, one has $v\u0303\xb7\u2207w\xaf=\u2212(\u2202x\psi \u0303)(\u2202y2U),v\xaf\xb7\u2207w\u0303=U\u2202xw\u0303$, and $v\u0303\xb7\u2207w\u0303\xaf=\u2212\u2202y2\u2009v\u0303xv\u0303y\xaf$. We will also assume $Q\u0303=\xi \u0303\u2212\mu dww\u0303$ and $Q\xaf=\u2212\mu zfw\xaf$. Here, $\xi \u0303$ 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 $w\u0303(x,t)$ of the form $\u222bd2x\u2032\u2009K(x,x\u2032,t)w\u0303(x\u2032,t)$. These transformations map $w\u0303(x,t)$ 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” $|w\u0303\u27e9$. (Analogous definitions will be assumed also for $|\psi \u0303\u27e9$ and $|\xi \u0303\u27e9$.) The original function $w\u0303(x,t)$ is then understood as a projection of $|w\u0303\u27e9$, namely, as its “coordinate representation” given by $w\u0303(x,t)=\u27e8x|w\u0303\u27e9$. Here, $|x\u27e9$ are the eigenstates of the position operator $x\u0302$ normalized such that $\u27e8x\u2032|x\u0302|x\u27e9=x\u27e8x\u2032|x\u27e9=x\u2009\delta (x\u2032\u2212x)$. 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 $|w\u0303\u27e9$ using a quantumlike formalism. This is done as follows.

In addition to the coordinate operator $x\u0302$, we introduce a momentum (wave-vector) operator $p\u0302$ such that, in the coordinate representation, $p\u0302\u2250\u2212i\u2207$. Accordingly, $|w\u0303\u27e9=\u2212p\u0302D2|\psi \u0303\u27e9$, where

Hence, Eq. (6a) can be represented in the following form:

The operator $H\u0302$ is given by

Also, $U\u0302\u2250U(y\u0302,t)$, and the prime above *U* henceforth denotes *∂ _{y}*; in particular, $U\u0302\u2033\u2250\u2202y2\u2009U(y\u0302,t)$.

### B. Generalized von Neumann equation

Let us express Eq. (9) as $H\u0302=H\u0302H+iH\u0302A$, where $H\u0302H\u2250(H\u0302+H\u0302\u2020)/2$ and $H\u0302A\u2250(H\u0302\u2212H\u0302\u2020)/(2i)$ are the Hermitian and anti-Hermitian parts of $H\u0302$, correspondingly. Explicitly, these operators can be written as

where $F\u0302\u2250|\xi \u0303\u27e9\u27e8w\u0303|+|w\u0303\u27e9\u27e8\xi \u0303|$. In particular, taking the trace of this equation also gives an equation for the “total number of DW quanta,” $N\u2250Tr\u2009W\u0302=\u222bd2x\u2009\u27e8x|W\u0302|x\u27e9=\u222bd2x\u2009w\u03032=\u27e8w\u0303|w\u0303\u27e9$; namely,

This indicates that $H\u0302A$ determines the loss of quanta, or dissipation of DWs. [In particular, the term $\mu dw$ in Eq. (10b) is responsible for DW dissipation to the external environment, whereas the term $[U\u0302\u2033,p\u0302xp\u0302D\u22122]\u2212/(2i)$ destroys DW quanta while conserving the total enstrophy, as will be discussed in Sec. III E.] Also, $H\u0302H$ determines the conservative dynamics of DWs and thus can be understood as the *drifton Hamiltonian*. (The non-Hermitian operator $H\u0302$ 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 $W\u0302$, i.e.,

which is real because $W\u0302$ 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 $w\u0303(x,t)$ 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*.

By applying the Weyl transform to Eq. (11), one obtains the following pseudo-differential equation:^{30,44}

Here ${{\xb7,\xb7}}$ and $[[\xb7,\xb7]]$ are Moyal's “sine bracket” [Eq. (A10)] and “cosine bracket” [Eq. (A12)]. The functions *H _{H}*(

*y*,

**p**,

*t*),

*H*(

_{A}*y*,

**p**,

*t*), and

*F*(

**x**,

**p**,

*t*) are the Weyl symbols of $H\u0302H,\u2009H\u0302A$, and $F\u0302$, respectively. In particular, using Eq. (A5) and the fact that

*U*is independent of

*x*, one obtains

where $pD2\u2250p2+LD\u22122$. By analogy with quantum mechanics, we call Eq. (16) a Wigner–Moyal equation.

Next, let us consider the zonal average of this equation

where $W\xaf=W\xaf(y,p,t)$. We adopt the ergodic assumption, namely, that the zonal average is equivalent to the ensemble average [denoted $\u27e8\u27e8\u2026\u27e9\u27e9$] over realizations of the random force $\xi \u0303$ (e.g., as done in Ref. 18). To calculate $F\xaf=\u27e8\u27e8F\u27e9\u27e9$, consider integrating Eq. (8) on a time interval (*t*_{0}, *t*). The result can be written as $|w\u0303t\u27e9=|w\u0303t0\u27e9+|\delta w\u0303t\u27e9+\u222bt0tdt\u2032|\xi \u0303t\u2032\u27e9,$ where the indexes denote the times at which functions are evaluated, and $|\delta w\u0303t\u27e9\u2250\u2212i\u222bt0tdt\u2032H\u0302|w\u0303t\u2032\u27e9$. We assume

where Ξ is a correlation function that is homogeneous in *x* but not necessarily in *y*.^{15,20} Since $|\delta w\u0303t\u27e9$ can be affected by $|\xi \u0303t\u2032\u27e9$ only if $t>t\u2032$, one has $\u27e8\u27e8|\xi \u0303t\u27e9\u27e8\delta w\u0303t|\u27e9\u27e9=0$. Hence,

where ‘h.c.’ denotes the “Hermitian conjugate.” In other words, once the correlation function Ξ of $\xi \u0303$ is specified, $F\xaf$ can be readily calculated as the Fourier image of Ξ.

This concludes the calculation of the functions that determine the evolution of $W\xaf$ 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

where we used Eq. (A3) in the last step. After introducing the averaged vorticity density $W\xaf$, Eq. (6b) becomes

### E. Main equations and conservation laws

Let us slightly change the notation and summarize the above equations in the following form:

*U*(

*y*,

*t*) is the ZF velocity. Also, $F\xaf=F\xaf(y,p)$ is determined by the correlation function of the external noise $\xi \u0303$ (Sec. III C). We have also introduced $H\u2250HH$ and $\Gamma \u2250HA+\mu dw$, or, explicitly

The function $H$ 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 ($F\xaf=0$ 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 $W(p)$. [As pointed out in Sec. III C, the corresponding $W(p)/(2\pi )2$ 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 $A\xb1q\u2250A(p\xb1eyq/2)$ for any *A*. Solving for $W\xafq$ in terms of *U _{q}* leads to

Then, Eq. (25b) yields

Due to Eq. (A16), this can be simplified as follows:

After substituting the expression for $W\xafq$, 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 $vgy\u22502\beta kxky/kD4$ is the DW group velocity in the absence of ZFs, $kD,\xb1q2\u2250kx2+(ky\xb1q/2)2+LD\u22122$, and $X\u2250(\gamma +2\mu dw)/(qvgy)$. Numerical solutions of Eq. (34) are presented in Fig. 1. As shown, solutions for *γ* are complex over some interval of *k _{y}*. This example shows the existence of “traveling” unstable ZF modes.

Additional insights on this phenomenon can be obtained by considering the limit, where $\mu dw,zf\u226a|\gamma |$ and *q* ≪ *k _{y}*. In this case, Eq. (34) simplifies to

where

One may consider this as the GO limit of Eq. (34). Equation (35) predicts four nontrivial solutions for *X*, which are given by $X2=\u22121+4\sigma \xb14[\sigma (\sigma \u22121)]1/2$. 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 $|ky|\u22720.33$. In the interval 0 < *σ* < 1 corresponding to $0.33\u2272|ky|\u22720.82$, *γ* is complex-valued. In the interval −1/8 ≤ *σ* ≤ 0, which corresponds to $0.82\u2272|ky|\u22721.07$, the solutions are purely imaginary. Finally, in the interval *σ* < −1/8 corresponding to $|ky|\u22731.07$, 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 $H$ and Γ. (The latter estimate is given for the *maximum* of $\u2202pH$, which is realized at $p\u223cLD\u22121$.^{37}) This gives

Then, using the lowest-order approximations of the Moyal products ( Appendix A), Eqs. (25) reduce to

*WKE limit*of the Wigner-Moyal formulation. Clearly, $H$ acts as the drifton ray Hamiltonian while Γ acts as the corresponding dissipation rate. [The factors of two in Eq. (41a) are due to the fact that $W\xaf$ is quadratic in the DW amplitude.] In other words, $\omega (y,p,t)\u2250H+i\Gamma \u2212i\mu dw$ can be viewed as the local complex frequency of DWs with a given wave vector

**p**.

Notice that our WKE differs from the tWKE, which assumes a simpler dispersion of DWs; namely,

*U*, these terms remain important for a number of reasons. For example, in the Hamiltonian $H,\u2009U\u2033$ can be comparable to

*β*(as is sometimes the case in geophysics

^{2}). Also, consider the following. In isolated systems, the tWKE is $\u2202tW\xaf={H,W\xaf}$, so it conserves the DW quanta, or, in other words, the DW enstrophy $Zdw$ [Eq. (28a)]. At the same time, the ZF enstrophy $Zzf$ [Eq. (28b)] generally evolves, so the total enstrophy $Z=Zdw+Zzf$ does too. This is in contradiction with the gHME, which conserves $Z$, and can lead to overestimating the ZF velocity and shear generated by DW turbulence.

^{38}In contrast to the tWKE, our formulation is free from such issues because Eqs. (41) and (42) exactly conserve both $Z$ and $E$ (Appendix C 2). Note that, in order to retain this conservation property, it is necessary to keep both $U\u2034$ and $U\u2033$ in Eqs. (42). In this sense, Eqs. (41) and (42) represent the simplest GO model that is physically meaningful in the nonlinear regime. This is in agreement with Ref. 11, where a similar conclusion was made based on comparing the linear zonostrophic instability rate predicted by the CE2. (As a note on terminology, Ref. 11 refers to the tWKE [Eqs. (41) and (43)] as the “Asymptotic WKE,” i.e., the limit obtained when one assumes the ZFs are asymptotically large scale. Also, Ref. 11 refers to Eqs. (41) and (42) as “CE2-GO.”)

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 $Zext$ that the external forcing $F\xaf$ injects into the DW-ZF system; namely, $Zext=(t/2)(2\pi )\u22122\u222bdy\u2009d2p\u2009F\xaf$. Within our model, the total enstrophy $Z$ remains always smaller than $Zext$, which is natural, since the simulation is done for *μ*_{dw,zf} > 0. In contrast, the tWKE model predicts that $Z$ can surpass $Zext$, 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 $Eext$ introduced by the external force, $Eext=(t/2)(2\pi )\u22122\u222bdy\u2009d2p\u2009F\xaf/pD2$. In both cases, $E(t)\u2264Eext(t)$, 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)].

## VI. CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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

This Appendix summarizes our conventions for the Weyl transform. (For more information, see the excellent reviews in Refs. 10 and 41.) The Weyl symbol *A*(**x**, **p**) of any given operator $A\u0302$ is defined as

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 $A\u0302$, its matrix elements in the coordinate representation, $A(x,x\u2032)\u2250\u27e8x|A\u0302|x\u2032\u27e9$, can be expressed as

so *A*(**x**, **p**) can be understood as a spectrum of $A(x,x\u2032)$. 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 $A\u0302$, its trace $Tr\u2009A\u0302\u2250\u222bdx\u2009\u27e8x|A\u0302|x\u27e9$ can be expressed as

If

*A*(**x**,**p**) is the Weyl symbol of $A\u0302$, then*A**(**x**,**p**) is the Weyl symbol of $A\u0302\u2020$. As a corollary, the Weyl symbol of a Hermitian operator is real.- For any $C\u0302=A\u0302B\u0302$, the corresponding Weyl symbols satisfy
^{23,24}(A5)$C(x,p)=A(x,p)\u22c6B(x,p).$Here, ⋆ is the*Moyal product*, which is given by(A6)$A(x,p)\u22c6B(x,p)\u2250A(x,p)eiL\u0302/2B(x,p),$and $L\u0302$ is the*Janus operator*, which is given by(A7)$L\u0302\u2250\u2202x\u2190\xb7\u2202p\u2192\u2212\u2202p\u2190\xb7\u2202x\u2192={\xb7,\xb7}.$The arrows indicate the directions in which the derivatives act, and $AL\u0302B={A,B}$ is the canonical Poisson bracket; namely,(A8)${A,B}\u2250(\u2202xA)\xb7(\u2202pB)\u2212(\u2202pA)\xb7(\u2202xB).$ The Moyal product is associative; i.e.,

- The anti-symmetrized Moyal product defines the so-called
*Moyal bracket*(A10)${{A,B}}\u2250\u2212i(A\u22c6B\u2212B\u22c6A)=2A\u2009sin(L\u0302/2)B.$Because of the latter equality, this bracket is also called the*sine bracket*. To lowest order,(A11)${{A,B}}\u2243AL\u0302B={A,B}.$ - The symmetrized Moyal product is defined as(A12)$[[A,B]]\u2250A\u22c6B+B\u22c6A=2A\u2009cos(L\u0302/2)B.$Because of the latter equality, this bracket is also called the
*cosine bracket*. To lowest order,(A13)$[[A,B]]\u22432AB.$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)$\u222bdnx\u2009dnp\u2009A\u22c6B=\u222bdnx\u2009dnp\u2009AB.$
As a corollary

(A15)(A15a)$\u222bdnx\u2009dnp\u2009{{A,B}}=0,$(A15b)$\u222bdnx\u2009dnp\u2009[[A,B]]=2\u222bdnx\u2009dnp\u2009AB.$ 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)$1\u0302\u2009\u21d4\u20091,\u2003x\u0302i\u2009\u21d4\u2009xi,\u2003p\u0302i\u2009\u21d4\u2009pi.$If
*f*and*g*are any two functions, then(A20)$f(x\u0302)\u2009\u21d4\u2009f(x),\u2003g(p\u0302)\u2009\u21d4\u2009g(p).$Similarly, using Eq. (A6), we have(A21)$f(x\u0302)p\u0302j\u2009\u21d4\u2009pjf(x)+(i/2)\u2202jf(x),$(A22)$p\u0302jf(x\u0302)\u2009\u21d4\u2009pjf(x)\u2212(i/2)\u2202jf(x).$One may also notice the connection between these relations and the commutation relation between operators; i.e., $[x\u0302j,p\u0302k]=x\u0302jp\u0302k\u2212p\u0302kx\u0302j=i\delta kj$.

### 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 $pD,\xb1q2\u2250pD2(p\xb1eyq/2)$. This leads to

Similarly,

so one obtains

Also, using Eq. (A18), we obtain

where $Ar,\xb1s\u2250Ar(p\xb1eys/2,t)$ for any *A _{r}*(

**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

First, consider the Wigner–Moyal model [Eqs. (25) and (26)]. To show conservation of enstrophy, we obtain

where we used Eqs. (A15). To evaluate the remaining terms, we use the Fourier representations of $W\xaf(y,p,t)$ and *U*(*y*, *t*) as defined via Eq. (B1). Specifically, after substituting Eq. (26b) for Γ, we obtain

The Moyal products and the brackets can be evaluated using Eqs. (A16) and (A17). Using the notation $A\xb1q\u2250A(p\xb1eyq/2)$ for any *A*(**p**), one then obtains

To show conservation of energy, we obtain

Here, we used the fact that the Taylor expansion of Eq. (A10) for the Moyal bracket ${{W\xaf,\beta px/pD2}}/pD2$ 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

Now let us consider the WKE model [Eqs. (41) and (42)]. To show the conservation of enstrophy, we obtain

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 $Zdw$, whereas $Zzf$ is obviously not conserved; thus, $Zdw+Zzf$ cannot be conserved either.

To show conservation of energy, we obtain

where the integral of ${W\xaf,\beta px/pD2}/pD2$ 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.

## References

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 $w\u0303(x,t)$ is real rather than complex. However, that is just a matter of initial conditions. The equation governing $|w\u0303\u27e9$ 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 $Zzf$ 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 (*L*_{D} → ∞), 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 $W\xaf=0$ for *p* < *p*_{0}, then we can estimate $\u2202pH\u2272H/p0$. Then, a GO approximation is possible if $(\lambda zfp0)\u22121\u226a1$.

The fundamental reason for this discrepancy is that the tWKE assumes $W\xaf$ to be the true distribution probability of driftons, but $W\xaf$ 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 $\u2202tU\u2033$ is negligible is discussed in Ref. 11.