Zonal-flow dynamics from a phase-space perspective

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⁹ and, in fact, is oversimplified even as a geometrical-optics¹⁰ (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).

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×∇ψ$ 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 ≐(∇2−LD−2α̂)ψ$⁠, where $α̂$ is an operator such that $α̂=1$ in parts of the spectrum corresponding to DWs and $α̂=0$ in those corresponding to ZFs. (The symbol $≐$ 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.²¹ 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¯≐∫dx g/Lx$⁠, where L_x, henceforth assumed equal to one, is the system length along x axis.) In particular, $w=w¯+w̃$⁠, where the two components of the generalized vorticity are¹⁸ 

where $∇D2 ≐ ∇2−LD−2$⁠. Equations for $w̃$ and $w¯$ are obtained by taking the zonal-average and fluctuating parts of Eq. (1). This gives

where $fNL≐ṽ·∇w̃−ṽ·∇w̃¯$ is a term nonlinear with respect to fluctuations. As discussed in Ref. 15, this term represents “eddy-eddy” interactions and is responsible for the Batchelor–Kraichnan inverse-energy cascade; however, it is inessential for the formation of ZFs. Since the main scope of this paper is to specifically study the interaction between eddies and ZFs, we ignore eddy-eddy interactions so f_NL will be neglected. Hence, Equations (4) compose the well-known quasilinear model.¹² 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¯=exU$⁠, whose only component U(y, t) is $U=−∂yψ¯$⁠. Specifically, one has $ṽ·∇w¯=−(∂xψ̃)(∂y2U),v¯·∇w̃=U∂xw̃$⁠, and $ṽ·∇w̃¯=−∂y2 ṽxṽy¯$⁠. We will also assume $Q̃=ξ̃−μdww̃$ and $Q¯=−μzfw¯$⁠. 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.

Consider a family of all reversible linear transformations of $w̃(x,t)$ of the form $∫d2x′ K(x,x′,t)w̃(x′,t)$⁠. These transformations map $w̃(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̃⟩$⁠. (Analogous definitions will be assumed also for $|ψ̃⟩$ and $|ξ̃⟩$⁠.) The original function $w̃(x,t)$ is then understood as a projection of $|w̃⟩$⁠, namely, as its “coordinate representation” given by $w̃(x,t)=⟨x|w̃⟩$⁠. Here, $|x⟩$ are the eigenstates of the position operator $x̂$ normalized such that $⟨x′|x̂|x⟩=x⟨x′|x⟩=x δ(x′−x)$⁠. This definition of a field is similar to that used in quantum mechanics for describing probability amplitudes.²⁷ Hence, it is convenient to describe the dynamics of $|w̃⟩$ using a quantumlike formalism. This is done as follows.

In addition to the coordinate operator $x̂$⁠, we introduce a momentum (wave-vector) operator $p̂$ such that, in the coordinate representation, $p̂≐−i∇$⁠. Accordingly, $|w̃⟩=−p̂D2|ψ̃⟩$⁠, where

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

The operator $Ĥ$ is given by

Also, $Û≐U(ŷ,t)$⁠, and the prime above U henceforth denotes ∂_y; in particular, $Û″≐∂y2 U(ŷ,t)$⁠.

Let us express Eq. (9) as $Ĥ=ĤH+iĤA$⁠, where $ĤH≐(Ĥ+Ĥ†)/2$ and $ĤA≐(Ĥ−Ĥ†)/(2i)$ are the Hermitian and anti-Hermitian parts of $Ĥ$⁠, correspondingly. Explicitly, these operators can be written as

where $[·,·]−$ denotes the commutator given by $[Â,B̂]−=ÂB̂−B̂Â$ and $[·,·]+$ denotes the anti-commutator given by $[Â,B̂]+=ÂB̂+B̂Â$⁠. Let us also introduce a Hermitian operator $Ŵ≐|w̃⟩⟨w̃|$⁠, which, by analogy with quantum mechanics, is interpreted as the “fluctuating-vorticity density” operator. It is seen from Eq. (8) that $Ŵ$ satisfies

where $F̂≐|ξ̃⟩⟨w̃|+|w̃⟩⟨ξ̃|$⁠. In particular, taking the trace of this equation also gives an equation for the “total number of DW quanta,” $N≐Tr Ŵ=∫d2x ⟨x|Ŵ|x⟩=∫d2x w̃2=⟨w̃|w̃⟩$⁠; namely,

This indicates that $ĤA$ determines the loss of quanta, or dissipation of DWs. [In particular, the term $μdw$ in Eq. (10b) is responsible for DW dissipation to the external environment, whereas the term $[Û″,p̂xp̂D−2]−/(2i)$ destroys DW quanta while conserving the total enstrophy, as will be discussed in Sec. III E.] Also, $ĤH$ 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.)

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,²⁸ 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π)² represents the phase-space probability density of driftons [the numerical coefficient (2π)² comes from Eq. (A4)], while beyond the GO limit, it can be considered as a generalization of this probability density.²⁹ Using the fact that $w̃(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 ${{·,·}}$ and $[[·,·]]$ 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, ĤA$⁠, and $F̂$⁠, respectively. In particular, using Eq. (A5) and the fact that U is independent of x, one obtains

where $pD2≐p2+LD−2$⁠. 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¯=W¯(y,p,t)$⁠. 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 $F¯=⟨⟨F⟩⟩$⁠, consider integrating Eq. (8) on a time interval (t₀, t). The result can be written as $|w̃t⟩=|w̃t0⟩+|δw̃t⟩+∫t0tdt′|ξ̃t′⟩,$ where the indexes denote the times at which functions are evaluated, and $|δw̃t⟩≐−i∫t0tdt′Ĥ|w̃t′⟩$⁠. We assume

where Ξ is a correlation function that is homogeneous in x but not necessarily in y.^15,20 Since $|δw̃t⟩$ can be affected by $|ξ̃t′⟩$ only if $t>t′$⁠, one has $⟨⟨|ξ̃t⟩⟨δw̃t|⟩⟩=0$⁠. Hence,

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

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

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¯$⁠, Eq. (6b) becomes

Since $W¯$ is independent of x and satisfies the condition (15), Eq. (23) can also be written as

The combination of Eqs. (19) and (24) forms a closed set of equations that can be used to calculate the dynamics of $W¯$ and U self-consistently.

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

As a reminder, $W¯(y,p,t)$ is the zonal-averaged spectral (or Wigner) function that describes the DW turbulence, and U(y, t) is the ZF velocity. Also, $F¯=F¯(y,p)$ is determined by the correlation function of the external noise $ξ̃$ (Sec. III C). We have also introduced $H≐HH$ and $Γ≐HA+μdw$⁠, or, explicitly In  Appendix B, we also present a spectral representation of these equations that can be used for a numerical implementation of the Wigner–Moyal formulation.

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¯=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:

where we used Eqs. (A4) and (A14) to derive the second set of equalities. According to Eqs. (28a) and (A4), note that the DW enstrophy $Zdw$ and the total number of DW quanta $N≐Tr Ŵ$ are the same up to a constant factor.

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.

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π)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 substituted Eqs. (26). The brackets can be calculated using Eqs. (A17). Hence, we obtain

where we assume the notation $A±q≐A(p±eyq/2)$ for any A. Solving for $W¯q$ 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¯q$⁠, one gets

As expected, this dispersion relation coincides with that obtained using the CE2 formalism.¹⁵ Notably, the dependence of the integrand on $W±q$ makes the expression similar to dispersion relations that emerge in quantum mechanics; for instance, cf. Ref. 31, Sec. 40.

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.³² 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≐2βkxky/kD4$ is the DW group velocity in the absence of ZFs, $kD,±q2≐kx2+(ky±q/2)2+LD−2$⁠, and $X≐(γ+2μ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 $μdw,zf≪|γ|$ 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=−1+4σ±4[σ(σ−1)]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|≲0.33$⁠. In the interval 0 < σ < 1 corresponding to $0.33≲|ky|≲0.82$⁠, γ is complex-valued. In the interval −1/8 ≤ σ ≤ 0, which corresponds to $0.82≲|ky|≲1.07$⁠, the solutions are purely imaginary. Finally, in the interval σ < −1/8 corresponding to $|ky|≳1.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.

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 $∂pH$⁠, which is realized at $p∼LD−1$⁠.³⁷) This gives

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

where ${·,·}$ is the canonical Poisson bracket (A8) and One may recognize Eq. (41a) as a variation of the WKE, so we attribute Eqs. (41) and (42) as the 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¯$ is quadratic in the DW amplitude.] In other words, $ω(y,p,t)≐H+iΓ−iμ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,

Although the difference is only in the high-order derivatives of U, these terms remain important for a number of reasons. For example, in the Hamiltonian $H, U″$ can be comparable to β (as is sometimes the case in geophysics²). Also, consider the following. In isolated systems, the tWKE is $∂tW¯={H,W¯}$⁠, 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.³⁸ 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‴$ and $U″$ 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¯$ injects into the DW-ZF system; namely, $Zext=(t/2)(2π)−2∫dy d2p F¯$⁠. 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π)−2∫dy d2p F¯/pD2$⁠. In both cases, $E(t)≤Eext(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)].

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.

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 $Â$ 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 $Â$⁠, its matrix elements in the coordinate representation, $A(x,x′)≐⟨x|Â|x′⟩$⁠, can be expressed as

so A(x, p) can be understood as a spectrum of $A(x,x′)$⁠. 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 $Tr Â≐∫dx ⟨x|Â|x⟩$ 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 $Ĉ=ÂB̂$⁠, the corresponding Weyl symbols satisfy^23,24

  $C(x,p)=A(x,p)⋆B(x,p).$

  (A5)
  Here, ⋆ is the Moyal product, which is given by

  $A(x,p)⋆B(x,p)≐A(x,p)eiL̂/2B(x,p),$

  (A6)
  and $L̂$ is the Janus operator, which is given by

  $L̂≐∂x←·∂p→−∂p←·∂x→={·,·}.$

  (A7)
  The arrows indicate the directions in which the derivatives act, and $AL̂B={A,B}$ is the canonical Poisson bracket; namely,

  ${A,B}≐(∂xA)·(∂pB)−(∂pA)·(∂xB).$

  (A8)

• The Moyal product is associative; i.e.,