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

TL;DR

This paper derives a wave kinetic equation for inhomogeneous drift-wave turbulence that includes wave-wave collisions, extending beyond the quasilinear approximation to better understand drift wave and zonal flow interactions.

Contribution

It introduces a new wave kinetic model with a collision operator that accounts for wave interactions, surpassing the limitations of previous quasilinear theories.

Findings

The model conserves enstrophy and energy.

The collision operator breaks down at the Rayleigh--Kuo threshold.

It highlights the limitations of homogeneous-turbulence theory.

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 literature.

Figures

Equations

$${\partial }_{t}J+\{J,\Omega \}=2\Gamma J,$$

$${\partial }_{t}w+\mathbf{v}\cdot \nabla w+\beta {\partial }_x\psi =Q.$$

$$w(t,\mathbf{x})\doteq ({\nabla }^2-L_{\mathrm{D}}^{-2}\hat{a})\psi ,$$

$$\mathcal{Z}(t)\doteq \frac{1}{2}\int {\mathrm{d}}^2\mathbf{x}{w}^2,\mathcal{E}(t)\doteq -\frac{1}{2}\int {\mathrm{d}}^2\mathbf{x}w\psi .$$

$$w=\overline{w}+{ϵ}_{\mathrm{nl}}\tilde{w},\mathbf{v}=U{\mathbf{e}}_x+{ϵ}_{\mathrm{nl}}\tilde{\mathbf{v}},$$

$$\overline{w}={\nabla }^2\overline{\psi },\tilde{w}={\nabla }_{\mathrm{D}}^2\tilde{\psi }.$$

$${\partial }_t\tilde{w}+U{\partial }_x\tilde{w}+[\beta -({\partial }_y^{2}U)]{\partial }_x\tilde{\psi }+{ϵ}_{\mathrm{nl}}{f}_{\mathrm{eddy}}=\tilde{Q},$$

$${\partial }_{t}U+{\mu }_{\mathrm{zf}}U+{ϵ}_{\mathrm{nl}}^2{\partial }_y\overline{{\tilde{v}}_x{\tilde{v}}_y}=\overline{Q},$$

$$ϵ\doteq \max \left(\frac{{\tau }_{\mathrm{dw}}}{T_{\mathrm{zf}}},\frac{{\lambda }_{\mathrm{dw}}}{L_{\mathrm{zf}}},\frac{L_{\mathrm{D}}}{L_{\mathrm{zf}}}\right)\ll 1.$$

$$ϵ{\partial }_t\tilde{w}+ϵU{\partial }_x\tilde{w}+ϵ[\beta -({\partial }_y^{2}U)]{\partial }_x\tilde{\psi }+ϵ{\mu }_{\mathrm{dw}}\tilde{w}+{ϵ}_{\mathrm{nl}}{f}_{\mathrm{eddy}}=ϵ\tilde{\xi },$$

$${\partial }_{t}U+{\mu }_{\mathrm{zf}}U+{ϵ}_{\mathrm{nl}}^2{\partial }_y\overline{{\tilde{v}}_x{\tilde{v}}_y}=0,$$

$$⟨\varphi \mid \psi ⟩=\int \mathrm{d}t{\mathrm{d}}^2\mathbf{x}{\varphi }^{\ast }(t,\mathbf{x})\psi (t,\mathbf{x}),$$

$${\hat{k}}_{\mathrm{D}}^2\doteq {\hat{k}}^2+1,{\hat{k}}^2\doteq \hat{\mathbf{k}}\cdot \hat{\mathbf{k}}.$$

$$\hat{\mathcal{D}}|\tilde{w}⟩=\frac{\mathrm{i}{ϵ}_{\mathrm{nl}}}{2}|f_{\mathrm{nl}}[\tilde{w},\tilde{w}]⟩-\frac{\mathrm{i}{ϵ}_{\mathrm{nl}}}{2}|\overline{f_{\mathrm{nl}}[\tilde{w},\tilde{w}]}⟩+\mathrm{i}ϵ|\tilde{\xi }⟩.$$

$$\hat{\mathcal{D}}\doteq \hat{\omega }-\hat{U}{\hat{k}}_x+(\beta -{\hat{U}}^{\prime \prime }){\hat{k}}_x{\hat{k}}_{\mathrm{D}}^{-2}+\mathrm{i}ϵ{\mu }_{\mathrm{dw}},$$

$$|f_{\mathrm{nl}}[\varphi ,\psi ]⟩$$

$$|\overline{f_{\mathrm{nl}}[\varphi ,\psi ]}⟩$$

$$\hat{\mathcal{K}}(\mathsf{x})\doteq {\hat{\mathcal{L}}}_j|\mathsf{x}⟩⟨\mathsf{x}|{\hat{\mathcal{R}}}_j+{\hat{\mathcal{R}}}_j|\mathsf{x}⟩⟨\mathsf{x}|{\hat{\mathcal{L}}}_j.$$

$${\hat{\mathcal{L}}}_j\doteq ({\mathbf{e}}_z\times \hat{\mathbf{k}}{)}_j{k}_{\mathrm{D}}^{-2},{\hat{\mathcal{R}}}_j\doteq {\hat{\mathbf{k}}}_j.$$

$$⟨\mathsf{x}\mid f_{\mathrm{nl}}[\tilde{w},\tilde{w}]⟩$$

$$=⟨\tilde{w}\mid {\hat{k}}_{\mathrm{D}}^{-2}({\mathbf{e}}_z\times \hat{\mathbf{k}}{)}_j\mid \mathsf{x}⟩⟨\mathsf{x}\mid {\hat{\mathbf{k}}}_j\mid \tilde{w}⟩+⟨\tilde{w}\mid {\hat{\mathbf{k}}}_j\mid \mathsf{x}⟩⟨\mathsf{x}\mid ({\mathbf{e}}_z\times \hat{\mathbf{k}}{)}_j{\tilde{k}}_{\mathrm{D}}^{-2}\mid \tilde{w}⟩$$

$$=-(⟨\mathsf{x}\mid ({\mathbf{e}}_z\times \hat{\mathbf{k}}{)}_j\mid \tilde{\psi }⟩{)}^{\ast }⟨\mathsf{x}\mid {\hat{\mathbf{k}}}_j\mid \tilde{w}⟩-(⟨\mathsf{x}\mid {\hat{\mathbf{k}}}_j\mid \tilde{w}⟩{)}^{\ast }⟨\mathsf{x}\mid ({\mathbf{e}}_z\times \hat{\mathbf{k}}{)}_j\mid \tilde{\psi }⟩$$

$$=-2({\mathbf{e}}_z\times ϵ\nabla \tilde{\psi })\cdot ϵ\nabla \tilde{w}$$

$$=-2\tilde{\mathbf{v}}\times ϵ\nabla \tilde{w},$$

$$|f_{\mathrm{nl}}[\varphi ,\psi ]⟩=|f_{\mathrm{nl}}[\psi ,\varphi ]⟩$$

$${\partial }_{t}U+{\mu }_{\mathrm{zf}}U={ϵ}_{\mathrm{nl}}^2{\partial }_y\overline{⟨\mathsf{x}\mid {\hat{k}}_x{\hat{k}}_{\mathrm{D}}^{-2}\mid \tilde{w}⟩⟨\tilde{w}\mid {\hat{k}}_{\mathrm{D}}^{-2}{\hat{k}}_y\mid \mathsf{x}⟩}.$$

$$\hat{\mathcal{W}}\doteq \overline{|\tilde{w}⟩⟨\tilde{w}|},$$

$${\partial }_{t}U+{\mu }_{\mathrm{zf}}U={ϵ}_{\mathrm{nl}}^2{\partial }_y⟨\mathsf{x}\mid {\hat{k}}_x{\hat{k}}_{\mathrm{D}}^{-2}\hat{\mathcal{W}}{\hat{k}}_{\mathrm{D}}^{-2}{\hat{k}}_y\mid \mathsf{x}⟩.$$

$$\hat{\mathcal{D}}\hat{\mathcal{W}}=$$

$$\displaystyle[\widehat{\mathscr{D}}_{\rm H},\widehat{\mathscr{W}}]_{-}+{\rm i}[\widehat{\mathscr{D}}_{\rm A},\widehat{\mathscr{W}}]_{+}={\rm i}\epsilon_{\rm nl}\,\big{[}\overline{\ket{\,f_{\rm nl}[\widetilde{w},\widetilde{w}]\,}\bra{\,\widetilde{w}\,}}\big{]}_{\rm H}+2{\rm i}\epsilon\,\big{[}\overline{\ket{\,\widetilde{\xi}\,}\bra{\,\widetilde{w}\,}}\big{]}_{\rm H},$$

$${\hat{\mathcal{D}}}_{\mathrm{H}}$$

$${\hat{\mathcal{D}}}_{\mathrm{A}}$$

$$\hat{\mathcal{D}}|\tilde{w}⟩\simeq \mathrm{i}ϵ|\tilde{\xi }⟩.$$

$$[\widehat{\mathscr{D}}_{\rm H},\widehat{\mathscr{W}}]_{-}+{\rm i}[\widehat{\mathscr{D}}_{\rm A},\widehat{\mathscr{W}}]_{+}=2{\rm i}\epsilon^{2}\,\big{[}\widehat{\mathscr{S}}(\widehat{\mathscr{D}}^{-1})^{\dagger}\big{]}_{\rm A},$$

Full text

1 Introduction

The interaction between drift-wave (DW) turbulence and zonal flows (ZFs) has been widely studied in plasma physics (Diamond et al., 2005; Fujisawa, 2009; Lin, 1998; Biglari et al., 1990; Dorland et al., 2000; Jenko et al., 2000; Connaughton et al., 2015). In the context of magnetic fusion (Fujisawa, 2009; Hillesheim et al., 2016; Conway et al., 2005), the spontaneous emergence of ZFs significantly affects the transport of energy, momentum, and particles. Understanding this phenomenon is critical to improving plasma confinement, but modeling 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, \egthe formation of ZFs (Parker, 2016; Ruiz et al., 2016; Zhu et al., 2018c; Parker, 2018; Zhu et al., 2018a, b; Diamond et al., 2005; Smolyakov & Diamond, 1999; Smolyakov et al., 2000; Malkov & Diamond, 2001; Malkov et al., 2001; Diamond et al., 1994; Kim & Diamond, 2003; Kaw et al., 2002; Trines et al., 2005; Singh et al., 2014). 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):

[TABLE]

where $\Omega$ is the local DW frequency, $\Gamma$ 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, Eq. (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) (Krommes & Kim, 2000; Smolyakov & Diamond, 1999), 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 Eqs. (57) and (67). In this final result, DWs are modeled in a similar manner as in Eq. (1). The difference is that Eq. (57) 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 (Connaughton et al., 2015; Krommes, 2002). 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 (Krommes & Kim, 2000; Smolyakov & Diamond, 1999), which in normalized units is given by

[TABLE]

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 $\beta$. In Eq. (2), time is measured in units of the inverse ion cyclotron frequency while length is measured in units of the ion sound radius. [See, \egZhu et al. (2018c) for more details.] The function $\psi(t,{\boldsymbol{x}})$ is the electric potential, ${\boldsymbol{v}}\doteq{\boldsymbol{e}}_{z}\times{\boldsymbol{\nabla}}\psi$ 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

[TABLE]

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_{\rm D}$ is the ion sound radius. ($L_{\rm D}=1$ in normalized units.) The term $Q(t,{\boldsymbol{x}})$ in Eq. (2) represents external forces and dissipation. Notably, Eq. (2) with $L_{\rm D}\to\infty$ also describes Rossby turbulence in planetary atmospheres (Farrell & Ioannou, 2003, 2007; Marston et al., 2008; Srinivasan & Young, 2012; Ait-Chaalal et al., 2016).

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

[TABLE]

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 $\psi$ 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\mathrm{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 $\epsilon_{\rm nl}\ll 1$. In particular,

[TABLE]

where $U(t,y)\doteq-\partial_{y}\overline{\psi}$ is the ZF velocity field and the two components of the generalized vorticity are related to $\psi$ as (Parker & Krommes, 2013)

[TABLE]

Here $\nabla_{\rm D}^{2}\doteq\nabla^{2}-1$. Also, $\overline{\widetilde{w}}=0$ and $\overline{\widetilde{{\boldsymbol{v}}}}=0$ by definition.

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

[TABLE]

where $\smash{f_{\rm eddy}(t,{\boldsymbol{x}})\doteq\widetilde{{\boldsymbol{v}}}\cdot{\boldsymbol{\nabla}}\widetilde{w}-\overline{\widetilde{{\boldsymbol{v}}}\cdot{\boldsymbol{\nabla}}\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 Eqs. (7), the parameter $\epsilon_{\rm nl}$ denotes the smallness of the nonlinear coupling terms. Also note that we considered the ordering $\smash{Q=\overline{Q}+\epsilon_{\rm 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 $\tau_{\rm dw}$ and $\lambda_{\rm 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_{\rm zf}$ and $L_{\rm 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

[TABLE]

In addition, we assume that the DWs are weakly damped and weakly dissipated so that $\smash{\widetilde{Q}=\epsilon\,\widetilde{\xi}-\epsilon\,\mu_{\rm dw}\widetilde{w}}$, where $\widetilde{\xi}$ is some white-noise external forcing with zero mean and $\mu_{\rm 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 $\smash{\overline{Q}=-\mu_{\rm zf}\,\overline{w}}$, where $\mu_{\rm zf}$ represents dissipation on the mean flows.

Upon using the orderings above, we write Eqs. (7) as

[TABLE]

where the GO parameter $\epsilon$ was inserted in front of the partial derivatives acting on fluctuating quantities in order to denote the small-scale ordering. (Here $\epsilon$ serves the same role as $\hbar$ in quantum mechanics.) Note that $\smash{\widetilde{{\boldsymbol{v}}}={\boldsymbol{e}}_{z}\times\epsilon{\boldsymbol{\nabla}}\widetilde{\psi}}$, $\smash{\widetilde{w}=(\epsilon^{2}\nabla^{2}-1)\widetilde{\psi}}$, and $f_{\rm eddy}=\smash{\widetilde{{\boldsymbol{v}}}\cdot\epsilon{\boldsymbol{\nabla}}\widetilde{w}-\overline{\widetilde{{\boldsymbol{v}}}\cdot\epsilon{\boldsymbol{\nabla}}\widetilde{w}}}$. Also, note $\smash{\partial_{y}\overline{\widetilde{v}_{x}\widetilde{v}_{y}}}\sim O(\epsilon^{0})$ since $\smash{\overline{\widetilde{v}_{x}\widetilde{v}_{y}}}$ is a mean quantity. Equations (9) 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)

[TABLE]

where the integrals are taken over $\mathbb{R}^{3}$. Let $\ket{\,t,{\boldsymbol{x}}\,}$ 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}})=\braket{\,t,{\boldsymbol{x}}\mid\widetilde{w}\,}$. Since $\widetilde{w}(t,{\boldsymbol{x}})$ is real, then $\braket{\,\widetilde{w}\mid t,{\boldsymbol{x}}\,}=\braket{\,t,{\boldsymbol{x}}\mid\widetilde{w}\,}$. In the following, we shall use the notation $\mathsf{x}\doteq(t,{\boldsymbol{x}})$ to denote spacetime coordinates. Likewise, $\ket{\,\mathsf{x}\,}\doteq\ket{\,t,{\boldsymbol{x}}\,}$ is the corresponding space-time eigenstate, and $\mathrm{d}^{3}\mathsf{x}\doteq\mathrm{d}t\,\mathrm{d}^{2}{\boldsymbol{x}}$ is the spacetime differential volume.

In addition to the time and position operators, we introduce the frequency operator $\widehat{\omega}$, such that $\widehat{\omega}\doteq{\rm i}\epsilon\partial_{t}$ in the coordinate representation. Likewise, the wavevector operator $\widehat{\mathscr{{\boldsymbol{k}}}}$ is defined as $\widehat{\mathscr{{\boldsymbol{k}}}}\doteq-{\rm i}\epsilon{\boldsymbol{\nabla}}$. In particular, these operators are Hermitian, and one has $\braket{\,\mathsf{x}\mid\widehat{\omega}\mid\widetilde{w}\,}={\rm i}\epsilon\partial_{t}\widetilde{w}$ and $\braket{\,\mathsf{x}\mid\widehat{{\boldsymbol{k}}}\mid\widetilde{w}\,}=-{\rm i}\epsilon{\boldsymbol{\nabla}}\widetilde{w}$. Using the relation between the fluctuating generalized vorticity and electric potential, one has $\ket{\widetilde{w}}=-\widehat{k}_{\rm D}^{2}\ket{\widetilde{\psi}}$, where

[TABLE]

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

[TABLE]

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

[TABLE]

where $\smash{\widehat{U}\doteq U(\widehat{t},\widehat{y})}$, and the prime above $U$ henceforth denotes $\partial_{y}$; in particular, $\smash{\widehat{U}^{\prime\prime}\doteq\partial_{y}^{2}\,U(\widehat{t},\widehat{y})}$. Also, $\smash{\ket{\,\widetilde{\xi}\,}}$ is the ket corresponding to the random forcing $\smash{\widetilde{\xi}}$. It is to be noted that the $\smash{\widehat{\mathscr{D}}}$ includes the nonlinear coupling between the ZFs and the DWs.

Additionally in Eq. (12), the kets $\ket{\,f_{\rm nl}[\phi,\psi]\,}$ and $\ket{\,\overline{f_{\rm nl}[\phi,\psi]}\,}$ are given by

[TABLE]

where $\widehat{\mathscr{K}}(\mathsf{x})$ is an operator describing the nonlinear eddy–eddy interactions

[TABLE]

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

[TABLE]

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

[TABLE]

which is one of the nonlinear terms appearing in Eq. (9a). Here we substituted the coordinate representation for the wavevector operator and used $\smash{\ket{\widetilde{w}}=-\widehat{k}_{\rm D}^{2}\ket{\widetilde{\psi}}}$. Note that the factor 2 appears because we have written the ket $\ket{\,f_{\rm nl}[\phi,\psi]\,}$ in a symmetric form so that

[TABLE]

for any two real fields $\phi$ and $\psi$.

Following Ruiz et al. (2016), the fluctuating terms appearing in Eq. (9b) can also be rewritten in the abstract representation. Upon noting that $\widetilde{v}_{x}\widetilde{v}_{y}=(-\epsilon\partial_{y}\widetilde{\psi})(\epsilon\partial_{x}\widetilde{\psi})=-\braket{\,\mathsf{x}\mid\widehat{k}_{x}\mid\widetilde{\psi}\,}\braket{\,\widetilde{\psi}\mid\widehat{k}_{y}\mid\mathsf{x}\,}$, one obtains

[TABLE]

3 Statistical closure

3.1 Statistical-closure problem

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

[TABLE]

where the Hermitian operator $\ket{\,\widetilde{w}\,}\bra{\,\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 $\smash{\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 Appendix A.2. By using the identity (90), one can write Eq. (20) in terms of the correlation operator:

[TABLE]

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

[TABLE]

Subtracting from Eq. (23) its Hermitian conjugate gives

[TABLE]

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

[TABLE]

Both $\widehat{\mathscr{D}}_{\rm H}$ and $\widehat{\mathscr{D}}_{\rm A}$ are Hermitian. Also, $[\cdot,\cdot]_{\mp}$ denote the commutators and anticommutators; i.e., $\smash{[\widehat{\mathscr{A}},\widehat{\mathscr{B}}]_{-}=\widehat{\mathscr{A}}\widehat{\mathscr{B}}-\widehat{\mathscr{B}}\widehat{\mathscr{A}}}$ and $\smash{[\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 Eq. (24), this result is similar to that obtained in Ruiz et al. (2016).

Equations (22) and (24) are not closed. The left-hand side of Eq. (24) is written in terms of $\smash{\widehat{\mathscr{W}}}$, which is bilinear in the fluctuating vorticity. However, the right-hand side of Eq. (24) 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 Eq. (24) in terms of the correlation operator $\smash{\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 Eq. (24). Then, one linearizes Eq. (12) so that

[TABLE]

After formally inverting $\widehat{\mathscr{D}}$, we have $\ket{\,\widetilde{w}\,}\simeq{\rm i}\epsilon\widehat{\mathscr{D}}^{-1}\ket{\,\widetilde{\xi}\,}$. Substituting into Eq. (24) leads to the closed equation

[TABLE]

where $\widehat{\mathscr{S}}\doteq\overline{\ket{\,\widetilde{\xi}\,}\bra{\,\widetilde{\xi}\,}}$ 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 $(-{\rm i}\,\widehat{\mathscr{A}})_{\rm H}=\widehat{\mathscr{A}}_{\rm A}$.

Equations (22) and (27) now form a closed system. In this approximation, the equation for the DW fluctuations was linearized in Eq. (26), but the DW nonlinearity is only kept in the equation for the mean field [Eq. (22)]. This constitutes the QL approximation. If one projects Eq. (27) on the double-physical coordinate space $(\mathsf{x},\mathsf{x^{\prime}})$ using multiplication by $\bra{\,t,{\boldsymbol{x}}\,}$ and $\ket{\,t,{\boldsymbol{x}}^{\prime}\,}$, 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 Eq. (27) on the ray space $(t,{\boldsymbol{x}},\omega,{\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_{\rm nl}$ in Eq. (24). Let us separate $\widetilde{w}$ into three components:

[TABLE]

Here $\ket{\,\widetilde{w}_{0}\,}=O(1)$ is chosen to satisfy the linear part of Eq. (12); i.e., $\widehat{\mathscr{D}}\ket{\,\widetilde{w}_{0}\,}=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{\xi}$. When formally inverting $\widehat{\mathscr{D}}$ in Eq. (12), one finds that, to lowest order in $\epsilon_{\rm nl}$ and $\epsilon$,

[TABLE]

Substituting Eqs. (28) and (29) into Eq. (24) leads to

[TABLE]

Here we neglected $O(\epsilon_{\rm nl}^{3},\epsilon_{\rm nl}^{2}\epsilon,\epsilon_{\rm nl}\epsilon^{2})$ terms, \eg$\epsilon_{\rm nl}^{3}\ket{\,f_{\rm nl}[\widetilde{w}_{0},\widetilde{\phi}]\,}\bra{\,\widetilde{\phi}\,}$ and $\epsilon_{\rm nl}^{2}\epsilon\ket{\,f_{\rm nl}[\widetilde{w}_{0},\widetilde{\phi}]\,}\bra{\,\widetilde{\varphi}\,}$. Note that the zonal averages of quantities involving both $\widetilde{w}_{0}$ and $\widetilde{\xi}$ (\eg$\smash{\epsilon_{\rm nl}\ket{\,\widetilde{\xi}\,}\bra{\,\widetilde{w}_{0}\,}}$ and $\smash{\epsilon_{\rm nl}\epsilon\ket{\,f_{\rm nl}[\widetilde{w}_{0},\widetilde{w}_{0}]\,}\bra{\,\widetilde{\xi}\,}}$) are zero since the statistics of $\widetilde{w}_{0}$ and $\widetilde{\xi}$ are independent. The factor of two in the second term on the right-hand side of Eq. (30) is due to the symmetry property (19).

Now, let us explicitly calculate the statistical average of the nonlinear terms in Eq. (30). To do this, we shall use the quasinormal approximation which expresses higher-order $(n\geq 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

[TABLE]

[TABLE]

where $\widehat{\mathscr{W}}_{0}\doteq\overline{\ket{\,\widetilde{w}_{0}\,}\bra{\,\widetilde{w}_{0}\,}}$ and $\mathsf{x}_{i}=(t_{i},{\boldsymbol{x}}_{i})$ denotes a spacetime coordinate.

Upon substituting Eq. (31), one finds that the first term appearing in the right-hand side of Eq. (30) vanishes; i.e., $\overline{\ket{\,f_{\rm nl}[\widetilde{w}_{0},\widetilde{w}_{0}]\,}\bra{\,\widetilde{w}_{0}\,}}=0$. After substituting Eqs. (14) and (29a), one obtains the following expression for the second term on the right-hand side of Eq. (30):

[TABLE]

where $\ket{\,\mathsf{x}\,}$ and $\ket{\,\mathsf{y}\,}$ are eigenstates of the coordinate operator. Inserting Eq. (32) in the abstract representation gives

[TABLE]

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\mathrm{d}^{3}\mathsf{x}^{\prime}\ket{\,\mathsf{x}^{\prime}\,}\bra{\,\mathsf{x}^{\prime}\,}$ and then using Eq. (32) explicitly. In the third line, we used the symmetry property given in Eq. (19). Substituting this result into Eq. (33) gives

[TABLE]

In a similar manner, substituting Eq. (29a) into the third term in Eq. (30) leads to

[TABLE]

As before, when using the quasinormal approximation, we obtain

[TABLE]

where we introduced the trace operation $\smash{\mathrm{Tr}\,\widehat{\mathscr{A}}\,=\int\mathrm{d}^{3}\mathsf{x}\braket{\,\mathsf{x}\mid\widehat{\mathscr{A}}\mid\mathsf{x}\,}}$. By using the identity operator $\smash{\widehat{1}=\int\mathrm{d}^{3}\mathsf{x}^{\prime}\ket{\,\mathsf{x}^{\prime}\,}\bra{\,\mathsf{x}^{\prime}\,}}$, we are able to write $\braket{\,w_{0}\mid\widehat{\mathscr{A}}\mid w_{0}\,}=\mathrm{Tr}(\,\widehat{\mathscr{A}}\ket{\,w_{0}\,}\bra{\,w_{0}\,}\,)$. Also, note that the trace $\mathrm{Tr}[\,\widehat{\mathscr{K}}(\mathsf{x})\,\widehat{\mathscr{W}}_{0}\,\widehat{\mathscr{K}}^{\dagger}(\mathsf{y})\,\widehat{\mathscr{W}}_{0}\,]$ is not an operator but rather a function of the spacetime coordinates $\mathsf{x}$ and $\mathsf{y}$. Upon inserting Eq. (39) into Eq. (38), we obtain

[TABLE]

We then substitute Eqs. (35) and (40) into Eq. (30). 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:

[TABLE]

where

[TABLE]

In summary, Eqs. (22), (41), and (42) form a complete set of equations for the zonal flow $U(t,y)$ and for the correlation operator $\smash{\widehat{\mathscr{W}}}$. Multiplying Eq. (41) by $\bra{\,t,{\boldsymbol{x}}\,}$ and $\ket{\,t^{\prime},{\boldsymbol{x}}^{\prime}\,}$ leads to the “quasinormal” equation for the two-point correlation function written in the double-physical coordinate space $(\mathsf{x},\mathsf{x^{\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 Appendix 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

[TABLE]

where $\mathsf{k}\doteq(\omega,{\boldsymbol{k}})$, $\mathsf{s}\doteq(\tau,{\boldsymbol{s}})$, $\mathsf{k}\cdot\mathsf{s}=\omega\tau-{\boldsymbol{k}}\cdot{\boldsymbol{s}}$, and $\mathrm{d}^{3}\mathsf{s}\doteq\mathrm{d}\tau\,\mathrm{d}^{2}{\boldsymbol{s}}$. Here the integrals span over $\mathbb{R}^{3}$. The Weyl symbol $A(\mathsf{x},\mathsf{k})$ of an operator $\smash{\widehat{\mathscr{A}}}$ is a function on the extended six-dimensional phase space $(\mathsf{x},\mathsf{k})=(t,{\boldsymbol{x}},\omega,{\boldsymbol{k}})$; i.e., $A=A(t,{\boldsymbol{x}},\omega,{\boldsymbol{k}})$. Physically, $A(\mathsf{x},\mathsf{k})$ can be interpreted as a local Fourier transform of the spacetime representation of the operator $\smash{\widehat{\mathscr{A}}}$ [see, for instance, Eq. (73)].

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

[TABLE]

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

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

[TABLE]

Here $D(t,y,\omega,{\boldsymbol{k}})$ is the Weyl symbol corresponding to $\smash{\widehat{\mathscr{D}}}$. From the properties of the Weyl transform, we observe that $D_{\rm H}=\mathrm{Re}\,D$ and $D_{\rm A}=\mathrm{Im}\,D$ are the (real) Weyl symbols corresponding to the operators $\smash{\widehat{\mathscr{D}}_{\rm H}}$ and $\smash{\widehat{\mathscr{D}}_{\rm A}}$, respectively. (“$\mathrm{Re}$” and “$\mathrm{Im}$” denote the real and imaginary parts, respectively.) These are given by

[TABLE]

In addition, $F(t,y,\omega,{\boldsymbol{k}})$, $\eta(t,y,\omega,{\boldsymbol{k}})$, and $[D^{-1}]^{*}(t,y,\omega,{\boldsymbol{k}})$ are the Weyl symbols corresponding to $\smash{\widehat{\mathscr{F}}}$, $\smash{\widehat{\mathscr{\eta}}}$, and $\smash{(\widehat{\mathscr{D}}^{-1})^{\dagger}}$, respectively. (These will be calculated explicitly later.) The Moyal product (75) is denoted by “$\star$”, and the brackets $\boldsymbol{\{}\!\{\cdot,\cdot\}\!\boldsymbol{\}}$ and $\boldsymbol{[}[\cdot,\cdot]\boldsymbol{]}$ are the Moyal brackets [Eqs. (79) and (80)]. Basic properties of the Moyal product and the Moyal brackets are given in Appendix A.1.

Modulo the statistical closure introduced in §3, Eq. (45) is an exact equation for the dynamics of the Wigner function. However, Eq. (45) is difficult to solve as is, both analytically and numerically. Thus, we shall reduce Eq. (45) 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 Eq. (45), we shall expand the Moyal products and brackets in terms of the ordering parameter $\epsilon$. As a reminder, the parameter $\epsilon$ was introduced in Eq. (8) in order to denote the spatio-temporal scale separation between the DW and ZF dynamics. From the definition of the Moyal product [Eq. (76)], one has

[TABLE]

where $\overleftrightarrow{\mathcal{L}}$ is the Janus operator [Eq. (77)], which basically serves as the canonical Poisson bracket in the extended six-dimensional phase space $(t,{\boldsymbol{x}},\omega,{\boldsymbol{k}})$. Likewise, the Moyal brackets in Eqs. (79) and (80) are approximated by

[TABLE]

For the purposes of this paper, higher-order corrections to Eqs. (47) and (48) 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 Eq. (23) can be written as $\smash{\widehat{\mathscr{D}}_{\rm H}\widehat{\mathscr{W}}=O(\epsilon,\epsilon_{\rm nl})}$, where we assumed small dissipation $[D_{\rm A}\sim O(\epsilon)]$. In the Weyl representation, this equation becomes $D_{\rm H}\star W=O(\epsilon,\epsilon_{\rm nl})$. Using Eq. (47) gives

[TABLE]

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

[TABLE]

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

To remain consistent with the asymptotic expansion in Eqs. (47) and (48), the symbols $D_{\rm H}$ and $D_{\rm A}$ in Eqs. (46) are approximated as well to lowest order in $\epsilon$:

[TABLE]

where we substituted Eqs. (48).

Let us now obtain an approximate expression for $[D^{-1}]^{*}$ in Eq. (45). As a reminder, $[D^{-1}]$ and $D^{-1}$ are not the equivalent; $[D^{-1}]$ is the Weyl symbol of $\smash{\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 $\smash{\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}]^{*}$ with its limiting form as $\epsilon$ tends to zero:

[TABLE]

where “$\mathcal{P}$” denotes the Cauchy principal value.

We then insert Eqs. (48)–(52) into Eq. (45) and integrate over the frequency variable $\omega$. Afterwards, we expand the Moyal products and brackets using Eqs. (47) and (48). Although the GO ansatz (50) and the expression for $[D^{-1}]$ in Eq. (52) involve Dirac delta functions whose derivatives in phase space are not smooth, in Appendix 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 $\epsilon$:

[TABLE]

Note that in Eq. (LABEL:eq:wke_aux) the Weyl symbols $W$, $F$, and $S$ are real because their associated operators are Hermitian.

For simplicity, we consider that the nonlinear scaling parameter $\epsilon_{\rm nl}$ and the GO parameter $\epsilon$ scale as $\epsilon\sim\epsilon_{\rm nl}$. Upon integrating Eq. (LABEL:eq:wke_aux) over the frequency and neglecting higher-order terms in $\epsilon$, we obtain the collisional wave kinetic equation (cWKE)

[TABLE]

where $\{\cdot,\cdot\}=\overleftarrow{\partial_{\boldsymbol{x}}}\cdot\overrightarrow{\partial_{\boldsymbol{k}}}-\overleftarrow{\partial_{\boldsymbol{k}}}\cdot\overrightarrow{\partial_{\boldsymbol{x}}}$. The wave frequency $\Omega$ and dissipation term $\Gamma$ are respectively obtained from the lowest-order expansion in $\epsilon$ of $D_{\rm H}$ and $D_{\rm A}$ in Eqs. (51):

[TABLE]

The external source term $S_{\rm ext}(t,y,{\boldsymbol{k}})$ appearing in Eq. (57) is given by

[TABLE]

Assuming that $\widetilde{\xi}(t,{\boldsymbol{x}})$ is white noise leads to

[TABLE]

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

[TABLE]

Here we used $\Xi(y,{\boldsymbol{s}})=\Xi(y,-{\boldsymbol{s}})$, which is due to the fact that $\widetilde{\xi}$ is real by definition.

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

[TABLE]

As shown in Appendix B.2, to the leading order in $\epsilon$, $\gamma_{\rm nl}[J]$ and $S_{\rm nl}[J,J]$ are given by

[TABLE]

Here $\Theta(t,y,{\boldsymbol{k}},{\boldsymbol{p}},{\boldsymbol{q}})\doteq\upi\delta(\Delta\Omega)$, and

[TABLE]

One can identify $\Delta\Omega=0$ as the frequency-resonance condition. Finally, the kernel $M({\boldsymbol{p}},{\boldsymbol{q}})$ in Eqs. (63) is

[TABLE]

4.3 Dynamics of the ZF velocity $U$

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

[TABLE]

where we used the Moyal product and substituted Eq. (50). Upon following Ruiz et al. (2016), we substitute $\epsilon\sim\epsilon_{\rm nl}$ and obtain to lowest order in $\epsilon$:

[TABLE]

5 Discussion

5.1 Main equations

Equations (57) and (67), as well as Eqs. (58), (62), and (63), 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 (57) governs the dynamics of the wave-action density $J$ for DWs. The left-hand side of Eq. (57) describes the wave refraction governed by the wave frequency $\Omega$ [Eq. (58a)], which serves as a Hamiltonian for the system. On the right-hand side, $\mu_{\rm dw}$ represents weak dissipation due to the external environment, and $\Gamma$ denotes linear dissipation caused by the ZFs (Parker, 2018; Ruiz et al., 2016). The term $S_{\rm ext}$ represents an external source term for the DW fluctuations.

The nonlinear term $C[J,J]$ in Eq. (57) plays the role of a wave scattering operator. It is composed of two terms, $\gamma_{\rm nl}$ and $S_{\rm nl}$, which arise from nonlinear wave–wave interactions. The nonlinear source term $S_{\rm nl}$ in (63b) 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 $\gamma_{\rm nl}$ in Eq. (63a) 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 $\gamma_{\rm nl}$ are called the coherent response (Krommes, 2002). The other terms in Eqs. (57) and (67) are the same as in the QL theory, and their physical meaning was already discussed elsewhere (Parker, 2016; Ruiz et al., 2016; Zhu et al., 2018c; Parker, 2018; Zhu et al., 2018a, b).

5.2 Conservation properties

Equations (57) and (67) inherit the same conservation laws of the original gHME (2). In other words, for isolated systems ($S_{\rm ext}=0$ and $\mu_{\rm zf,dw}=0$), Eqs. (57) and (67) conserve the total enstrophy and total energy

[TABLE]

where we used $\epsilon\sim\epsilon_{\rm nl}$. Also, the expressions for the DW and ZF components of the enstrophy and energy are

[TABLE]

These conservation laws are proven in Appendix 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., 2005; Smolyakov & Diamond, 1999; Smolyakov et al., 2000; Malkov & Diamond, 2001; Malkov et al., 2001; Diamond et al., 1994; Kim & Diamond, 2003; Kaw et al., 2002; Trines et al., 2005; Singh et al., 2014) with $\Omega=k_{x}U-\beta k_{y}/k_{\rm D}^{2}$ and $\Gamma=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 in the structure of the DW phase space and its relation to the Rayleigh–Kuo parameter $U^{\prime\prime}/\beta$ (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, \egKrommes (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, Eqs. (57) and (67) make a distinction between the statistics of the DWs and the ZF velocity field: the DW component is modeled 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 Eq. (62) coincides with those derived previously within homogeneous WTT [see, \egConnaughton et al. (2015)]. Hence, the cWKE (57) 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, \egSec. 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 nonzero $U$. In particular, the frequency-resonance condition in WTT depends on a simplified expression for the DW frequency: $\Omega(t,y,{\boldsymbol{k}})\simeq-\beta k_{x}/k_{\rm D}^{2}$, which neglects the ZF velocity. Hence, one has

[TABLE]

In contrast, in the cWKE, $\Omega$ is given by Eq. (58a), so

[TABLE]

where the contribution of the Doppler shifts canceled out due to the spatial-resonance condition $({\boldsymbol{k}}={\boldsymbol{p}}+{\boldsymbol{q}})$. The factor $|\beta-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 (Zhu et al., 2018c; Kuo, 1949; Numata et al., 2007). From this result, it seems that the present theory breaks down in regions where $\beta-U^{\prime\prime}\simeq\epsilon$. 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 et al. (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 (Parker, 2018; Ruiz et al., 2016; Zhu et al., 2018b; Parker, 2016; Zhu et al., 2018a) 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.

The authors thank H. Zhu, Y. Zhou, and J. B. Parker for helpful discussions. This work was supported by the U.S. DOE through Contract DE-AC02-09CH11466 and by Sandia National Laboratories. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. DOE National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. DOE or the U.S. Government.

Appendix A Conventions and definitions

A.1 Weyl calculus

This appendix summarizes our conventions used for the Weyl calculus. For more information, see the reviews by Ruiz (2017), Tracy et al. (2014), Imre et al. (1967), Baker Jr. (1958), and McDonald (1988).

Let $\widehat{\mathscr{A}}$ be an operator defined on the Hilbert space $L^{2}(\mathbb{R}^{3})$ with the inner product (10). The Weyl symbol $A(\mathsf{x},\mathsf{k})$ is defined as the Weyl transform $\mathsf{W}[\widehat{\mathscr{A}}]$ of $\widehat{\mathscr{A}}$; namely,

[TABLE]

where $\ket{\,\mathsf{x}\,}$ are the eigenstates of the position operator. Also, $\mathsf{k}\doteq(\omega,{\boldsymbol{k}})$, $\mathsf{s}\doteq(\tau,{\boldsymbol{s}})$, $\mathsf{k}\cdot\mathsf{s}=\omega\tau-{\boldsymbol{k}}\cdot{\boldsymbol{s}}$, and $\mathrm{d}^{3}\mathsf{s}\doteq\mathrm{d}\tau\,\mathrm{d}^{2}{\boldsymbol{s}}$. The integrals span $\mathbb{R}^{3}$. This description of the operators is known as the phase-space representation since the Weyl symbols are functions of the six-dimensional phase space. Conversely, the inverse Weyl transform is given by

[TABLE]

where $\mathrm{d}^{3}\mathsf{k}\doteq\mathrm{d}\omega\,\mathrm{d}^{2}{\boldsymbol{k}}$. The coordinate representation $\mathscr{A}(\mathsf{x},\mathsf{x}^{\prime})=\braket{\,\mathsf{x}\mid\widehat{\mathscr{A}}\mid\mathsf{x}^{\prime}\,}$ of $\widehat{\mathscr{A}}$ is

[TABLE]

In the following, we shall outline a number of useful properties of the Weyl transform.

(i) The trace $\mathrm{Tr}\,\widehat{\mathscr{A}}\,\doteq\int\mathrm{d}^{3}\mathsf{x}\,\braket{\,\mathsf{x}\mid\widehat{\mathscr{A}}\mid\mathsf{x}\,}$ of $\widehat{\mathscr{A}}$ is

[TABLE]

(ii) If $A(\mathsf{x},\mathsf{k})$ is the Weyl symbol of $\widehat{\mathscr{A}}$, then $A^{*}(\mathsf{x},\mathsf{k})$ is the Weyl symbol of $\widehat{\mathscr{A}}^{\dagger}$. As a corollary, if the operator $\widehat{\mathscr{A}}$ is Hermitian $(\widehat{\mathscr{A}}=\widehat{\mathscr{A}}^{\dagger})$, then $A(\mathsf{x},\mathsf{k})$ is real.

(iii) For linear operators $\widehat{\mathscr{A}}$, $\widehat{\mathscr{B}}$ and $\widehat{\mathscr{C}}$ where $\widehat{\mathscr{C}}=\widehat{\mathscr{A}}\widehat{\mathscr{B}}$, the corresponding Weyl symbols satisfy

[TABLE]

Here “$\star$” refers to the Moyal product (Moyal, 1949), which is

[TABLE]

Also, $\overleftrightarrow{\mathcal{L}}$ is the Janus operator

[TABLE]

The arrows indicate the direction in which the derivatives act. Note that $A\overleftrightarrow{\mathcal{L}}B$ serves as the canonical Poisson bracket in the extended six-dimensional phase space $(t,{\boldsymbol{x}},\omega,{\boldsymbol{k}})$.

(iv) The Moyal product is associative; i.e., for arbitrary symbols $A$, $B$, and $C$, one has

[TABLE]

(v) The anti-symmetrized Moyal product defines the so-called Moyal bracket, namely,

[TABLE]

Likewise, the symmetrized Moyal product is defined as

[TABLE]

(vi) For fields that vanish rapidly enough at infinity, when integrated over all phase space, the Moyal product of two symbols equals the regular product; i.e.,

[TABLE]

(vii) Now we tabulate some Weyl transforms of various operators. First, the Weyl transform of the identity is

[TABLE]

The Weyl transforms of the time and position operators are given by

[TABLE]

where $\widehat{t}=t$ and $\widehat{{\boldsymbol{x}}}={\boldsymbol{x}}$ in the coordinate representation. Likewise,

[TABLE]

where $\widehat{\omega}={\rm i}\epsilon\partial_{t}$ and $\widehat{{\boldsymbol{k}}}=-{\rm i}\epsilon{\boldsymbol{\nabla}}$ in the coordinate representation. For any two operators $f(\widehat{t},\widehat{{\boldsymbol{x}}})$ and $g(\widehat{\omega},\widehat{{\boldsymbol{k}}})$,

[TABLE]

Upon using the Moyal product (76), one has

[TABLE]

A.2 Zonal average of operators

The zonal average $\overline{\widehat{\mathscr{A}}}$ of any given operator $\widehat{\mathscr{A}}$ is defined through the Weyl calculus (Appendix A.1) and is schematically shown in Fig. 1. First, by using the Weyl transform (71), one calculates the Weyl symbol $A(t,{\boldsymbol{x}},\omega,{\boldsymbol{k}})$ [Eq. (71)] corresponding to the operator $\smash{\widehat{\mathscr{A}}}$. Then, one calculates the zonal average defined in §2.2 on the Weyl symbol $A(t,{\boldsymbol{x}},\omega,{\boldsymbol{k}})$. This leads to

[TABLE]

The zonal-averaged operator $\smash{\overline{\widehat{\mathscr{A}}}}$ is obtained by applying the inverse Weyl transform (72) on $\overline{A}(t,y,\omega,{\boldsymbol{k}})$.

From Eq. (73), the coordinate representation of $\overline{\widehat{\mathscr{A}}}$ is

[TABLE]

where $\overline{A}_{\rm F}$ is the inverse Fourier transform on the frequency and wavevector variables of the Weyl symbol $\overline{A}(t,y,\omega,{\boldsymbol{k}})$. Also, $\overline{t}\doteq(t+t^{\prime})/2$, $\overline{y}\doteq(y+y^{\prime})/2$, $\Delta t\doteq t-t^{\prime}$, and $\Delta{\boldsymbol{x}}\doteq{\boldsymbol{x}}-{\boldsymbol{x}}^{\prime}$.

Another useful identity that we shall use is

[TABLE]

where we substituted Eq. (73) in the first and last lines.

Appendix B Auxiliary calculations

B.1 Simplifying the Moyal products

To derive the WKE (57), we need to approximate the Moyal products appearing in Eq. (45). The most difficult terms to approximate are those involving derivatives of Dirac delta functions. As an example, in this appendix we calculate the integral

[TABLE]

where $A(z)$ is an arbitrary function, $W(z)=2\upi\epsilon\delta\boldsymbol{(}D_{\rm H}(z)\boldsymbol{)}J(t,{\boldsymbol{x}},{\boldsymbol{k}})$ is the GO ansatz (50), and $z\doteq(t,{\boldsymbol{x}},\omega,{\boldsymbol{k}})$. From Eqs. (51a) and (58a), the GO dispersion relation is $D_{\rm H}(z)\simeq\omega-\Omega(t,{\boldsymbol{x}},{\boldsymbol{k}})=0$. Substituting Eq. (50) into Eq. (91) leads to

[TABLE]

where $\overleftrightarrow{\mathcal{L}}$ is the Janus operator (77) and

[TABLE]

Now, let us calculate each of the terms appearing in Eq. (92). The $n=0$ term is simply given by

[TABLE]

For the $n=1$ term in Eq. (92), we write the Janus operator as $\overleftrightarrow{\mathcal{L}}\doteq\overleftarrow{\partial_{\mu}}\mathcal{J}^{\mu\nu}\overrightarrow{\partial_{\nu}}$, where $\mathcal{J}^{\mu\nu}$ is the canonical Poisson tensor in $z$ space. We then obtain

[TABLE]

Here we integrated by parts in $\omega$ and used the fact that $J$ and $\Omega$ are independent of $\omega$. Explicitly writing the derivatives in terms of phase-space coordinates leads to

[TABLE]

where $\{\cdot,\cdot\}\doteq\overleftarrow{\partial_{\boldsymbol{x}}}\cdot\overrightarrow{\partial_{\boldsymbol{k}}}-\overleftarrow{\partial_{\boldsymbol{k}}}\cdot\overrightarrow{\partial_{\boldsymbol{x}}}$ is the canonical Poisson bracket.

After integrating by parts, we have been able to express $\mathcal{I}_{1}$ in Eq. (98) as a sum of derivatives acting on smooth functions. Since no singularities are present, then $\mathcal{I}_{1}=O(\epsilon^{0})$. By reiterating the same calculation in Eq. (95), one can show that the higher $\mathcal{I}_{n}$ terms in Eq. (92) are also smooth. Thus, one can approximate integrals as in Eq. (91) by their leading-order expansion in $\epsilon$. With these results, one can simplify the Moyal products and brackets in Eq. (45) in order to obtain the WKE (57).

In particular, let us calculate the special case where $A=D_{H}=\omega-\Omega$. Upon using Eq. (98), we find

[TABLE]

which is the advection term appearing in the WKE (57).

B.2 Calculation of the nonlinear dissipation rate and the nonlinear source term

We begin by calculating the Weyl symbol $F(t,{\boldsymbol{x}},\omega,{\boldsymbol{k}})$ corresponding to $\widehat{\mathscr{F}}$ in Eq. (42b). Upon substituting Eq. (16), we first note that the trace operator appearing in Eq. (42b) can be written as

[TABLE]

where the operators $\widehat{\mathscr{L}}_{j}$ and $\widehat{\mathscr{R}}_{j}$ are defined in Eq. (17).

After substituting Eq. (100) into Eq. (42b) and applying the Weyl transform (71), we obtain integrals of the form

[TABLE]

where $\widehat{\mathscr{A}}$ and $\widehat{\mathscr{B}}$ represent the terms appearing in Eq. (100). To evaluate these integrals, we shall use the identity

[TABLE]

(This property is analogous to the convolution theorem frequently used in the Fourier transform.) Hence, we have

[TABLE]

where we used the Moyal product (75). Also, $L_{j}({\boldsymbol{k}})=-({\boldsymbol{e}}_{z}\times{\boldsymbol{k}})_{j}k_{D}^{-2}$ and $R_{j}({\boldsymbol{k}})={\boldsymbol{k}}_{j}$ are respectively the Weyl symbols corresponding to the operators $\smash{\widehat{\mathscr{L}}_{j}}$ and $\smash{\widehat{\mathscr{R}}_{j}}$ in Eqs. (17). The wave–wave nonlinearities are considered to be weak; hence, the Moyal products in Eq. (105) can be replaced by ordinary products. Hence,

[TABLE]

where we used the reality property of $\widetilde{\psi}$ so $W(\mathsf{x},\mathsf{q})=W(\mathsf{x},-\mathsf{q})$ and introduced $M({\boldsymbol{p}},{\boldsymbol{q}})\doteq L_{j}({\boldsymbol{p}})R_{j}({\boldsymbol{q}})+L_{k}({\boldsymbol{q}})R_{k}({\boldsymbol{p}})={\boldsymbol{e}}_{z}\cdot({\boldsymbol{p}}\times{\boldsymbol{q}})\left(q_{\rm D}^{-2}-p_{\rm D}^{-2}\right)$. [The Wigner function $W$ is a function of $(t,y,\omega,{\boldsymbol{k}})$ only. However, to simplify our notation, we used the arguments $(\mathsf{x},\mathsf{k})$ to denote the dependence on the phase-space variables.] Substituting the GO ansatz in Eq. (50) into Eq. (106) and integrating in the frequency variables leads to

[TABLE]

Substituting Eq. (109) into Eq. (LABEL:eq:wke_aux) and integrating in $\omega$ leads to the term $S_{\rm nl}(t,y,{\boldsymbol{k}})$ reported in Eq. (63b).

Now, let us compute the Weyl symbol $\eta(\mathsf{x},\mathsf{k})$ of $\widehat{\mathscr{\eta}}$ in Eq. (42a). Substituting Eq. (16) into Eq. (42a) leads to

[TABLE]

To calculate the Weyl transform of the above, we shall use the following result:

[TABLE]

where in the third line we substituted Eq. (101). We then calculate $\eta(\mathsf{x},\mathsf{k})$ in Eq. (110) by using Eq. (111). Similarly as in Eq. (109), we later approximate the Moyal products by ordinary products. To leading order, we obtain

[TABLE]

From Eq. (52), we approximate $\mathrm{Im}\{[D^{-1}]^{*}(\mathsf{x},-\mathsf{q})\}\simeq\upi\delta\boldsymbol{(}q_{0}-\Omega(t,y,{\boldsymbol{q}})\boldsymbol{)}$. When substituting $W(t,y,\omega,{\boldsymbol{k}})=2\upi\epsilon\delta\boldsymbol{(}\omega-\Omega(t,y,{\boldsymbol{k}})\boldsymbol{)}J(t,y,{\boldsymbol{k}})$ and $M({\boldsymbol{p}},-{\boldsymbol{q}})=-M({\boldsymbol{p}},{\boldsymbol{q}})$, one obtains the following:

[TABLE]

Finally, substituting Eq. (117) into Eq. (LABEL:eq:wke_aux) and integrating in $\omega$ leads to the term $\gamma_{\rm nl}(t,y,{\boldsymbol{k}})$ reported in Eq. (63a).

B.3 Conservation of the total enstrophy $\mathcal{Z}$ and the total energy $\mathcal{E}$

The time derivatives of the total enstrophy $\mathcal{Z}$ and total energy $\mathcal{E}$ are

[TABLE]

where we substituted Eqs. (57), (67), and (69). As shown by Ruiz et al. (2016), the first integrals on the right-hand side of the equations above are zero. Hence, in order to show conservation of total enstrophy and energy, one needs to show that

[TABLE]

where $\sigma({\boldsymbol{k}})=1$ for the case of enstrophy and $\sigma({\boldsymbol{k}})=k_{\rm D}^{-2}$ for the case of energy. Since $C[J,J]$ is of the canonical form of the scattering operators found in homogeneous-turbulence theories, the proof of its conservation properties closely follows that sketched in §4.2.4 of Krommes (2002). Substituting Eq. (62) leads to

[TABLE]

One can then exchange the momentum variables since they are simply integration variables. After using the symmetry property of $M({\boldsymbol{p}},{\boldsymbol{q}})$ and the identity $\Theta(t,y,{\boldsymbol{q}},-{\boldsymbol{p}},{\boldsymbol{k}})=\Theta(t,y,{\boldsymbol{k}},{\boldsymbol{p}},{\boldsymbol{q}})$, one can rewrite $\mathcal{G}$ as follows:

[TABLE]

where

[TABLE]

When substituting the spatial-resonance condition ${\boldsymbol{k}}={\boldsymbol{p}}+{\boldsymbol{q}}$, one shows that $\chi({\boldsymbol{k}},{\boldsymbol{p}},{\boldsymbol{q}})=0$ for both $\sigma({\boldsymbol{k}})=1$ and $\sigma({\boldsymbol{k}})=k_{\rm D}^{-2}$. Hence, $\mathcal{G}=0$, which proves that Eqs. (57) and (67) conserve total enstrophy $\mathcal{Z}$ and total energy $\mathcal{E}$.