Wave kinetic equation for inhomogeneous drift-wave turbulence beyond the quasilinear approximation

Abstract

The formation of zonal flows from inhomogeneous drift-wave (DW) turbulence is often described using statistical theories derived within the quasilinear approximation. However, this approximation neglects wave–wave collisions. Hence, some important effects such as the Batchelor–Kraichnan inverse-energy cascade are not captured within this approach. Here we derive a wave kinetic equation that includes a DW collision operator in the presence of zonal flows. Our derivation makes use of the Weyl calculus, the quasinormal statistical closure and the geometrical-optics approximation. The obtained model conserves both the total enstrophy and energy of the system. The derived DW collision operator breaks down at the Rayleigh–Kuo threshold. This threshold is missed by homogeneous-turbulence theory but expected from a full-wave quasilinear analysis. In the future, this theory might help better understand the interactions between drift waves and zonal flows, including the validity domain of the quasilinear approximation that is commonly used in the literature.

1 Introduction

The interaction between drift-wave (DW) turbulence and zonal flows (ZFs) has been widely studied in plasma physics (Biglari, Diamond & Terry 1990; Lin 1998; Dorland et al. 2000; Jenko et al. 2000; Diamond et al. 2005; Fujisawa 2009; Connaughton, Nazarenko & Quinn 2015). In the context of magnetic fusion (Conway et al. 2005; Fujisawa 2009; Hillesheim et al. 2016), the spontaneous emergence of ZFs significantly affects the transport of energy, momentum and particles. Understanding this phenomenon is critical to improving plasma confinement, but modelling the underlying physics is difficult. For example, direct numerical simulations of interacting DWs and ZFs strongly depend on the initial conditions and the external random forcing. Thus, statistical methods have been useful and are widely applied in the DW turbulence research, even at the cost of introducing approximations.

One particular statistical approach is the so-called quasilinear (QL) approximation (Farrell & Ioannou 2003), where the ZF equation is kept nonlinear and the equation for DWs is linearized. Among statistical QL theories, the wave kinetic equation (WKE) is a popular model that captures the essential basic physics of DW turbulence, e.g. the formation of ZFs (Diamond et al. 1994; Smolyakov & Diamond 1999; Smolyakov, Diamond & Malkov 2000; Malkov & Diamond 2001; Malkov, Diamond & Rosenbluth 2001; Kaw, Singh & Diamond 2002; Kim & Diamond 2003; Trines et al. 2005; Diamond et al. 2005; Singh et al. 2014; Parker 2016; Ruiz et al. 2016; Parker 2018; Zhu, Zhou & Dodin 2018a , b ; Zhu et al. 2018c ). The WKE has the intuitive form of the Liouville equation for the DW action density $J$ in the ray phase space (Parker 2016; Ruiz et al. 2016; Zhu et al. 2018c ; Parker 2018; Zhu et al. 2018a , b ):

(1.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}J+\{J,\unicode[STIX]{x1D6FA}\}=2\unicode[STIX]{x1D6E4}J, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}$ is the local DW frequency, $\unicode[STIX]{x1D6E4}$ is a dissipation rate due to interactions with ZFs and $\{\cdot ,\cdot \}$ is the canonical Poisson bracket. (For the sake of clarity, terms related to external forcing and dissipation are omitted here.) However, as with all QL models, equation (1.1) neglects nonlinear wave–wave scattering, and in consequence, is not able to capture the Batchelor–Kraichnan inverse-energy cascade (Srinivasan & Young 2012) or produce the Kolmogorov–Zakharov spectra for DWs (Connaughton et al. 2015). Hence, a question remains as to whether the existing WKE for inhomogeneous turbulence can be complemented with a wave–wave collision operator $C[J,J]$ . The goal of this work is to calculate $C[J,J]$ explicitly.

Starting from the generalized Hasegawa–Mima equations (gHME) (Smolyakov & Diamond 1999; Krommes & Kim 2000), we derive a WKE with a DW collision operator $C[J,J]$ for inhomogeneous DW turbulence. Our derivation is based on the Weyl calculus (Weyl 1931), which makes our approach similar to that in Ruiz et al. (2016) where the QL approximation was used. However, in contrast to Ruiz et al. (2016), we do not rely on the QL approximation here but instead account for DW collisions perturbatively using the quasinormal approximation. The main result of this work are (4.12) and (4.22). In this final result, DWs are modelled in a similar manner as in (1.1). The difference is that (4.12) includes nonlinear DW scattering, which is described by a wave–wave collision operator $C[J,J]$ that is bilinear in the DW wave-action density $J$ . The resulting model conserves the two nonlinear invariants of the gHME, which are the total enstrophy and the total energy.

The present formulation is fundamentally different from the previously reported homogeneous weak-wave turbulence models for DW turbulence (Krommes 2002; Connaughton et al. 2015). While DWs are described as an incoherent fluctuating field as usual, ZFs are now treated as coherent structures, which are missed in homogeneous-turbulence theory. The obtained model motivates future investigations of the effects of nonlinear wave–wave scattering on DW–ZF interactions, in particular, the spontaneous emergence of ZFs and the eventual saturation of the ZFs and the DW spectra. This theory might also help better understand the validity domain of the quasilinear approach to DW turbulence that has been commonly used in the literature.

The present work is organized as follows. In § 2, we introduce the gHME and obtain the governing equations for the mean and fluctuating components of the fields. In § 3, we introduce the quasinormal statistical closure to obtain a closed equation for the correlation operator describing the vorticity fluctuations. In § 4, we project the equations into the DW-ray phase space and obtain the collisional WKE. In § 5, we discuss the obtained equations, their conservation properties and their relation to previous models. Final conclusions and remarks are given in § 6. In appendix A, we give a brief introduction to the Weyl calculus and also define the zonal average of an arbitrary operator. In appendix B, we present some auxiliary calculations.

2 Basic equations

2.1 The generalized Hasegawa–Mima model

We consider a magnetized plasma in a uniform magnetic field in the $z$ direction and with an equilibrium local gradient of the plasma density. Upon assuming a quasi-adiabatic response for the electrons and a fluid description for the ion dynamics, we model electrostatic two-dimensional (2-D) turbulent flows in the $xy$ plane using the gHME (Smolyakov & Diamond 1999; Krommes & Kim 2000), which in normalized units is given by

(2.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}w+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}w+\unicode[STIX]{x1D6FD}\,\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713}=Q. & & \displaystyle\end{eqnarray}$$

Here $\boldsymbol{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 background plasma density, which is measured by the constant $\unicode[STIX]{x1D6FD}$ . In (2.1), time is measured in units of the inverse ion cyclotron frequency while length is measured in units of the ion sound radius. (See, e.g. Zhu et al. (2018c ) for more details.) The function $\unicode[STIX]{x1D713}(t,\boldsymbol{x})$ is the electric potential, $\boldsymbol{v}\doteq \boldsymbol{e}_{z}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ is the ion fluid velocity on the $xy$ plane and $\boldsymbol{e}_{z}$ is a unit vector normal to this plane. (The symbol $\doteq$ denotes definitions.) The generalized vorticity $w(t,\boldsymbol{x})$ is given by

(2.2) $$\begin{eqnarray}\displaystyle w(t,\boldsymbol{x})\doteq (\unicode[STIX]{x1D6FB}^{2}-L_{\text{D}}^{-2}\widehat{a})\unicode[STIX]{x1D713}, & & \displaystyle\end{eqnarray}$$

where $\widehat{a}$ is an operator such that $\widehat{a}=1$ in parts of the spectrum corresponding to DWs and $\widehat{a}=0$ in those corresponding to ZFs. The constant $L_{\text{D}}$ is the ion sound radius. ( $L_{\text{D}}=1$ in normalized units.) The term $Q(t,\boldsymbol{x})$ in (2.1) represents external forces and dissipation. Notably, equation (2.1) with $L_{\text{D}}\rightarrow \infty$ also describes Rossby turbulence in planetary atmospheres (Farrell & Ioannou 2003, 2007; Marston, Conover & Schneider 2008; Srinivasan & Young 2012; Ait-Chaalal et al. 2016).

For isolated systems, where $Q=0$ , equation (2.1) conserves two nonlinear invariants: the enstrophy ${\mathcal{Z}}$ and the energy ${\mathcal{E}}$ (strictly speaking, free energy). These are defined as

(2.3a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{Z}}(t)\doteq \frac{1}{2}\int \text{d}^{2}\boldsymbol{x}\,w^{2},\quad {\mathcal{E}}(t)\doteq -\frac{1}{2}\int \text{d}^{2}\boldsymbol{x}\,w\unicode[STIX]{x1D713}. & & \displaystyle\end{eqnarray}$$

The statistical model that we shall derive conserves both of these nonlinear invariants.

2.2 Separating the mean and fluctuating components of the fields

Let us decompose the fields $\unicode[STIX]{x1D713}$ and $w$ into their mean and fluctuating components, which will be denoted by bars and tildes, respectively. For any arbitrary field $g(t,\boldsymbol{x})$ , the mean part is defined as $\overline{g}(t,y)\doteq \int \text{d}x\,\boldsymbol{\langle }\!\langle g\rangle \!\boldsymbol{\rangle }/L_{x}$ , where $L_{x}$ is the system’s length along $x$ and $\boldsymbol{\langle }\!\langle \cdot \rangle \!\boldsymbol{\rangle }$ denotes the statistical average over realizations of initial conditions or of the external random forcing. The mean part will describe the coherent ZF dynamics. In contrast, the fluctuating quantities will describe the incoherent DW dynamics.

We consider the fluctuating quantities to be small in amplitude. The small-amplitude ordering is denoted by using the small dimensionless parameter $\unicode[STIX]{x1D716}_{\text{nl}}\ll 1$ . In particular,

(2.4a,b ) $$\begin{eqnarray}\displaystyle w=\overline{w}+\unicode[STIX]{x1D716}_{\text{nl}}\,\widetilde{w},\quad \boldsymbol{v}=U\boldsymbol{e}_{x}+\unicode[STIX]{x1D716}_{\text{nl}}\widetilde{\boldsymbol{v}}, & & \displaystyle\end{eqnarray}$$

where $U(t,y)\doteq -\unicode[STIX]{x2202}_{y}\overline{\unicode[STIX]{x1D713}}$ is the ZF velocity field and the two components of the generalized vorticity are related to $\unicode[STIX]{x1D713}$ as (Parker & Krommes 2013)

(2.5a,b ) $$\begin{eqnarray}\displaystyle \overline{w}=\unicode[STIX]{x1D6FB}^{2}\overline{\unicode[STIX]{x1D713}},\quad \widetilde{w}=\unicode[STIX]{x1D6FB}_{\text{D}}^{2}\widetilde{\unicode[STIX]{x1D713}}. & & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D6FB}_{\text{D}}^{2}\doteq \unicode[STIX]{x1D6FB}^{2}-1$ . Also, $\overline{\widetilde{w}}=0$ and $\overline{\widetilde{\boldsymbol{v}}}=0$ by definition.

From (2.1), we then derive the governing equations for the fluctuating and the mean fields (Srinivasan & Young 2012):

(2.6a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\widetilde{w}+U\unicode[STIX]{x2202}_{x}\widetilde{w}+[\unicode[STIX]{x1D6FD}-(\unicode[STIX]{x2202}_{y}^{2}U)]\unicode[STIX]{x2202}_{x}\widetilde{\unicode[STIX]{x1D713}}+\unicode[STIX]{x1D716}_{\text{nl}}\,f_{\text{eddy}}=\widetilde{Q}, & \displaystyle\end{eqnarray}$$

(2.6b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}U+\unicode[STIX]{x1D707}_{\text{zf}}U+\unicode[STIX]{x1D716}_{\text{nl}}^{2}\,\unicode[STIX]{x2202}_{y}\overline{\widetilde{v}_{x}\widetilde{v}_{y}}=\overline{Q}, & \displaystyle\end{eqnarray}$$

where $f_{\text{eddy}}(t,\boldsymbol{x})\doteq \widetilde{\boldsymbol{v}}\cdot \unicode[STIX]{x1D735}\widetilde{w}-\overline{\widetilde{\boldsymbol{v}}\cdot \unicode[STIX]{x1D735}\widetilde{w}}$ is a nonlinear term representing eddy–eddy interactions and is responsible for the Batchelor–Kraichnan inverse-energy cascade (Srinivasan & Young 2012). As seen from (2.6), the parameter $\unicode[STIX]{x1D716}_{\text{nl}}$ denotes the smallness of the nonlinear coupling terms. Also note that we considered the ordering $Q=\overline{Q}+\unicode[STIX]{x1D716}_{\text{nl}}\,\widetilde{Q}$ .

2.3 Temporal and spatial scale separation of fluctuating and mean fields

In the model derived below, we shall assume that there is a temporal and spatial scale separation between the DW and ZF fields. Specifically, let $\unicode[STIX]{x1D70F}_{\text{dw}}$ and $\unicode[STIX]{x1D706}_{\text{dw}}$ respectively denote the characteristic period and wavelength of the DWs. In a similar manner, the characteristic time and length scales of the ZFs are given by $T_{\text{zf}}$ and $L_{\text{zf}}$ , respectively. Upon following the discussion in Ruiz et al. (2016), we shall characterize the scale separation between DWs and ZFs by introducing the geometrical-optics (GO) parameter

(2.7) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}\doteq \text{max}\left(\frac{\unicode[STIX]{x1D70F}_{\text{dw}}}{T_{\text{zf}}},\frac{\unicode[STIX]{x1D706}_{\text{dw}}}{L_{\text{zf}}},\frac{L_{\text{D}}}{L_{\text{zf}}}\right)\ll 1. & & \displaystyle\end{eqnarray}$$

In addition, we assume that the DWs are weakly damped and weakly dissipated so that $\widetilde{Q}=\unicode[STIX]{x1D716}\,\widetilde{\unicode[STIX]{x1D709}}-\unicode[STIX]{x1D716}\,\unicode[STIX]{x1D707}_{\text{dw}}\widetilde{w}$ , where $\widetilde{\unicode[STIX]{x1D709}}$ is some white-noise external forcing with zero mean and $\unicode[STIX]{x1D707}_{\text{dw}}$ is intended to emulate the dissipation of DWs caused by the external environment. The external dissipation and random forcing are scaled using the GO parameter (McDonald & Kaufman 1985). For the mean quantities, we consider $\overline{Q}=-\unicode[STIX]{x1D707}_{\text{zf}}\,\overline{w}$ , where $\unicode[STIX]{x1D707}_{\text{zf}}$ represents dissipation on the mean flows.

Upon using the orderings above, we write (2.6) as

(2.8a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}\unicode[STIX]{x2202}_{t}\widetilde{w}+\unicode[STIX]{x1D716}U\unicode[STIX]{x2202}_{x}\widetilde{w}+\unicode[STIX]{x1D716}[\unicode[STIX]{x1D6FD}-(\unicode[STIX]{x2202}_{y}^{2}U)]\unicode[STIX]{x2202}_{x}\widetilde{\unicode[STIX]{x1D713}}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D707}_{\text{dw}}\widetilde{w}+\unicode[STIX]{x1D716}_{\text{nl}}f_{\text{eddy}}=\unicode[STIX]{x1D716}\widetilde{\unicode[STIX]{x1D709}}, & \displaystyle\end{eqnarray}$$

(2.8b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}U+\unicode[STIX]{x1D707}_{\text{zf}}U+\unicode[STIX]{x1D716}_{\text{nl}}^{2}\,\unicode[STIX]{x2202}_{y}\overline{\widetilde{v}_{x}\widetilde{v}_{y}}=0, & \displaystyle\end{eqnarray}$$

where the GO parameter $\unicode[STIX]{x1D716}$ was inserted in front of the partial derivatives acting on fluctuating quantities in order to denote the small-scale ordering. (Here $\unicode[STIX]{x1D716}$ serves the same role as $\hbar$ in quantum mechanics.) Note that $\widetilde{\boldsymbol{v}}=\boldsymbol{e}_{z}\times \unicode[STIX]{x1D716}\unicode[STIX]{x1D735}\widetilde{\unicode[STIX]{x1D713}}$ , $\widetilde{w}=(\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FB}^{2}-1)\widetilde{\unicode[STIX]{x1D713}}$ and $f_{\text{eddy}}=\widetilde{\boldsymbol{v}}\cdot \unicode[STIX]{x1D716}\unicode[STIX]{x1D735}\widetilde{w}-\overline{\widetilde{\boldsymbol{v}}\cdot \unicode[STIX]{x1D716}\unicode[STIX]{x1D735}\widetilde{w}}$ . Also, note $\unicode[STIX]{x2202}_{y}\overline{\widetilde{v}_{x}\widetilde{v}_{y}}\sim O(\unicode[STIX]{x1D716}^{0})$ since $\overline{\widetilde{v}_{x}\widetilde{v}_{y}}$ is a mean quantity. Equations (2.8) govern the dynamics of the fluctuating vorticity $\widetilde{w}(t,\boldsymbol{x})$ and of the ZF velocity $U(t,y)$ . Although we have presented our ordering assumptions in § 2.2 and § 2.3, no approximations have been adopted yet, so these equations are equivalent to the gHME.

2.4 Abstract vector representation

Let us write the fluctuating fields as elements of an abstract Hilbert space $L^{2}(\mathbb{R}^{3})$ of wave states with inner product (Dodin 2014; Littlejohn & Winston 1993)

(2.9) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D719}\mid \unicode[STIX]{x1D713}\rangle =\int \text{d}t\,\text{d}^{2}\boldsymbol{x}\,\unicode[STIX]{x1D719}^{\ast }(t,\boldsymbol{x})\,\unicode[STIX]{x1D713}(t,\boldsymbol{x}), & & \displaystyle\end{eqnarray}$$

where the integrals are taken over $\mathbb{R}^{3}$ . Let $|t,\boldsymbol{x}\rangle$ be the eigenstates of the time and position operators $\widehat{t}$ and $\widehat{\boldsymbol{x}}$ . (Considering $\widehat{t}$ as an operator will allow us to write the statistical closure in § 3 in a simple manner.) Hence, $\widetilde{w}(t,\boldsymbol{x})$ is written as $\widetilde{w}(t,\boldsymbol{x})=\langle t,\boldsymbol{x}\mid \widetilde{w}\rangle$ . Since $\widetilde{w}(t,\boldsymbol{x})$ is real, then $\langle \widetilde{w}\mid t,\boldsymbol{x}\rangle =\langle t,\boldsymbol{x}\mid \widetilde{w}\rangle$ . In the following, we shall use the notation $\unicode[STIX]{x1D639}\doteq (t,\boldsymbol{x})$ to denote space–time coordinates. Likewise, $|\unicode[STIX]{x1D639}\rangle \doteq |t,\boldsymbol{x}\rangle$ is the corresponding space–time eigenstate, and $\text{d}^{3}\unicode[STIX]{x1D639}\doteq \text{d}t\,\text{d}^{2}\boldsymbol{x}$ is the space–time differential volume.

In addition to the time and position operators, we introduce the frequency operator $\widehat{\unicode[STIX]{x1D714}}$ , such that $\widehat{\unicode[STIX]{x1D714}}\doteq \text{i}\unicode[STIX]{x1D716}\unicode[STIX]{x2202}_{t}$ in the coordinate representation. Likewise, the wavevector operator $\widehat{\boldsymbol{k}}$ is defined as $\widehat{\boldsymbol{k}}\doteq -\text{i}\unicode[STIX]{x1D716}\unicode[STIX]{x1D735}$ . In particular, these operators are Hermitian, and one has $\langle \unicode[STIX]{x1D639}\mid \widehat{\unicode[STIX]{x1D714}}\mid \widetilde{w}\rangle =\text{i}\unicode[STIX]{x1D716}\unicode[STIX]{x2202}_{t}\widetilde{w}$ and $\langle \unicode[STIX]{x1D639}\mid \widehat{\boldsymbol{k}}\mid \widetilde{w}\rangle =-\text{i}\unicode[STIX]{x1D716}\unicode[STIX]{x1D735}\widetilde{w}$ . Using the relation between the fluctuating generalized vorticity and electric potential, one has $|\widetilde{w}\rangle =-\widehat{k}_{\text{D}}^{2}|\widetilde{\unicode[STIX]{x1D713}}\rangle$ , where

(2.10a,b ) $$\begin{eqnarray}\displaystyle \widehat{k}_{\text{D}}^{2}\doteq \widehat{k}^{2}+1,\quad \widehat{k}^{2}\doteq \widehat{\boldsymbol{k}}\cdot \widehat{\boldsymbol{k}}. & & \displaystyle\end{eqnarray}$$

Let us write the fluctuating quantities appearing in (2.8) in terms of the abstract states and operators. One finds that (2.8a ) can be written as

(2.11) $$\begin{eqnarray}\displaystyle \widehat{\mathscr{D}}|\widetilde{w}\rangle =\frac{\text{i}\unicode[STIX]{x1D716}_{\text{nl}}}{2}|f_{\text{nl}}[\widetilde{w},\widetilde{w}]\rangle -\frac{\text{i}\unicode[STIX]{x1D716}_{\text{nl}}}{2}|\overline{f_{\text{nl}}[\widetilde{w},\widetilde{w}]}\rangle +\text{i}\unicode[STIX]{x1D716}|\widetilde{\unicode[STIX]{x1D709}}\rangle . & & \displaystyle\end{eqnarray}$$

Here $\widehat{\mathscr{D}}$ is the DW dispersion operator

(2.12) $$\begin{eqnarray}\displaystyle \widehat{\mathscr{D}}\doteq \widehat{\unicode[STIX]{x1D714}}-\widehat{U}\widehat{k}_{x}+(\unicode[STIX]{x1D6FD}-\widehat{U}^{\prime \prime })\widehat{k}_{x}\widehat{k}_{\text{D}}^{-2}+\text{i}\unicode[STIX]{x1D716}\unicode[STIX]{x1D707}_{\text{dw}}, & & \displaystyle\end{eqnarray}$$

where $\widehat{U}\doteq U(\widehat{t},\widehat{y})$ , and the prime above $U$ henceforth denotes $\unicode[STIX]{x2202}_{y}$ ; in particular, $\widehat{U}^{\prime \prime }\doteq \unicode[STIX]{x2202}_{y}^{2}\,U(\widehat{t},\widehat{y})$ . Also, $|\widetilde{\unicode[STIX]{x1D709}}\rangle$ is the ket corresponding to the random forcing $\widetilde{\unicode[STIX]{x1D709}}$ . It is to be noted that the $\widehat{\mathscr{D}}$ includes the nonlinear coupling between the ZFs and the DWs.

Additionally in (2.11), the kets $|\,f_{\text{nl}}[\unicode[STIX]{x1D719},\unicode[STIX]{x1D713}]\rangle$ and $|\,\overline{f_{\text{nl}}[\unicode[STIX]{x1D719},\unicode[STIX]{x1D713}]}\rangle$ are given by

(2.13) $$\begin{eqnarray}\displaystyle & \displaystyle |f_{\text{nl}}[\unicode[STIX]{x1D719},\unicode[STIX]{x1D713}]\rangle \doteq \int \text{d}^{3}\unicode[STIX]{x1D639}\,|\unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D719}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \unicode[STIX]{x1D713}\rangle , & \displaystyle\end{eqnarray}$$

(2.14) $$\begin{eqnarray}\displaystyle & \displaystyle |\overline{f_{\text{nl}}[\unicode[STIX]{x1D719},\unicode[STIX]{x1D713}]}\rangle \doteq \int \text{d}^{3}\unicode[STIX]{x1D639}\,|\unicode[STIX]{x1D639}\rangle \overline{\langle \unicode[STIX]{x1D719}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \unicode[STIX]{x1D713}\rangle }, & \displaystyle\end{eqnarray}$$

where $\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})$ is an operator describing the nonlinear eddy–eddy interactions

(2.15) $$\begin{eqnarray}\displaystyle \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\doteq \widehat{\mathscr{L}}_{j}|\unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D639}|\widehat{\mathscr{R}}_{j}+\widehat{\mathscr{R}}_{j}|\unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D639}|\widehat{\mathscr{L}}_{j}. & & \displaystyle\end{eqnarray}$$

(The summation over repeating indices is assumed.) The operators $\widehat{\mathscr{L}}_{j}$ and $\widehat{\mathscr{R}}_{j}$ are given by

(2.16a,b ) $$\begin{eqnarray}\displaystyle \widehat{\mathscr{L}}_{j}\doteq (\boldsymbol{e}_{z}\times \widehat{\boldsymbol{k}})_{j}k_{\text{D}}^{-2},\quad \widehat{\mathscr{R}}_{j}\doteq \widehat{\boldsymbol{k}}_{j}. & & \displaystyle\end{eqnarray}$$

Indeed, we can verify that $|f_{\text{nl}}[\widetilde{w},\widetilde{w}]\rangle$ represents nonlinear eddy–eddy interactions. When projecting onto the coordinate eigenstates, one obtains

(2.17) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D639}\mid f_{\text{nl}}[\widetilde{w},\widetilde{w}]\rangle & = & \displaystyle \langle \widetilde{w}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}\rangle \nonumber\\ \displaystyle & = & \displaystyle \langle \widetilde{w}\mid \widehat{k}_{\text{D}}^{-2}(\boldsymbol{e}_{z}\times \widehat{\boldsymbol{k}})_{j}\mid \unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D639}\mid \widehat{\boldsymbol{k}}_{j}\mid \widetilde{w}\rangle +\langle \widetilde{w}\mid \widehat{\boldsymbol{k}}_{j}\mid \unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D639}\mid (\boldsymbol{e}_{z}\times \widehat{\boldsymbol{k}})_{j}\widetilde{k}_{\text{D}}^{-2}\mid \widetilde{w}\rangle \nonumber\\ \displaystyle & = & \displaystyle -(\langle \unicode[STIX]{x1D639}\mid (\boldsymbol{e}_{z}\times \widehat{\boldsymbol{k}})_{j}\mid \widetilde{\unicode[STIX]{x1D713}}\rangle )^{\ast }\langle \unicode[STIX]{x1D639}\mid \widehat{\boldsymbol{k}}_{j}\mid \widetilde{w}\rangle -(\langle \unicode[STIX]{x1D639}\mid \widehat{\boldsymbol{k}}_{j}\mid \widetilde{w}\rangle )^{\ast }\langle \unicode[STIX]{x1D639}\mid (\boldsymbol{e}_{z}\times \widehat{\boldsymbol{k}})_{j}\mid \widetilde{\unicode[STIX]{x1D713}}\rangle \nonumber\\ \displaystyle & = & \displaystyle -2(\boldsymbol{e}_{z}\times \unicode[STIX]{x1D716}\unicode[STIX]{x1D735}\widetilde{\unicode[STIX]{x1D713}})\cdot \unicode[STIX]{x1D716}\unicode[STIX]{x1D735}\widetilde{w}\nonumber\\ \displaystyle & = & \displaystyle -2\widetilde{\boldsymbol{v}}\times \unicode[STIX]{x1D716}\unicode[STIX]{x1D735}\widetilde{w},\end{eqnarray}$$

which is one of the nonlinear terms appearing in (2.8a ). Here we substituted the coordinate representation for the wavevector operator and used $|\widetilde{w}\rangle =-\widehat{k}_{\text{D}}^{2}|\widetilde{\unicode[STIX]{x1D713}}\rangle$ . Note that the factor 2 appears because we have written the ket $|f_{\text{nl}}[\unicode[STIX]{x1D719},\unicode[STIX]{x1D713}]\rangle$ in a symmetric form so that

(2.18) $$\begin{eqnarray}\displaystyle |f_{\text{nl}}[\unicode[STIX]{x1D719},\unicode[STIX]{x1D713}]\rangle =|f_{\text{nl}}[\unicode[STIX]{x1D713},\unicode[STIX]{x1D719}]\rangle & & \displaystyle\end{eqnarray}$$

for any two real fields $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D713}$ .

Following Ruiz et al. (2016), the fluctuating terms appearing in (2.8b ) can also be rewritten in the abstract representation. Upon noting that $\widetilde{v}_{x}\widetilde{v}_{y}=(-\unicode[STIX]{x1D716}\unicode[STIX]{x2202}_{y}\widetilde{\unicode[STIX]{x1D713}})(\unicode[STIX]{x1D716}\unicode[STIX]{x2202}_{x}\widetilde{\unicode[STIX]{x1D713}})=-\langle \unicode[STIX]{x1D639}\mid \widehat{k}_{x}\mid \widetilde{\unicode[STIX]{x1D713}}\rangle \langle \widetilde{\unicode[STIX]{x1D713}}\mid \widehat{k}_{y}\mid \unicode[STIX]{x1D639}\rangle$ , one obtains

(2.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}U+\unicode[STIX]{x1D707}_{\text{zf}}U=\unicode[STIX]{x1D716}_{\text{nl}}^{2}\,\unicode[STIX]{x2202}_{y}\overline{\langle \unicode[STIX]{x1D639}\mid \widehat{k}_{x}\widehat{k}_{\text{D}}^{-2}\mid \widetilde{w}\rangle \langle \widetilde{w}\mid \widehat{k}_{\text{D}}^{-2}\widehat{k}_{y}\mid \unicode[STIX]{x1D639}\rangle }. & & \displaystyle\end{eqnarray}$$

3 Statistical closure

3.1 Statistical-closure problem

Let us now introduce the correlation operator for the fluctuating vorticity field:

(3.1) $$\begin{eqnarray}\displaystyle \widehat{\mathscr{W}}\doteq \overline{|\widetilde{w}\rangle \langle \widetilde{w}|}, & & \displaystyle\end{eqnarray}$$

where the Hermitian operator $|\widetilde{w}\rangle \langle \widetilde{w}|$ can be interpreted as the fluctuating-vorticity density operator by analogy with quantum mechanics (Ruiz et al. 2016). It is to be noted that $\widehat{\mathscr{W}}$ is defined as the zonal average of the fluctuating-vorticity density operator. The zonal average of an abstract operator is discussed in § A.2. By using the identity (A 20), one can write (2.19) in terms of the correlation operator:

(3.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}U+\unicode[STIX]{x1D707}_{\text{zf}}U=\unicode[STIX]{x1D716}_{\text{nl}}^{2}\,\unicode[STIX]{x2202}_{y}\langle \unicode[STIX]{x1D639}\mid \widehat{k}_{x}\widehat{k}_{\text{D}}^{-2}\widehat{\mathscr{W}}\widehat{k}_{\text{D}}^{-2}\widehat{k}_{y}\mid \unicode[STIX]{x1D639}\rangle . & & \displaystyle\end{eqnarray}$$

Now, let us obtain the governing equation for $\widehat{\mathscr{W}}$ . Multiplying (2.11) by $\langle \widetilde{w}|$ from the right and zonal averaging leads to an equation for the correlation operator:

(3.3) $$\begin{eqnarray}\displaystyle \widehat{\mathscr{D}}\widehat{\mathscr{W}}=\frac{\text{i}\unicode[STIX]{x1D716}_{\text{nl}}}{2}\,\overline{|f_{\text{nl}}[\widetilde{w},\widetilde{w}]\rangle \langle \widetilde{w}|}+\text{i}\unicode[STIX]{x1D716}\,\overline{|\widetilde{\unicode[STIX]{x1D709}}\rangle \langle \widetilde{w}|}. & & \displaystyle\end{eqnarray}$$

Subtracting from (3.3) its Hermitian conjugate gives

(3.4) $$\begin{eqnarray}\displaystyle [\widehat{\mathscr{D}}_{\text{H}},\widehat{\mathscr{W}}]_{-}+\text{i}[\widehat{\mathscr{D}}_{\text{A}},\widehat{\mathscr{W}}]_{+}=\text{i}\unicode[STIX]{x1D716}_{\text{nl}}\,[\overline{|f_{\text{nl}}[\widetilde{w},\widetilde{w}]\rangle \langle \widetilde{w}|}]_{\text{H}}+2\text{i}\unicode[STIX]{x1D716}\,[\overline{|\widetilde{\unicode[STIX]{x1D709}}\rangle \langle \widetilde{w}|}]_{\text{H}}, & & \displaystyle\end{eqnarray}$$

where the subscripts ‘H’ and ‘A’ denote the Hermitian and anti-Hermitian parts of an arbitrary operator; i.e. $\widehat{\mathscr{A}}_{\text{H}}\doteq (\widehat{\mathscr{A}}+\widehat{\mathscr{A}}^{\dagger })/2$ and $\widehat{\mathscr{A}}_{\text{A}}\doteq (\widehat{\mathscr{A}}-\widehat{\mathscr{A}}^{\dagger })/(2\text{i})$ . Specifically, the operator $\widehat{\mathscr{D}}$ is decomposed as $\widehat{\mathscr{D}}=\widehat{\mathscr{D}}_{\text{H}}+\text{i}\widehat{\mathscr{D}}_{\text{A}}$ , where

(3.5a ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\mathscr{D}}_{\text{H}}=\widehat{\unicode[STIX]{x1D714}}-\widehat{U}\widehat{k}_{x}+[\unicode[STIX]{x1D6FD}-\widehat{U}^{\prime \prime },\widehat{k}_{x}\widehat{k}_{\text{D}}^{-2}]_{+}/2, & \displaystyle\end{eqnarray}$$

(3.5b ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\mathscr{D}}_{\text{A}}=\text{i}[\widehat{U}^{\prime \prime },\widehat{k}_{x}\widehat{k}_{\text{D}}^{-2}]_{-}/2+\unicode[STIX]{x1D716}\unicode[STIX]{x1D707}_{\text{dw}}. & \displaystyle\end{eqnarray}$$

Both $\widehat{\mathscr{D}}_{\text{H}}$ and $\widehat{\mathscr{D}}_{\text{A}}$ are Hermitian. Also, $[\cdot ,\cdot ]_{\mp }$ denote the commutators and anticommutators; i.e. $[\widehat{\mathscr{A}},\widehat{\mathscr{B}}]_{-}=\widehat{\mathscr{A}}\widehat{\mathscr{B}}-\widehat{\mathscr{B}}\widehat{\mathscr{A}}$ and $[\widehat{\mathscr{A}},\widehat{\mathscr{B}}]_{+}=\widehat{\mathscr{A}}\widehat{\mathscr{B}}+\widehat{\mathscr{B}}\widehat{\mathscr{A}}$ . Aside from the additional first term on the right-hand side of (3.4), this result is similar to that obtained in Ruiz et al. (2016).

Equations (3.2) and (3.4) are not closed. The left-hand side of (3.4) is written in terms of $\widehat{\mathscr{W}}$ , which is bilinear in the fluctuating vorticity. However, the right-hand side of (3.4) contains terms that are linear and cubic with respect to $\widetilde{w}$ . This is the fundamental statistical closure problem (Frisch & Kolmogorov 1995; Kraichnan 2013; Krommes 2002). The next step is to introduce a statistical closure in order to express (3.4) in terms of the correlation operator $\widehat{\mathscr{W}}$ only.

3.2 A comment on the quasilinear approximation

One possibility is to neglect the first term on the right-hand side of (3.4). Then, one linearizes (2.11) so that

(3.6) $$\begin{eqnarray}\displaystyle \widehat{\mathscr{D}}|\widetilde{w}\rangle \simeq \text{i}\unicode[STIX]{x1D716}|\widetilde{\unicode[STIX]{x1D709}}\rangle . & & \displaystyle\end{eqnarray}$$

After formally inverting $\widehat{\mathscr{D}}$ , we have $|\widetilde{w}\rangle \simeq \text{i}\unicode[STIX]{x1D716}\widehat{\mathscr{D}}^{-1}|\widetilde{\unicode[STIX]{x1D709}}\rangle$ . Substituting into (3.4) leads to the closed equation

(3.7) $$\begin{eqnarray}\displaystyle [\widehat{\mathscr{D}}_{\text{H}},\widehat{\mathscr{W}}]_{-}+\text{i}[\widehat{\mathscr{D}}_{\text{A}},\widehat{\mathscr{W}}]_{+}=2\text{i}\unicode[STIX]{x1D716}^{2}\,[\widehat{\mathscr{S}}(\widehat{\mathscr{D}}^{-1})^{\dagger }]_{\text{A}}, & & \displaystyle\end{eqnarray}$$

where $\widehat{\mathscr{S}}\doteq \overline{|\widetilde{\unicode[STIX]{x1D709}}\rangle \langle \widetilde{\unicode[STIX]{x1D709}}|}$ is the zonal-averaged density operator associated with the random external forcing. We also used the fact that, for any operator $\widehat{\mathscr{A}}$ , one has $(-\text{i}\,\widehat{\mathscr{A}})_{\text{H}}=\widehat{\mathscr{A}}_{\text{A}}$ .

Equations (3.2) and (3.7) now form a closed system. In this approximation, the equation for the DW fluctuations was linearized in (3.6), but the DW nonlinearity is only kept in the equation for the mean field (3.2). This constitutes the QL approximation. If one projects (3.7) on the double-physical coordinate space $(\unicode[STIX]{x1D639},\unicode[STIX]{x1D639}^{\prime })$ using multiplication by $\langle t,\boldsymbol{x}|$ and $|t,\boldsymbol{x}^{\prime }\rangle$ , then one obtains the so-called CE2 equations (Farrell & Ioannou 2003, 2007; Marston et al. 2008; Srinivasan & Young 2012; Ait-Chaalal et al. 2016). Likewise, if one projects (3.7) on the ray space $(t,\boldsymbol{x},\unicode[STIX]{x1D714},\boldsymbol{k})$ using the Weyl transform (Weyl 1931), then one obtains the Wigner–Moyal model discussed in Ruiz et al. (2016), Parker (2018) and Zhu et al. (2018c ).

3.3 A statistical closure beyond the quasilinear approximation

We extend our theory beyond the QL approximation in order to retain wave–wave scattering, namely, the term $f_{\text{nl}}$ in (3.4). Let us separate $\widetilde{w}$ into three components:

(3.8) $$\begin{eqnarray}\displaystyle |\widetilde{w}\rangle =|\widetilde{w}_{0}\rangle +\unicode[STIX]{x1D716}_{\text{nl}}|\widetilde{\unicode[STIX]{x1D719}}\rangle +\unicode[STIX]{x1D716}|\widetilde{\unicode[STIX]{x1D711}}\rangle . & & \displaystyle\end{eqnarray}$$

Here $|\widetilde{w}_{0}\rangle =O(1)$ is chosen to satisfy the linear part of (2.11); i.e. $\widehat{\mathscr{D}}|\widetilde{w}_{0}\rangle =0$ . The fluctuations in $\widetilde{w}_{0}$ are due to random initial conditions, whose statistics are considered to be uncorrelated to those of the random forcing $\widetilde{\unicode[STIX]{x1D709}}$ . When formally inverting $\widehat{\mathscr{D}}$ in (2.11), one finds that, to lowest order in $\unicode[STIX]{x1D716}_{\text{nl}}$ and $\unicode[STIX]{x1D716}$ ,

(3.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle |\widetilde{\unicode[STIX]{x1D719}}\rangle \simeq \frac{\text{i}}{2}\,\widehat{\mathscr{D}}^{-1}\{|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{w}_{0}]\rangle -\overline{|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{w}_{0}]\rangle }\}, & \displaystyle\end{eqnarray}$$

(3.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle |\widetilde{\unicode[STIX]{x1D711}}\rangle \simeq \text{i}\widehat{\mathscr{D}}^{-1}|\widetilde{\unicode[STIX]{x1D709}}\rangle . & \displaystyle\end{eqnarray}$$

Substituting (3.8) and (3.9) into (3.4) leads to

(3.10) $$\begin{eqnarray}\displaystyle [\widehat{\mathscr{D}}_{\text{H}},\widehat{\mathscr{W}}]_{-}+\text{i}[\widehat{\mathscr{D}}_{\text{A}},\widehat{\mathscr{W}}]_{+} & = & \displaystyle \text{i}\unicode[STIX]{x1D716}_{\text{nl}}\{\overline{|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{w}_{0}]\rangle \langle \widetilde{w}_{0}|}\}_{\text{H}}+2\text{i}\unicode[STIX]{x1D716}_{\text{nl}}^{2}\{\overline{|f_{\text{nl}}[\widetilde{\unicode[STIX]{x1D719}},\widetilde{w}_{0}]\rangle \langle \widetilde{w}_{0}|}\}_{\text{H}}\nonumber\\ \displaystyle & & \displaystyle +\,\text{i}\unicode[STIX]{x1D716}_{\text{nl}}^{2}\{\overline{|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{w}_{0}]\rangle \langle \widetilde{\unicode[STIX]{x1D719}}|}\}_{\text{H}}+2\text{i}\unicode[STIX]{x1D716}^{2}[\widehat{\mathscr{S}}(\widehat{\mathscr{D}}^{-1})^{\dagger }]_{\text{A}}\nonumber\\ \displaystyle & & \displaystyle +\,O(\unicode[STIX]{x1D716}_{\text{nl}}^{3},\unicode[STIX]{x1D716}_{\text{nl}}^{2}\unicode[STIX]{x1D716},\unicode[STIX]{x1D716}_{\text{nl}}\unicode[STIX]{x1D716}^{2}).\end{eqnarray}$$

Here we neglected $O(\unicode[STIX]{x1D716}_{\text{nl}}^{3},\unicode[STIX]{x1D716}_{\text{nl}}^{2}\unicode[STIX]{x1D716},\unicode[STIX]{x1D716}_{\text{nl}}\unicode[STIX]{x1D716}^{2})$ terms, e.g. $\unicode[STIX]{x1D716}_{\text{nl}}^{3}|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{\unicode[STIX]{x1D719}}]\rangle \langle \widetilde{\unicode[STIX]{x1D719}}|$ and $\unicode[STIX]{x1D716}_{\text{nl}}^{2}\unicode[STIX]{x1D716}|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{\unicode[STIX]{x1D719}}]\rangle \langle \widetilde{\unicode[STIX]{x1D711}}|$ . Note that the zonal averages of quantities involving both $\widetilde{w}_{0}$ and $\widetilde{\unicode[STIX]{x1D709}}$ (e.g. $\unicode[STIX]{x1D716}_{\text{nl}}|\widetilde{\unicode[STIX]{x1D709}}\rangle \langle \widetilde{w}_{0}|$ and $\unicode[STIX]{x1D716}_{\text{nl}}\unicode[STIX]{x1D716}|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{w}_{0}]\rangle \langle \widetilde{\unicode[STIX]{x1D709}}|$ ) are zero since the statistics of $\widetilde{w}_{0}$ and $\widetilde{\unicode[STIX]{x1D709}}$ are independent. The factor of two in the second term on the right-hand side of (3.10) is due to the symmetry property (2.18).

Now, let us explicitly calculate the statistical average of the nonlinear terms in (3.10). To do this, we shall use the quasinormal approximation which expresses higher-order $(n\geqslant 3)$ statistical moments of $\widetilde{w}_{0}$ in terms of the lower-order moments. A further discussion on the validity of this approximation and its relation to the widely used random-phase approximation in homogeneous wave turbulence theory is given in § 5. When adopting the quasinormal approximation, one specifically obtains

(3.11) $$\begin{eqnarray}\displaystyle \overline{\widetilde{w}_{0}(\unicode[STIX]{x1D639}_{1})\widetilde{w}_{0}(\unicode[STIX]{x1D639}_{2})\widetilde{w}_{0}(\unicode[STIX]{x1D639}_{3})}=0, & & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle \overline{\widetilde{w}_{0}(\unicode[STIX]{x1D639}_{1})\widetilde{w}_{0}(\unicode[STIX]{x1D639}_{2})\widetilde{w}_{0}(\unicode[STIX]{x1D639}_{3})\widetilde{w}_{0}(\unicode[STIX]{x1D639}_{4})} & = & \displaystyle \langle \unicode[STIX]{x1D639}_{1}\mid \widehat{\mathscr{W}}_{0}\mid \unicode[STIX]{x1D639}_{2}\rangle \,\langle \unicode[STIX]{x1D639}_{3}\mid \widehat{\mathscr{W}}_{0}\mid \unicode[STIX]{x1D639}_{4}\rangle \nonumber\\ \displaystyle & & \displaystyle +\,\langle \unicode[STIX]{x1D639}_{1}\mid \widehat{\mathscr{W}}_{0}\mid \unicode[STIX]{x1D639}_{3}\rangle \,\langle \unicode[STIX]{x1D639}_{2}\mid \widehat{\mathscr{W}}_{0}\mid \unicode[STIX]{x1D639}_{4}\rangle \nonumber\\ \displaystyle & & \displaystyle +\,\langle \unicode[STIX]{x1D639}_{1}\mid \widehat{\mathscr{W}}_{0}\mid \unicode[STIX]{x1D639}_{4}\rangle \,\langle \unicode[STIX]{x1D639}_{2}\mid \widehat{\mathscr{W}}_{0}\mid \unicode[STIX]{x1D639}_{3}\rangle ,\end{eqnarray}$$

where $\widehat{\mathscr{W}}_{0}\doteq \overline{|\widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|}$ and $\unicode[STIX]{x1D639}_{i}=(t_{i},\boldsymbol{x}_{i})$ denotes a space–time coordinate.

Upon substituting (3.11), one finds that the first term appearing in the right-hand side of (3.10) vanishes; i.e. $\overline{|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{w}_{0}]\rangle \langle \widetilde{w}_{0}|}=0$ . After substituting (2.13) and (3.9a ), one obtains the following expression for the second term on the right-hand side of (3.10):

(3.13) $$\begin{eqnarray}\displaystyle & & \displaystyle \overline{|f_{\text{nl}}[\widetilde{\unicode[STIX]{x1D719}},\widetilde{w}_{0}]\rangle \langle \widetilde{w}_{0}|}=\int \text{d}^{3}\unicode[STIX]{x1D639}\,|\unicode[STIX]{x1D639}\rangle \overline{\langle \widetilde{\unicode[STIX]{x1D719}}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|}=-\frac{\text{i}}{2}\int \text{d}^{3}\unicode[STIX]{x1D639}\,\text{d}^{3}\unicode[STIX]{x1D63A}\,|\unicode[STIX]{x1D639}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad \times \,\overline{[\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle -\overline{\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle }]\langle \unicode[STIX]{x1D63A}\mid (\widehat{\mathscr{D}}^{-1})^{\dagger }\,\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|},\end{eqnarray}$$

where $|\unicode[STIX]{x1D639}\rangle$ and $|\unicode[STIX]{x1D63A}\rangle$ are eigenstates of the coordinate operator. Inserting (3.12) in the abstract representation gives

(3.14) $$\begin{eqnarray}\displaystyle & & \displaystyle \overline{[\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle -\overline{\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle }]\langle \unicode[STIX]{x1D63A}\mid (\widehat{\mathscr{D}}^{-1})^{\dagger }\,\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|}\nonumber\\ \displaystyle & & \displaystyle \quad =\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle \langle \unicode[STIX]{x1D63A}\mid (\widehat{\mathscr{D}}^{-1})^{\dagger }\,\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle \langle \unicode[STIX]{x1D63A}\mid (\widehat{\mathscr{D}}^{-1})^{\dagger }\,\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|\nonumber\\ \displaystyle & & \displaystyle \quad =2\langle \unicode[STIX]{x1D63A}\mid (\widehat{\mathscr{D}}^{-1})^{\dagger }\,\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|\nonumber\\ \displaystyle & & \displaystyle \quad =2\langle \unicode[STIX]{x1D63A}|\,(\widehat{\mathscr{D}}^{-1})^{\dagger }\,\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\,\overline{|\widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|}\,\widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\,\overline{|\widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|}\nonumber\\ \displaystyle & & \displaystyle \quad =2\langle \unicode[STIX]{x1D63A}|\,(\widehat{\mathscr{D}}^{-1})^{\dagger }\,\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\,\widehat{\mathscr{W}}_{0}\,\widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\,\widehat{\mathscr{W}}_{0},\end{eqnarray}$$

where the connecting overlines indicate the correlations of the terms that contribute to the final result. (This notation is also commonly used in quantum field theory (Peskin & Schroeder 1995).) It is to be noted that this result can also be obtained by inserting the completeness relation $\widehat{1}=\int \text{d}^{3}\unicode[STIX]{x1D639}^{\prime }|\unicode[STIX]{x1D639}^{\prime }\rangle \langle \unicode[STIX]{x1D639}^{\prime }|$ and then using (3.12) explicitly. In the third line, we used the symmetry property given in (2.18). Substituting this result into (3.13) gives

(3.15) $$\begin{eqnarray}\displaystyle \overline{|f_{\text{nl}}[\widetilde{\unicode[STIX]{x1D719}},\widetilde{w}_{0}]\rangle \langle \widetilde{w}_{0}|}=-\text{i}\int \text{d}^{3}\unicode[STIX]{x1D639}\,\text{d}^{3}\unicode[STIX]{x1D63A}\,|\unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D63A}|(\widehat{\mathscr{D}}^{-1})^{\dagger }\,\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\,\widehat{\mathscr{W}}_{0}\,\widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\,\widehat{\mathscr{W}}_{0}. & & \displaystyle\end{eqnarray}$$

In a similar manner, substituting (3.9a ) into the third term in (3.10) leads to

(3.16) $$\begin{eqnarray}\displaystyle \overline{|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{w}_{0}]\rangle \langle \widetilde{\unicode[STIX]{x1D719}}|} & = & \displaystyle \int \text{d}^{3}\unicode[STIX]{x1D639}\,|\unicode[STIX]{x1D639}\rangle \overline{\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{\unicode[STIX]{x1D719}}|}=-\frac{\text{i}}{2}\int \text{d}^{3}\unicode[STIX]{x1D639}\,\text{d}^{3}\unicode[STIX]{x1D63A}\,|\unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D63A}|(\widehat{\mathscr{D}}^{-1})^{\dagger }\nonumber\\ \displaystyle & & \displaystyle \times \,\overline{\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle [\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle -\overline{\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle }]}.\end{eqnarray}$$

As before, when using the quasinormal approximation, we obtain

(3.17) $$\begin{eqnarray}\displaystyle & & \displaystyle \overline{\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \,[\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle -\overline{\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle }]}\nonumber\\ \displaystyle & & \displaystyle \quad =\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle +\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =2\langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\mid \widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}\mid \widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\mid \widetilde{w}_{0}\rangle \nonumber\\ \displaystyle & & \displaystyle \quad =2\,\text{Tr}[\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\,\overline{|\widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|}\,\widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\,\overline{|\widetilde{w}_{0}\rangle \langle \widetilde{w}_{0}|}]\nonumber\\ \displaystyle & & \displaystyle \quad =2\,\text{Tr}[\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\,\widehat{\mathscr{W}}_{0}\,\widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\,\widehat{\mathscr{W}}_{0}],\end{eqnarray}$$

where we introduced the trace operation $\text{Tr}\,\widehat{\mathscr{A}}=\int \text{d}^{3}\unicode[STIX]{x1D639}\langle \unicode[STIX]{x1D639}\mid \widehat{\mathscr{A}}\mid \unicode[STIX]{x1D639}\rangle$ . By using the identity operator $\widehat{1}=\int \text{d}^{3}\unicode[STIX]{x1D639}^{\prime }|\unicode[STIX]{x1D639}^{\prime }\rangle \langle \unicode[STIX]{x1D639}^{\prime }|$ , we are able to write $\langle w_{0}\mid \widehat{\mathscr{A}}\mid w_{0}\rangle =\text{Tr}(\widehat{\mathscr{A}}|w_{0}\rangle \langle w_{0}|)$ . Also, note that the trace $\text{Tr}[\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\,\widehat{\mathscr{W}}_{0}\,\widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\,\widehat{\mathscr{W}}_{0}]$ is not an operator but rather a function of the space–time coordinates $\unicode[STIX]{x1D639}$ and $\unicode[STIX]{x1D63A}$ . Upon inserting (3.17) into (3.16), we obtain

(3.18) $$\begin{eqnarray}\displaystyle \overline{|f_{\text{nl}}[\widetilde{w}_{0},\widetilde{w}_{0}]\rangle \langle \widetilde{\unicode[STIX]{x1D719}}|}=-\text{i}\int \text{d}^{3}\unicode[STIX]{x1D639}\,\text{d}^{3}\unicode[STIX]{x1D63A}\,|\unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D63A}|(\widehat{\mathscr{D}}^{-1})^{\dagger }\,\text{Tr}[\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\,\widehat{\mathscr{W}}_{0}\,\widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\,\widehat{\mathscr{W}}_{0}]. & & \displaystyle\end{eqnarray}$$

We then substitute (3.15) and (3.18) into (3.10). Approximating $\widehat{\mathscr{W}}_{0}\simeq \widehat{\mathscr{W}}$ , which is valid to the leading order of accuracy of the theory, gives a closed equation for the correlation operator:

(3.19) $$\begin{eqnarray}\displaystyle [\widehat{\mathscr{D}}_{\text{H}},\widehat{\mathscr{W}}]_{-}+\text{i}[\widehat{\mathscr{D}}_{\text{A}},\widehat{\mathscr{W}}]_{+}=2\text{i}\unicode[STIX]{x1D716}_{\text{nl}}^{2}[\widehat{\mathscr{F}}\,(\widehat{\mathscr{D}}^{-1})^{\dagger }]_{\text{A}}-2\text{i}\unicode[STIX]{x1D716}_{\text{nl}}^{2}[\widehat{\unicode[STIX]{x1D702}}\,\widehat{\mathscr{W}}]_{\text{A}}+2\text{i}\unicode[STIX]{x1D716}^{2}[\widehat{\mathscr{S}}(\widehat{\mathscr{D}}^{-1})^{\dagger }]_{\text{A}}, & & \displaystyle\end{eqnarray}$$

where

(3.20a ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\unicode[STIX]{x1D702}}\doteq -\int \text{d}^{3}\unicode[STIX]{x1D639}\,\text{d}^{3}\unicode[STIX]{x1D63A}\,|\unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D63A}|(\widehat{\mathscr{D}}^{-1})^{\dagger }\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\,\widehat{\mathscr{W}}\,\widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A}), & \displaystyle\end{eqnarray}$$

(3.20b ) $$\begin{eqnarray}\displaystyle & \displaystyle \widehat{\mathscr{F}}\doteq \frac{1}{2}\int \text{d}^{3}\unicode[STIX]{x1D639}\,\text{d}^{3}\unicode[STIX]{x1D63A}\,|\unicode[STIX]{x1D639}\rangle \langle \unicode[STIX]{x1D63A}|\text{Tr}[\widehat{\mathscr{K}}(\unicode[STIX]{x1D639})\,\widehat{\mathscr{W}}\,\widehat{\mathscr{K}}^{\dagger }(\unicode[STIX]{x1D63A})\,\widehat{\mathscr{W}}]. & \displaystyle\end{eqnarray}$$

In summary, equations (3.2), (3.19), and (3.20) form a complete set of equations for the zonal flow $U(t,y)$ and for the correlation operator $\widehat{\mathscr{W}}$ . Multiplying (3.19) by $\langle t,\boldsymbol{x}|$ and $|t^{\prime },\boldsymbol{x}^{\prime }\rangle$ leads to the ‘quasinormal’ equation for the two-point correlation function written in the double-physical coordinate space $(\unicode[STIX]{x1D639},\unicode[STIX]{x1D639}^{\prime })$ . However, we shall instead project these equations into the phase space by using the Weyl transform. Then, by adopting approximations that are consistent with the GO description of eikonal fields, we shall obtain a WKE model describing the vorticity fluctuations. (Readers who are not familiar with the Weyl calculus are encouraged to read § A.1 before continuing further.)

4 Dynamics in the ray phase space

4.1 Wigner–Moyal equation

The Weyl transform is a mapping from operators in a Hilbert space into functions of phase space (Tracy et al. 2014). In this work, the Weyl transform is defined as

(4.1) $$\begin{eqnarray}\displaystyle A(\unicode[STIX]{x1D639},\unicode[STIX]{x1D62C})\doteq \unicode[STIX]{x1D61E}[\widehat{\mathscr{A}}]=\int \text{d}^{3}\unicode[STIX]{x1D634}\,\text{e}^{\text{i}\unicode[STIX]{x1D62C}\cdot \unicode[STIX]{x1D634}/\unicode[STIX]{x1D716}}\langle \unicode[STIX]{x1D639}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D634}\mid \widehat{\mathscr{A}}\mid \unicode[STIX]{x1D639}-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D634}\rangle , & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D62C}\doteq (\unicode[STIX]{x1D714},\boldsymbol{k})$ , $\unicode[STIX]{x1D634}\doteq (\unicode[STIX]{x1D70F},\boldsymbol{s})$ , $\unicode[STIX]{x1D62C}\cdot \unicode[STIX]{x1D634}=\unicode[STIX]{x1D714}\unicode[STIX]{x1D70F}-\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{s}$ and $\text{d}^{3}\unicode[STIX]{x1D634}\doteq \text{d}\unicode[STIX]{x1D70F}\,\text{d}^{2}\boldsymbol{s}$ . Here the integrals span over $\mathbb{R}^{3}$ . The Weyl symbol $A(\unicode[STIX]{x1D639},\unicode[STIX]{x1D62C})$ of an operator $\widehat{\mathscr{A}}$ is a function on the extended six-dimensional phase space $(\unicode[STIX]{x1D639},\unicode[STIX]{x1D62C})=(t,\boldsymbol{x},\unicode[STIX]{x1D714},\boldsymbol{k})$ ; i.e. $A=A(t,\boldsymbol{x},\unicode[STIX]{x1D714},\boldsymbol{k})$ . Physically, $A(\unicode[STIX]{x1D639},\unicode[STIX]{x1D62C})$ can be interpreted as a local Fourier transform of the space–time representation of the operator $\widehat{\mathscr{A}}$ (see, for instance, equation (A 3)).

The Weyl symbol $W(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})$ corresponding to $\widehat{\mathscr{W}}$ is referred to as the Wigner function of the vorticity fluctuations (Wigner 1932). Since $W(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})$ is a Weyl symbol of the zonal-averaged operator (3.1), it does not depend on the $x$ coordinate. Upon following § A.2, one can write $W(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})$ explicitly as

(4.2) $$\begin{eqnarray}\displaystyle W(t,y,\unicode[STIX]{x1D714},\boldsymbol{k}) & = & \displaystyle \int \text{d}^{3}\unicode[STIX]{x1D634}\,\text{e}^{\text{i}\unicode[STIX]{x1D62C}\cdot \unicode[STIX]{x1D634}/\unicode[STIX]{x1D716}}\langle \unicode[STIX]{x1D639}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D634}\mid \widehat{\mathscr{W}}\mid \unicode[STIX]{x1D639}-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D634}\rangle \nonumber\\ \displaystyle & = & \displaystyle \int \text{d}^{3}\unicode[STIX]{x1D634}\,\text{d}x\,\text{e}^{\text{i}\unicode[STIX]{x1D62C}\cdot \unicode[STIX]{x1D634}/\unicode[STIX]{x1D716}}\boldsymbol{\langle }\!\langle \widetilde{w}(\unicode[STIX]{x1D639}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D634})\,\widetilde{w}(\unicode[STIX]{x1D639}-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D634})\rangle \!\boldsymbol{\rangle }/L_{x}.\end{eqnarray}$$

Since $\widetilde{w}$ is real, then $W(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})=W(t,y,-\unicode[STIX]{x1D714},-\boldsymbol{k})$ . Also, $W(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})$ is a real function because $\widehat{\mathscr{W}}$ is Hermitian. Similar arguments apply to the Weyl symbol $S(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})$ corresponding to the operator $\widehat{\mathscr{S}}$ .

Applying the Weyl transform to (3.19) leads to the Wigner–Moyal formulation of DW–ZF dynamics:

(4.3) $$\begin{eqnarray}\displaystyle & & \displaystyle \mathbf{\{}\!\{D_{\text{H}},W\}\!\mathbf{\}}+\text{i}\,\mathbf{[}[D_{\text{A}},W]\mathbf{]}\nonumber\\ \displaystyle & & \displaystyle \quad =2\text{i}\unicode[STIX]{x1D716}_{\text{nl}}^{2}\,\text{Im}\{F\star [D^{-1}]^{\ast }\}-2\text{i}\unicode[STIX]{x1D716}_{\text{nl}}^{2}\,\text{Im}(\unicode[STIX]{x1D702}\star W)+2\text{i}\unicode[STIX]{x1D716}^{2}\,\text{Im}\{S\star [D^{-1}]^{\ast }\}.\end{eqnarray}$$

Here $D(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})$ is the Weyl symbol corresponding to $\widehat{\mathscr{D}}$ . From the properties of the Weyl transform, we observe that $D_{\text{H}}=\text{Re}\,D$ and $D_{\text{A}}=\text{Im}\,D$ are the (real) Weyl symbols corresponding to the operators $\widehat{\mathscr{D}}_{\text{H}}$ and $\widehat{\mathscr{D}}_{\text{A}}$ , respectively. (‘ $\text{Re}$ ’ and ‘ $\text{Im}$ ’ denote the real and imaginary parts, respectively.) These are given by

(4.4a ) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\text{H}}(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})=\unicode[STIX]{x1D714}-k_{x}U+{\textstyle \frac{1}{2}}\mathbf{[}[(\unicode[STIX]{x1D6FD}-U^{\prime \prime }),k_{x}/k_{\text{D}}^{2}]\mathbf{]}, & \displaystyle\end{eqnarray}$$

(4.4b ) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\text{A}}(t,y,\boldsymbol{k})=\unicode[STIX]{x1D716}\unicode[STIX]{x1D707}_{\text{dw}}-{\textstyle \frac{1}{2}}\mathbf{\{}\!\{U^{\prime \prime },k_{x}/k_{\text{D}}^{2}\}\!\mathbf{\}}. & \displaystyle\end{eqnarray}$$

In addition, $F(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})$ , $\unicode[STIX]{x1D702}(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})$ and $[D^{-1}]^{\ast }(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})$ are the Weyl symbols corresponding to $\widehat{\mathscr{F}}$ , $\widehat{\mathscr{\{}}$ and $(\widehat{\mathscr{D}}^{-1})^{\dagger }$ , respectively. (These will be calculated explicitly later.) The Moyal product (A 5) is denoted by ‘ $\star$ ’, and the brackets $\mathbf{\{}\!\{\cdot ,\cdot \}\!\mathbf{\}}$ and $\mathbf{[}[\cdot ,\cdot ]\mathbf{]}$ are the Moyal brackets ((A 9) and (A 10)). Basic properties of the Moyal product and the Moyal brackets are given in § A.1.

Modulo the statistical closure introduced in § 3, equation (4.3) is an exact equation for the dynamics of the Wigner function. However, equation (4.3) is difficult to solve as is, both analytically and numerically. Thus, we shall reduce (4.3) to a first-order partial-differential equation (PDE) in phase space by using the GO approximation.

4.2 Collisional wave kinetic equation

In order to simplify (4.3), we shall expand the Moyal products and brackets in terms of the ordering parameter $\unicode[STIX]{x1D716}$ . As a reminder, the parameter $\unicode[STIX]{x1D716}$ was introduced in (2.7) in order to denote the spatio-temporal scale separation between the DW and ZF dynamics. From the definition of the Moyal product (A 6), one has

(4.5) $$\begin{eqnarray}\displaystyle A\star B=AB+\frac{\text{i}\unicode[STIX]{x1D716}}{2}(A\overleftrightarrow{{\mathcal{L}}}B)+O(\unicode[STIX]{x1D716}^{2}), & & \displaystyle\end{eqnarray}$$

where $\overleftrightarrow{{\mathcal{L}}}$ is the Janus operator (A 7), which basically serves as the canonical Poisson bracket in the extended six-dimensional phase space $(t,\boldsymbol{x},\unicode[STIX]{x1D714},\boldsymbol{k})$ . Likewise, the Moyal brackets in (A 9) and (A 10) are approximated by

(4.6a,b ) $$\begin{eqnarray}\displaystyle \mathbf{\{}\!\{A,B\}\!\mathbf{\}}=\text{i}\unicode[STIX]{x1D716}(A\overleftrightarrow{{\mathcal{L}}}B)+O(\unicode[STIX]{x1D716}^{3}),\quad \mathbf{[}[A,B]\mathbf{]}=2AB+O(\unicode[STIX]{x1D716}^{2}). & & \displaystyle\end{eqnarray}$$

For the purposes of this paper, higher-order corrections to (4.5) and (4.6) will not be needed. It is to be noted that the asymptotic expansions above are valid as long as the functions involved are smooth.

In addition to asymptotically expanding the Moyal products and brackets, let us find a proper ansatz for the Wigner function. Note that (3.3) can be written as $\widehat{\mathscr{D}}_{\text{H}}\widehat{\mathscr{W}}=O(\unicode[STIX]{x1D716},\unicode[STIX]{x1D716}_{\text{nl}})$ , where we assumed small dissipation $[D_{\text{A}}\sim O(\unicode[STIX]{x1D716})]$ . In the Weyl representation, this equation becomes $D_{\text{H}}\star W=O(\unicode[STIX]{x1D716},\unicode[STIX]{x1D716}_{\text{nl}})$ . Using (4.5) gives

(4.7) $$\begin{eqnarray}\displaystyle D_{\text{H}}(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})\,W(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})\simeq O(\unicode[STIX]{x1D716},\unicode[STIX]{x1D716}_{\text{nl}}). & & \displaystyle\end{eqnarray}$$

To satisfy this equation, we adopt the GO ansatz (McDonald & Kaufman 1985; McDonald 1991)

(4.8) $$\begin{eqnarray}\displaystyle W(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})=2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}\,\unicode[STIX]{x1D6FF}\boldsymbol{(}D_{\text{H}}(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})\boldsymbol{)}\,J(t,y,\boldsymbol{k}). & & \displaystyle\end{eqnarray}$$

Here $J(t,y,\boldsymbol{k})$ is interpreted as the wave-action density for the vorticity fluctuations, and $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}$ is added to ensure the proper normalization.

To remain consistent with the asymptotic expansion in (4.5) and (4.6), the symbols $D_{\text{H}}$ and $D_{\text{A}}$ in (4.4) are approximated as well to lowest order in $\unicode[STIX]{x1D716}$ :

(4.9a ) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\text{H}}(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})\simeq \unicode[STIX]{x1D714}-k_{x}U+(\unicode[STIX]{x1D6FD}-U^{\prime \prime })k_{x}/k_{\text{D}}^{2}, & \displaystyle\end{eqnarray}$$

(4.9b ) $$\begin{eqnarray}\displaystyle & \displaystyle D_{\text{A}}(t,y,\boldsymbol{k})\simeq \unicode[STIX]{x1D716}\unicode[STIX]{x1D707}_{\text{dw}}+\unicode[STIX]{x1D716}k_{x}k_{y}U^{\prime \prime \prime }/k_{\text{D}}^{4}, & \displaystyle\end{eqnarray}$$

where we substituted (4.6).

Let us now obtain an approximate expression for $[D^{-1}]^{\ast }$ in (4.3). As a reminder, $[D^{-1}]$ and $D^{-1}$ are not equivalent; $[D^{-1}]$ is the Weyl symbol of $\widehat{\mathscr{D}}^{-1}$ while $D^{-1}$ is simply the inverse of the Weyl symbol $D$ corresponding to $\widehat{\mathscr{D}}$ . However, note that the Weyl representation of $\widehat{\mathscr{D}}\widehat{\mathscr{D}}^{-1}=\widehat{1}$ is $D\star [D^{-1}]=1$ . By replacing the Moyal product with an ordinary product, we obtain $[D^{-1}]\simeq D^{-1}$ . Then, by using the Sokhotski–Plemelj theorem, we replace $[D^{-1}]^{\ast }$ with its limiting form as $\unicode[STIX]{x1D716}$ tends to zero:

(4.10) $$\begin{eqnarray}\displaystyle [D^{-1}]^{\ast }\simeq \lim _{\unicode[STIX]{x1D716}\rightarrow 0}\,\frac{1}{D_{\text{H}}-\text{i}D_{\text{A}}}=\text{i}\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}(D_{\text{H}})+{\mathcal{P}}\,\frac{1}{D_{\text{H}}}, & & \displaystyle\end{eqnarray}$$

where ‘ ${\mathcal{P}}$ ’ denotes the Cauchy principal value.

We then insert (4.6)–(4.10) into (4.3) and integrate over the frequency variable $\unicode[STIX]{x1D714}$ . Afterwards, we expand the Moyal products and brackets using (4.5) and (4.6). Although the GO ansatz (4.8) and the expression for $[D^{-1}]$ in (4.10) involve Dirac delta functions whose derivatives in phase space are not smooth, in § B.1 we show that these singularities can be removed via integration by parts. After some calculations, we obtain the following to lowest order in $\unicode[STIX]{x1D716}$ :

(4.11) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D716}\int \text{d}\unicode[STIX]{x1D714}\,[\unicode[STIX]{x1D716}(D_{\text{H}}\overleftrightarrow{{\mathcal{L}}}J)\,\unicode[STIX]{x1D6FF}(D_{\text{H}})+2D_{\text{A}}\unicode[STIX]{x1D6FF}(D_{\text{H}})J][1+O(\unicode[STIX]{x1D716})]\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D716}_{\text{nl}}^{2}\int \text{d}\unicode[STIX]{x1D714}\,[F\unicode[STIX]{x1D6FF}(D_{\text{H}})-2\text{Im}(\unicode[STIX]{x1D702})J\unicode[STIX]{x1D6FF}(D_{\text{H}})][1+O(\unicode[STIX]{x1D716})]+\unicode[STIX]{x1D716}^{2}\int \text{d}\unicode[STIX]{x1D714}\,S\unicode[STIX]{x1D6FF}(D_{\text{H}})[1+O(\unicode[STIX]{x1D716})].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Note that in (4.11) the Weyl symbols $W$ , $F$ and $S$ are real because their associated operators are Hermitian.

For simplicity, we consider that the nonlinear scaling parameter $\unicode[STIX]{x1D716}_{\text{nl}}$ and the GO parameter $\unicode[STIX]{x1D716}$ scale as $\unicode[STIX]{x1D716}\sim \unicode[STIX]{x1D716}_{\text{nl}}$ . Upon integrating (4.11) over the frequency and neglecting higher-order terms in $\unicode[STIX]{x1D716}$ , we obtain the collisional wave kinetic equation (cWKE)

(4.12) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}J+\{J,\unicode[STIX]{x1D6FA}\}=-2\unicode[STIX]{x1D707}_{\text{dw}}J+2\unicode[STIX]{x1D6E4}J+S_{\text{ext}}+\unicode[STIX]{x1D716}\,C[J,J], & & \displaystyle\end{eqnarray}$$

where $\{\cdot ,\cdot \}=\overleftarrow{\unicode[STIX]{x2202}_{\boldsymbol{x}}}\cdot \overrightarrow{\unicode[STIX]{x2202}_{\boldsymbol{k}}}-\overleftarrow{\unicode[STIX]{x2202}_{\boldsymbol{k}}}\cdot \overrightarrow{\unicode[STIX]{x2202}_{\boldsymbol{x}}}$ . The wave frequency $\unicode[STIX]{x1D6FA}$ and dissipation term $\unicode[STIX]{x1D6E4}$ are respectively obtained from the lowest-order expansion in $\unicode[STIX]{x1D716}$ of $D_{\text{H}}$ and $D_{\text{A}}$ in (4.9):

(4.13a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}(t,y,\boldsymbol{k})\doteq k_{x}U-(\unicode[STIX]{x1D6FD}-U^{\prime \prime })k_{x}/k_{\text{D}}^{2}, & \displaystyle\end{eqnarray}$$

(4.13b ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6E4}(t,y,\boldsymbol{k})\doteq -U^{\prime \prime \prime }k_{x}k_{y}/k_{\text{D}}^{4}. & \displaystyle\end{eqnarray}$$

The external source term $S_{\text{ext}}(t,y,\boldsymbol{k})$ appearing in (4.12) is given by

(4.14) $$\begin{eqnarray}\displaystyle S_{\text{ext}}(t,y,\boldsymbol{k})=\int \text{d}^{3}\unicode[STIX]{x1D634}\,\text{d}x\,\text{e}^{\text{i}\unicode[STIX]{x1D62C}\cdot \unicode[STIX]{x1D634}/\unicode[STIX]{x1D716}}\boldsymbol{\langle }\!\langle \widetilde{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D639}+{\textstyle \frac{1}{2}}\unicode[STIX]{x1D634})\,\widetilde{\unicode[STIX]{x1D709}}(\unicode[STIX]{x1D639}-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D634})\rangle \!\boldsymbol{\rangle }/L_{x}. & & \displaystyle\end{eqnarray}$$

Assuming that $\widetilde{\unicode[STIX]{x1D709}}(t,\boldsymbol{x})$ is white noise leads to

(4.15) $$\begin{eqnarray}\displaystyle \int \text{d}x\,\boldsymbol{\langle }\!\langle \widetilde{\unicode[STIX]{x1D709}}(t,\boldsymbol{x})\widetilde{\unicode[STIX]{x1D709}}(t^{\prime },\boldsymbol{x}^{\prime })\rangle \!\boldsymbol{\rangle }/L_{x}=\unicode[STIX]{x1D6FF}(t-t^{\prime })\,\unicode[STIX]{x1D6EF}\boldsymbol{(}{\textstyle \frac{1}{2}}(y+y^{\prime }),\boldsymbol{x}-\boldsymbol{x}^{\prime }\boldsymbol{)}. & & \displaystyle\end{eqnarray}$$

Thus, the Weyl symbol of the zonal-averaged operator for the stochastic forcing is

(4.16) $$\begin{eqnarray}\displaystyle S_{\text{ext}}(t,y,\boldsymbol{k})=2\int \text{d}^{2}\boldsymbol{s}\,\unicode[STIX]{x1D6EF}(y,\boldsymbol{s})\cos (\boldsymbol{p}\cdot \boldsymbol{s}/\unicode[STIX]{x1D716}). & & \displaystyle\end{eqnarray}$$

Here we used $\unicode[STIX]{x1D6EF}(y,\boldsymbol{s})=\unicode[STIX]{x1D6EF}(y,-\boldsymbol{s})$ , which is due to the fact that $\widetilde{\unicode[STIX]{x1D709}}$ is real by definition.

The term $C[J,J](t,y,\boldsymbol{k})$ in (4.12) represents wave–wave collisions and is given by

(4.17) $$\begin{eqnarray}\displaystyle C[J,J](t,y,\boldsymbol{k})\doteq S_{\text{nl}}[J,J]-2\unicode[STIX]{x1D6FE}_{\text{nl}}[J]J. & & \displaystyle\end{eqnarray}$$

As shown in § B.2, to the leading order in $\unicode[STIX]{x1D716}$ , $\unicode[STIX]{x1D6FE}_{\text{nl}}[J]$ and $S_{\text{nl}}[J,J]$ are given by

(4.18a ) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FE}_{\text{nl}}[J](t,y,\boldsymbol{k})\doteq \int \frac{\text{d}^{2}\boldsymbol{p}\,\text{d}^{2}\boldsymbol{q}}{(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716})^{2}}\,\unicode[STIX]{x1D6FF}^{2}(\boldsymbol{k}-\boldsymbol{p}-\boldsymbol{q})\,\unicode[STIX]{x1D6E9}(t,y,\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})\,M(\boldsymbol{p},\boldsymbol{q})M(\boldsymbol{p},\boldsymbol{k})\,J(t,y,\boldsymbol{p}), & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

(4.18b ) $$\begin{eqnarray}\displaystyle & \displaystyle S_{\text{nl}}[J,J](t,y,\boldsymbol{k})\doteq \int \frac{\text{d}^{2}\boldsymbol{p}\,\text{d}^{2}\boldsymbol{q}}{(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716})^{2}}\,\unicode[STIX]{x1D6FF}^{2}(\boldsymbol{k}-\boldsymbol{p}-\boldsymbol{q})\,\unicode[STIX]{x1D6E9}(t,y,\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})\,|M(\boldsymbol{p},\boldsymbol{q})|^{2}\,J(t,y,\boldsymbol{p})J(t,y,\boldsymbol{q}). & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Here $\unicode[STIX]{x1D6E9}(t,y,\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})\doteq \unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}(\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA})$ , and

(4.19) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}(t,y,\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})\doteq \unicode[STIX]{x1D6FA}(t,y,\boldsymbol{k})-\unicode[STIX]{x1D6FA}(t,y,\boldsymbol{p})-\unicode[STIX]{x1D6FA}(t,y,\boldsymbol{q}). & & \displaystyle\end{eqnarray}$$

One can identify $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FA}=0$ as the frequency-resonance condition. Finally, the kernel $M(\boldsymbol{p},\boldsymbol{q})$ in (4.18) is

(4.20) $$\begin{eqnarray}\displaystyle M(\boldsymbol{p},\boldsymbol{q})\doteq \boldsymbol{e}_{z}\cdot (\boldsymbol{p}\times \boldsymbol{q})(q_{\text{D}}^{-2}-p_{\text{D}}^{-2}). & & \displaystyle\end{eqnarray}$$

4.3 Dynamics of the ZF velocity $U$

Returning to (2.19) for the zonal flow velocity, the term on the right-hand side can be rewritten in terms of the Wigner function $W$ by using (A 4). One has

(4.21) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}U+\unicode[STIX]{x1D707}_{\text{zf}}U & = & \displaystyle \unicode[STIX]{x1D716}_{\text{nl}}^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\int \frac{\text{d}\unicode[STIX]{x1D714}\,\text{d}^{2}\boldsymbol{k}}{(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716})^{3}}\,\frac{k_{x}}{k_{\text{D}}^{2}}\star W(t,y,\unicode[STIX]{x1D714},\boldsymbol{k})\star \frac{k_{y}}{k_{\text{D}}^{2}},\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D716}_{\text{nl}}^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\int \frac{\text{d}^{2}\boldsymbol{k}}{(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716})^{2}}\,\frac{k_{x}}{k_{\text{D}}^{2}}\star J(t,y,\boldsymbol{k})\star \frac{k_{y}}{k_{\text{D}}^{2}},\end{eqnarray}$$

where we used the Moyal product and substituted (4.8). Upon following Ruiz et al. (2016), we substitute $\unicode[STIX]{x1D716}\sim \unicode[STIX]{x1D716}_{\text{nl}}$ and obtain to lowest order in $\unicode[STIX]{x1D716}$ :

(4.22) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x2202}_{t}U+\unicode[STIX]{x1D707}_{\text{zf}}U=\unicode[STIX]{x1D716}^{2}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}y}\int \frac{\text{d}^{2}\boldsymbol{k}}{(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716})^{2}}\,\frac{k_{x}k_{y}}{k_{\text{D}}^{4}}J(t,y,\boldsymbol{k}). & & \displaystyle\end{eqnarray}$$

5 Discussion

5.1 Main equations

Equations (4.12) and (4.22), as well as (4.13), (4.17) and (4.18), are the main result of our work. These equations describe the coupled interaction between an incoherent wave bath of DWs and a coherent ZF velocity field. The cWKE (4.12) governs the dynamics of the wave-action density $J$ for DWs. The left-hand side of (4.12) describes the wave refraction governed by the wave frequency $\unicode[STIX]{x1D6FA}$ (4.13a ), which serves as a Hamiltonian for the system. On the right-hand side, $\unicode[STIX]{x1D707}_{\text{dw}}$ represents weak dissipation due to the external environment, and $\unicode[STIX]{x1D6E4}$ denotes linear dissipation caused by the ZFs (Parker 2018; Ruiz et al. 2016). The term $S_{\text{ext}}$ represents an external source term for the DW fluctuations.

The nonlinear term $C[J,J]$ in (4.12) plays the role of a wave scattering operator. It is composed of two terms, $\unicode[STIX]{x1D6FE}_{\text{nl}}$ and $S_{\text{nl}}$ , which arise from nonlinear wave–wave interactions. The nonlinear source term $S_{\text{nl}}$ in (4.18b ) is a bilinear functional on the action density $J$ . It is always positive and represents contributions to $J$ coming from waves with wavevectors $\boldsymbol{p}$ and $\boldsymbol{q}$ different from $\boldsymbol{k}$ . This term is also known as (the variance of) incoherent noise (Krommes 2002). The nonlinear damping-rate term $\unicode[STIX]{x1D6FE}_{\text{nl}}$ in (4.18a ) linearly depends on $J$ and represents a sink term where the wave action in the $\boldsymbol{k}$ wavevector is transferred to other modes with different wavevectors. The effects described by $\unicode[STIX]{x1D6FE}_{\text{nl}}$ are called the coherent response (Krommes 2002). The other terms in (4.12) and (4.22) are the same as in the QL theory, and their physical meaning was already discussed elsewhere (Parker 2016; Ruiz et al. 2016; Parker 2018; Zhu et al. 2018a , b , c ).

5.2 Conservation properties

Equations (4.12) and (4.22) inherit the same conservation laws of the original gHME (2.1). In other words, for isolated systems ( $S_{\text{ext}}=0$ and $\unicode[STIX]{x1D707}_{\text{zf},\text{dw}}=0$ ), equations (4.12) and (4.22) conserve the total enstrophy and total energy

(5.1a,b ) $$\begin{eqnarray}\displaystyle {\mathcal{Z}}={\mathcal{Z}}_{\text{zf}}+\unicode[STIX]{x1D716}^{2}{\mathcal{Z}}_{\text{dw}},\quad {\mathcal{E}}={\mathcal{E}}_{\text{zf}}+\unicode[STIX]{x1D716}^{2}{\mathcal{E}}_{\text{dw}}, & & \displaystyle\end{eqnarray}$$

where we used $\unicode[STIX]{x1D716}\sim \unicode[STIX]{x1D716}_{\text{nl}}$ . Also, the expressions for the DW and ZF components of the enstrophy and energy are

(5.2a ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{Z}}_{\text{zf}}\doteq \frac{1}{2}\int \text{d}y\,\overline{w}^{2}=\frac{1}{2}\int \text{d}y\,(U^{\prime })^{2}, & \displaystyle\end{eqnarray}$$

(5.2b ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{Z}}_{\text{dw}}\doteq \frac{1}{2}\int \text{d}^{2}\boldsymbol{x}\,\boldsymbol{\langle }\!\langle \widetilde{w}^{2}\rangle \!\boldsymbol{\rangle }=\frac{1}{2}\int \frac{\text{d}y\,\text{d}^{2}\boldsymbol{k}}{(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716})^{2}}\,J, & \displaystyle\end{eqnarray}$$

(5.2c ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{E}}_{\text{zf}}\doteq -\frac{1}{2}\int \text{d}y\,\overline{w}\overline{\unicode[STIX]{x1D713}}=\frac{1}{2}\int \text{d}y\,U^{2}, & \displaystyle\end{eqnarray}$$

(5.2d ) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{E}}_{\text{dw}}\doteq -\frac{1}{2}\int \text{d}^{2}\boldsymbol{x}\,\boldsymbol{\langle }\!\langle \widetilde{w}\widetilde{\unicode[STIX]{x1D713}}\rangle \!\boldsymbol{\rangle }=\frac{1}{2}\int \frac{\text{d}y\,\text{d}^{2}\boldsymbol{k}}{(2\unicode[STIX]{x03C0}\unicode[STIX]{x1D716})^{2}}\,\frac{J}{k_{\text{D}}^{2}}. & \displaystyle\end{eqnarray}$$

These conservation laws are proven in § B.3.

5.3 Comparison with quasilinear models and weak turbulence theory

There exists a vast literature of QL WKE-based models of DW–ZF interactions. These models, which are called ‘improved WKE’ (iWKE) by Zhu et al. (2018a , b ) and ‘CE2-GO’ by Parker (2016, 2018), neglect wave–wave collisions and only differ from our derived model by letting $C[J,J]=0$ . Earlier works also reported a simpler QL WKE (Diamond et al. 1994; Smolyakov & Diamond 1999; Smolyakov et al. 2000; Malkov & Diamond 2001; Malkov et al. 2001; Kaw et al. 2002; Kim & Diamond 2003; Diamond et al. 2005; Trines et al. 2005; Singh et al. 2014) with $\unicode[STIX]{x1D6FA}=k_{x}U-\unicode[STIX]{x1D6FD}k_{y}/k_{\text{D}}^{2}$ and $\unicode[STIX]{x1D6E4}=0$ . This simpler WKE-based model is referred as the ‘traditional WKE’ (tWKE) in Ruiz et al. (2016), Zhu et al. (2018c , a , b ) and Parker (2018). It is also referred as ‘WKE’ in Parker (2016). It was shown in Parker (2016) that the tWKE does not conserve the total enstrophy (in contradiction with the underlying gHME model) and also leads to unphysical growth rates of the zonostrophic instability. In contrast, the iWKE conserves total enstrophy and total energy and shows good qualitative agreement with direct numerical simulations of the QL gHME at least in certain regimes (Parker 2016, 2018). Also, the iWKE was useful to obtain clarifying insights into the structure of the DW phase space and its relation to the Rayleigh–Kuo parameter $U^{\prime \prime }/\unicode[STIX]{x1D6FD}$ (Zhu et al. 2018b ). For a more detailed discussion on the iWKE, see Parker (2018) and Zhu et al. (2018c ).

The WKE is based on the assumption of a scale separation between the DW and ZF dynamics. There also exist QL statistical models that do not assume scale separation. One notable example is the second-order quasilinear expansion, or CE2 (Farrell & Ioannou 2003, 2007; Marston et al. 2008; Srinivasan & Young 2012; Ait-Chaalal et al. 2016), whose applications to DW–ZF physics were pursued in Farrell & Ioannou (2009), Parker & Krommes (2013, 2014) and Parker (2014). Another QL theory is the Wigner–Moyal (WM) formulation (Ruiz et al. 2016; Zhu et al. 2018c ) which describes the DW dynamics in ray phase space. Both the CE2 and WM theories are mathematically equivalent. However, the CE2 is based on the double-physical space representation (§ 3.2) so it is not very intuitive. For the same reason, its robustness with respect to further approximations remains obscure. The WM formulation involves an infinite-order PDE in phase space, which arguably makes it more tractable.

It is worth mentioning that DW scattering has been widely studied in the context of homogeneous weak turbulence theory (WTT). (For recent reviews, see, e.g. Krommes (2002), Nazarenko (2011) and Connaughton et al. (2015).) Previous works in WTT have derived wave–wave collision operators for DWs. However, in these models, both the DW and ZF components of the fields were considered as incoherent. In other words, the quasinormal approximation, or the related random phase approximation, where applied to both the mean and fluctuating quantities. However, direct statistical simulations based the QL approximation have shown that ZFs are not incoherent (Farrell & Ioannou 2009; Srinivasan & Young 2012; Tobias & Marston 2013; Parker & Krommes 2013, 2014). In light of this, equations (4.12) and (4.22) make a distinction between the statistics of the DWs and the ZF velocity field: the DW component is modelled as an incoherent wave ensemble (or wave bath) described by the cWKE, while ZFs are treated as coherent structures. Such partitioning of the statistics is expected to lead to different results.

In the absence of ZFs ( $U=0$ ), the collision operator $C[J,J]$ in (4.17) coincides with those derived previously within homogeneous WTT (see, e.g. Connaughton et al. (2015)). Hence, the cWKE (4.12) is expected to yield the same stationary, homogeneous, non-equilibrium solutions found in WTT. In other words, one would obtain the well-known Kolmogorov–Zakharov spectra for DW turbulence. In this work, we shall not further discuss this topic. The interested reader can refer to, e.g. § 3.3 of Connaughton et al. (2015).

Upon using the wave–wave collision operator that we propose here, future work will be devoted to studying the possible modifications to the DW–ZF dynamics caused by $C[J,J]$ at non-zero $U$ . In particular, the frequency-resonance condition in WTT depends on a simplified expression for the DW frequency: $\unicode[STIX]{x1D6FA}(t,y,\boldsymbol{k})\simeq -\unicode[STIX]{x1D6FD}k_{x}/k_{\text{D}}^{2}$ , which neglects the ZF velocity. Hence, one has

$$\begin{eqnarray}\unicode[STIX]{x1D6E9}_{\text{WTT}}(t,y,\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})=\frac{\unicode[STIX]{x03C0}}{\unicode[STIX]{x1D6FD}}\,\unicode[STIX]{x1D6FF}\left(\frac{k_{x}}{k_{\text{D}}^{2}}-\frac{p_{x}}{p_{\text{D}}^{2}}-\frac{q_{x}}{q_{\text{D}}^{2}}\right).\end{eqnarray}$$

In contrast, in the cWKE, $\unicode[STIX]{x1D6FA}$ is given by (4.13a ), so

(5.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E9}(t,y,\boldsymbol{k},\boldsymbol{p},\boldsymbol{q})=\frac{\unicode[STIX]{x03C0}}{|\unicode[STIX]{x1D6FD}-U^{\prime \prime }(t,y)|}\,\unicode[STIX]{x1D6FF}\left(\frac{k_{x}}{k_{\text{D}}^{2}}-\frac{p_{x}}{p_{\text{D}}^{2}}-\frac{q_{x}}{q_{\text{D}}^{2}}\right), & & \displaystyle\end{eqnarray}$$

where the contribution of the Doppler shifts cancelled out due to the spatial-resonance condition $(\boldsymbol{k}=\boldsymbol{p}+\boldsymbol{q})$ . The factor $|\unicode[STIX]{x1D6FD}-U^{\prime \prime }|$ appearing in the denominator is related to the Rayleigh–Kuo threshold. The relevance of this threshold to the GO dynamics of DWs was shown by Zhu et al. (2018b ), and in a broader context, it also marks the onset of the tertiary instability (Kuo 1949; Numata, Ball & Dewar 2007; Zhu et al. 2018c ). From this result, it seems that the present theory breaks down in regions where $\unicode[STIX]{x1D6FD}-U^{\prime \prime }\simeq \unicode[STIX]{x1D716}$ . In such regions, the DW wave frequency tends to zero, and the interaction time between the waves becomes long. In general, WTT is based on the assumption that the interaction time between the waves are small, so this violates the assumed orderings used to close the equations and leads to the breakdown of the model.

5.4 Comparison with the generalized quasilinear approximation

Let us also mention that the present work is somewhat similar in spirit to the generalized quasilinear (GQL) approximation proposed by Marston, Chini & Tobias (2016), where the dynamical fields are separated into small and large zonal scales via a spectral filter that depends on a wavenumber cutoff. In order to incorporate the nonlinear energy transfer due to eddy–eddy interactions, Marston et al. (2016) includes nonlinear large-scale ZF–ZF interactions but neglects small-scale DW–DW interactions. In simulations of zonal-jet formation (Marston et al. 2016) and rotating three-dimensional Couette flow (Tobias & Marston 2016), it was shown that the GQL gives an improvement in accuracy on calculations of the mean flows, spectra and two-point correlation functions over models that are QL.

However, the GQL framework is not the same as the one presented here, namely, due to the following. First, in our work we decompose the fields into mean (i.e. independent of $x$ ) and fluctuating quantities. Within the gHME, nonlinear large-scale ZF–ZF interactions do not appear when using such decomposition. Second, we allow for small-scale DW–DW interactions. In consequence, the governing equation for the fluctuations remains nonlinear. Studying the differences between these models and developing a new framework that takes advantages of each model could be a topic of future research.

5.5 Realizability of the cWKE model

Since it is known that the quasinormal approximation can lead to unrealizable statistics for the Euler equations (Ogura 1963; Leslie 1973), one may wonder if such issues regarding realizability could present themselves in the cWKE model proposed here. We anticipate that this is not the case for the following two reasons. First, the turbulence that we address is of different nature. Unlike the isotropic turbulence governed by the Euler equations, which is dominated by eddy–eddy interactions, the turbulence considered here is mainly determined by collective wave–wave interactions. Second, as mentioned in § 5.3, the cWKE model is identical to homogeneous WTT theory (Krommes 2002; Nazarenko 2011; Connaughton et al. 2015) when no ZFs are present ( $U=0$ ). In homogeneous WTT, both the quasinormal approximation and the random-phase approximation lead to the same wave–wave scattering operators and show realizable statistics. Since the structures of the equations and of the collisional operator in our cWKE model are the same as those found in homogeneous WTT, we believe that the statistical cWKE model proposed here leads to meaningful statistics. However, this remains to be checked.

6 Conclusions

In this work, we present a nonlinear wave kinetic equation for studying DW–ZF interactions. In contrast with previous works that used the quasilinear approximation, here we perturbatively include nonlinear wave–wave collisions. Our derivation makes use of the Weyl calculus, in conjunction with the quasinormal statistical closure and the well-known geometrical-optics assumptions. The obtained model is similar to previous reported works (Ruiz et al. 2016; Parker 2016, 2018; Zhu et al. 2018a , b ) but also includes a nonlinear term describing wave–wave scattering. Unlike in WTT, the collision operator depends on the local ZF velocity and breaks down at the Rayleigh–Kuo threshold. Our model conserves both the total enstrophy and the total energy of the system.

The model presented here might allow us to investigate the effects of wave–wave scattering on the DW–ZF system. Several questions that could be addressed are the following. How will the spontaneous emergence of ZFs (also called the zonostrophic instability) be modified in the presence of wave–wave scattering? Regarding the saturated states, how will the partition of enstrophy and energy between DWs and ZFs change? Will the stationary Kolmogorov–Zakharov spectra for DWs be modified in the presence of ZFs? These and other questions will be subject to future investigations.