2.2.1. Generalized (broadened) kinetic equation

In the limit of small nonlinearity, one develops a perturbation expansion in the nonlinearity strength, which leads, under certain assumptions, to the wave turbulence kinetic equation. The derivation of the resonant kinetic equation is well understood and studied, see Zakharov et al. (Reference Zakharov, Lvov and Falkovich1992) and Nazarenko (Reference Nazarenko2011). Taking non-resonant interactions leads to a different version of the kinetic equation with the frequency delta functions being replaced by a Lorentzian, see Lvov et al. (Reference Lvov, Lvov, Newell and Zakharov1997) and Lvov, Polzin & Yokoyama (Reference Lvov, Polzin and Yokoyama2012). This derivation also hinges on the assessment that fourth-order cumulants are a subleading term compared with the product of two double correlators (Deng & Hani Reference Deng and Hani2023). For the three-wave Hamiltonian (2.14), the kinetic equation is (2.16), describing general internal waves interacting in both rotating and non-rotating environments

(2.16) \begin{align} \frac{\partial}{\partial t} n_{{\boldsymbol{p}}} &= \int \int {\rm d}{\boldsymbol{p}}_1 \,{\rm d}{\boldsymbol{p}}_2 \big( |V_{{\boldsymbol{p}}_1,{\boldsymbol{p}}_2}^{{\boldsymbol{p}}}|^2\delta(\kern1pt {\boldsymbol{p}}-{\boldsymbol{p}}_1-{\boldsymbol{p}}_2){\mathcal{L}} (\Delta \sigma_{p12}, \varGamma_{p12}) [n_{{\boldsymbol{p}}_1}n_{{\boldsymbol{p}}_2}-n_{{\boldsymbol{p}}}[n_{{\boldsymbol{p}}_1}+n_{{\boldsymbol{p}}_2}]]\nonumber\\ &\quad- |V_{{\boldsymbol{p}}_2,{\boldsymbol{p}}}^{{\boldsymbol{p}}_1} |^2 \delta(\kern1pt {\boldsymbol{p}}-{\boldsymbol{p}}_1+{\boldsymbol{p}}_2) {\mathcal{L}}(\Delta \sigma_{12p}, \varGamma_{p12}) [n_{{\boldsymbol{p}}_2}n_{{\boldsymbol{p}}}-n_{{\boldsymbol{p}}_1}[n_{{\boldsymbol{p}}_2}+n_{{\boldsymbol{p}}}]]\nonumber\\ &\quad -\,|V_{{\boldsymbol{p}},{\boldsymbol{p}}_1}^{{\boldsymbol{p}}_2} |^2 \delta(\kern1pt {\boldsymbol{p}}+{\boldsymbol{p}}_1-{\boldsymbol{p}}_2){\mathcal{L}}(\Delta \sigma_{2p1}, \varGamma_{p12})[n_{{\boldsymbol{p}}}n_{{\boldsymbol{p}}_1}-n_{{\boldsymbol{p}}_2}[n_{{\boldsymbol{p}}}+n_{{\boldsymbol{p}}_1}]] \big) , \end{align}

where Lorentzian $\mathcal {L}$ is given by ${\mathcal {L}} = {\varGamma _{p12}}/({(\Delta \sigma )^2 + \varGamma _{p12}^2})$ and $\Delta \sigma _{p12}=\sigma _p-\sigma _{p_1}-\sigma _{p_2}$ represents the distance from the resonant surface. The resonant manifold is defined by

(2.17) \begin{equation} \left.\begin{gathered} \sigma = \sigma_1 + \sigma_2 ; \quad \boldsymbol{p} = \boldsymbol{p}_1 + \boldsymbol{p}_2,\\ \sigma = \sigma_1 - \sigma_2 ; \quad \boldsymbol{p} = \boldsymbol{p}_1- \boldsymbol{p}_2,\\ \sigma = \sigma_2 - \sigma_1 ; \quad \boldsymbol{p} = \boldsymbol{p}_2- \boldsymbol{p}_1, \end{gathered}\right\} \end{equation}

and appears in figure 1. In wave turbulence theory of internal waves the importance of special extreme scale-separated triads was recognized in McComas & Bretherton (Reference McComas and Bretherton1977). These extreme scale-separated limits are called induced diffusion(ID), elastic scattering (ES) and parametric subharmonic instability (PSI). For an explanation of these triads see also McComas & Müller (Reference McComas and Müller1981a).

Figure 1. The resonant manifold (2.17) in the situation where the three horizontal wavevectors are either parallel or anti-parallel, plotted in a vertical wavenumber–frequency space, for a wave at the centre of the green circle. With rotation, extreme scale separations in horizontal wavenumber lead to the extreme scale-separated triads mentioned in the introduction. These triads are Bragg scattering (also called elastic scattering, or ES, and a phase velocity equals ground velocity resonance condition, called induced diffusion, or ID, being located at the Coriolis frequency $f$. This study focuses on the latter class, with scale separation in both horizontal and vertical wavenumbers. Near-resonant ID conditions are depicted in green, bandwidth-limited non-resonant ID forcing in cyan. The third type of extreme scale separated triads, called parametric subharmonic instability, or PSI, does not play a role in this manuscript.

The total resonance width associated with a specific triad is given by $\varGamma _{p12} = \gamma _p + \gamma _1 +\gamma _2$ and the equation for the individual resonance widths is given by

(2.18) \begin{align} \gamma_p &= \displaystyle \int \int {\rm d}{\boldsymbol{p}}_1\,{\rm d}{\boldsymbol{p}}_2\big( |V_{{\boldsymbol{p}}_1,{\boldsymbol{p}}_2}^{{\boldsymbol{p}}} |^2 \delta(\kern1pt {\boldsymbol{p}}-{\boldsymbol{p}}_{1}-{\boldsymbol{p}}_2) {\mathcal{L}} (\sigma-\sigma_1-\sigma_2) [n_{{\boldsymbol{p}}_1}+n_{{\boldsymbol{p}}_2}] \nonumber\\ &\quad + |V_{{\boldsymbol{p}}_2,{\boldsymbol{p}}}^{{\boldsymbol{p}}_1} |^2 \delta(\kern1pt {\boldsymbol{p}}-{\boldsymbol{p}}_1+{\boldsymbol{p}}_2){\mathcal{L}}(\sigma-\sigma_1+\sigma_2) [n_{{\boldsymbol{p}}_2}-n_{{\boldsymbol{p}}_1}] \nonumber\\ &\quad +\, |V_{{\boldsymbol{p}},{\boldsymbol{p}}_1}^{{\boldsymbol{p}}_2} |^2 \delta(\kern1pt {\boldsymbol{p}}+{\boldsymbol{p}}_1-{\boldsymbol{p}}_2){\mathcal{L}}(\sigma+\sigma_1-\sigma_2) [n_{{\boldsymbol{p}}_1}-n_{{\boldsymbol{p}}_2}]\big). \end{align}

Physically, this $\gamma _p$ represents the fast time scale of decay of a narrow perturbation to the otherwise stationary spectrum (Lvov et al. Reference Lvov, Lvov, Newell and Zakharov1997; Polzin & Lvov Reference Polzin and Lvov2017). It coincides with Langevin rates estimated by Pomphrey, Meiss & Watson (Reference Pomphrey, Meiss and Watson1980) and the decay rate of the McComas (Reference McComas1977) spike experiments. The replacement of the frequency-conserving delta function by the Lorentzian takes into account not only resonant but also near-resonant and non-resonant interactions. Non-resonant interactions appear as a result of the Lorentzian decaying slowly. The role of the non-resonant interactions have to be investigated separately for each particular problem.

2.2.2. Bandwidth estimates

The broadened kinetic equation provides access to a time scale as the inverse of the bandwidth $\gamma _{p}$ (2.18).

Considered in isolation from other interactions and in the resonant limit, the extreme scale-separated ID mechanism gives rise to (Polzin & Lvov Reference Polzin and Lvov2017)

(2.19)\begin{equation} \gamma_p^{ID} = \frac{2}{\rm \pi} \frac{k e_o m_{{\ast}}}{N} \frac{\sigma^2}{f^2} , \end{equation}

where $k$ is horizontal wavenumber magnitude, $e_o$ is the total energy density in the GM model and $m_{\ast }$ is the vertical wavenumber bandwidth parameter; see Appendix A. In the context of a broadened kinetic equation, $\gamma _p^{ID}$ represents the decay of a spike introduced into an otherwise smooth spectrum (Lvov et al. Reference Lvov, Lvov, Newell and Zakharov1997) and is identified as such by McComas & Müller (Reference McComas and Müller1981a) in the numerical experiments of McComas (Reference McComas1977). The associated time scale is far smaller than a wave period and at loggerheads with the intent of having (2.16) describe the slow time evolution of the spectral density.

At finite, but yet small-amplitude, analysis of (2.18) reveals that the bandwidth is composed of both resonant $\gamma _{res}$ and non-resonant terms

(2.20)\begin{equation} \left.\begin{gathered} \text{resonant:} \ \gamma_{res} = \frac{2}{\rm \pi} \frac{k e_0 m_{{\ast}}}{N} \big(\frac{\sigma}{f}\big)^2 \\ \text{non-resonant:} \ \frac{2}{\rm \pi} \frac{k^2 e_0}{f^2} \,\gamma_{res} \end{gathered}\right\}. \end{equation}

The non-resonant contributions come from vertical scales that are larger than those of the high-frequency wave and come from the tails of the Lorentzian admitting non-resonant interactions with the highly energetic near-inertial wave field.

The non-resonant contributions increase at a greater rate with respect to increasing amplitude than resonant contributions. At realistic oceanic amplitudes, numerical evaluation (Polzin & Lvov Reference Polzin and Lvov2017) demonstrates that $\gamma _p^{ID}$ is proportional to the root-mean-square (r.m.s.) Doppler shift. The difficulty of placing a physical interpretation on such a short time decay time scale (2.19) and suspicion that the broadened kinetic equation would suffer from the Galilean invariance that plagued Kraichnan's 1959 DIA model for 3-D turbulence led to scepticism and extensive commentary in the literature (Holloway Reference Holloway1980, Reference Holloway1982). Addressing these issues required development of a broadened kinetic equation, which in turn required waiting for the derivation in canonical coordinates presented at the beginning of this section.

2.2.3. Fokker–Planck diffusion limit

Following McComas & Bretherton (Reference McComas and Bretherton1977) for the resonant kinetic equation, we start from (2.16) and pick off the interactions having ${\boldsymbol {p}}$ nearly parallel to ${\boldsymbol {p}}_1$ with ${\boldsymbol {p}}_2$ small, or nearly parallel to ${\boldsymbol {p}}_2$ with ${\boldsymbol {p}}_1$ small, which selects the ID class triads. For a sufficiently red spectrum, this permits discarding of the small $n_{\boldsymbol {p}} n_{\boldsymbol {p} 1}$ ($n_{\boldsymbol {p}} n_{\boldsymbol {p} 2}$, respectively) terms. We then rewrite (2.16) as

(2.21)\begin{equation} \frac{\partial n_{{\boldsymbol{p}}}}{\partial t} = \int {\rm d} {\boldsymbol{q}} \left( {\mathcal{B}}(\kern1pt {\boldsymbol{p}})-{\mathcal{B}}(\kern1pt {\boldsymbol{p}}+{\boldsymbol{q}})\right)\simeq{-}\int {\rm d} {\boldsymbol{q}}\left({\boldsymbol{q}}\,\boldsymbol{\cdot}\,\frac{\partial}{\partial {\boldsymbol{p}}}\right) {\mathcal{B}}(\kern1pt {\boldsymbol{p}}), \end{equation}

where we introduced

(2.22)\begin{equation} {\mathcal{B}}(\kern1pt {\boldsymbol{p}}) = 8{\rm \pi} \int |V_{{\boldsymbol{p}}_1,{\boldsymbol{q}}}^{{\boldsymbol{p}}}|^2 \, n_q(n_{\boldsymbol {p}1}-n_{\boldsymbol{p}}) \, \delta_{{{\boldsymbol{p}} - {\boldsymbol{p}}_1-{\boldsymbol{q}}}} \, {\mathcal{L}}({\sigma_{{\boldsymbol{p}}} -\sigma_{{{\boldsymbol{p}}_1}}-\sigma_{{{\boldsymbol{q}}}}}) \,{\rm d} {\boldsymbol{p}}_{1}, \end{equation}

and expanded the difference $( {\mathcal {B}}({\boldsymbol {p}})-{\mathcal {B}}({\boldsymbol {p}}+{\boldsymbol {q}}))$ in a Taylor series using ${\boldsymbol {q}}$ to represent the small difference in wavenumber between the two high-frequency waves. Expanding the difference $n_{{\boldsymbol {p}}_1}-n_{\boldsymbol {p}}$ for small ${\boldsymbol {q}}$ gives

(2.23)\begin{equation} {\mathcal{B}}(\kern1pt {\boldsymbol{p}}) \simeq{-}8{\rm \pi} \left({\boldsymbol{q}}\,\boldsymbol{\cdot}\,\frac{\partial n_{\boldsymbol{p}}}{\partial {\boldsymbol{p}}}\right) \int |V_{{\boldsymbol{p}}_1,{\boldsymbol{q}}}^{{\boldsymbol{p}}}|^2 n_q \delta_{{{\boldsymbol{p}} - {\boldsymbol{p}}_1-{\boldsymbol{q}}}} {\mathcal{L}}({\sigma_{{\boldsymbol{p}}} -\sigma_{{{\boldsymbol{p}}_1}}-\sigma_{{{\boldsymbol{q}}}}}) \,{\rm d} {\boldsymbol{p}}_{1}. \end{equation}

Combining (2.21) with (2.23) we obtain

(2.24)\begin{align} \frac{\partial n_{{\boldsymbol{p}}}}{\partial t} = 8{\rm \pi} \int {\rm d} {\boldsymbol{q}}\left({\boldsymbol{q}}\,\boldsymbol{\cdot}\,\frac{\partial}{\partial {\boldsymbol{p}}}\right) \left({\boldsymbol{q}}\,\boldsymbol{\cdot}\,\frac{\partial n_{\boldsymbol{p}}}{\partial {\boldsymbol{p}}}\right) \left[ \int |V_{{\boldsymbol{p}}_1,{\boldsymbol{q}}}^{{\boldsymbol{p}}}|^2 n_q \delta_{{{\boldsymbol{p}} - {\boldsymbol{p}}_1-{\boldsymbol{q}}}} {\mathcal{L}}({\sigma_{{\boldsymbol{p}}} -\sigma_{{{\boldsymbol{p}}_1}}-\sigma_{{{\boldsymbol{q}}}}}) \,{\rm d} {\boldsymbol{p}}_{1}\right], \end{align}

or

(2.25)\begin{equation} \left.\begin{gathered} \frac{\partial n_{{\boldsymbol{p}}}}{\partial t} = \frac{\partial}{\partial p_i} D_{ij}(\kern1pt {\boldsymbol{p}},{\boldsymbol{q}}) \frac{\partial}{\partial p_j} n_{{\boldsymbol{p}}},\\ D_{ij}(\kern1pt {\boldsymbol{p}},{\boldsymbol{q}}) = 8{\rm \pi} \int {\rm d} {\boldsymbol{q}} \left(q_i q_j\right) |V_{{\boldsymbol{p}}_1,{\boldsymbol{q}}}^{{\boldsymbol{p}}}|^2 n_q \delta_{{{\boldsymbol{p}} - {\boldsymbol{p}}_1-{\boldsymbol{q}}}} {\mathcal{L}}({\sigma_{{\boldsymbol{p}}} -\sigma_{{{\boldsymbol{p}}_1}}-\sigma_{{{\boldsymbol{q}}}}}) \,{\rm d} {\boldsymbol{p}}_{1}. \end{gathered}\right\}\end{equation}

This is a Fokker–Planck diffusion equation describing the diffusion of wave action in the system dominated by non-local-in-wavenumber interactions. Comparison between the Fokker–Planck equation (2.25) obtained here, and the similar Fokker–Planck equation (obtained using Wentzel–Kramers–Brillouin (WKB) theory (§ 3.2.4 below) will lead to critical insights into the spectral energy transfers in internal wave systems and ultimately to the parametrization of the energy supply to internal wave breaking processes.

Extended versions of the Garrett and Munk model (Appendix A) can then be inserted into (2.25) to arrive at families of stationary states (e.g. McComas & Müller Reference McComas and Müller1981b; Polzin & Lvov Reference Polzin and Lvov2011; Dematteis et al. Reference Dematteis, Polzin and Lvov2022) and estimate downscale transport. We will require transport estimates the vertical–vertical component of the diffusivity tensor, $D_{33}$. For the Garrett and Munk model (GM76),

(2.26)\begin{equation} D_{33} = \frac{2}{\rm \pi} \frac{k m^2 e_0 m_{{\ast}}}{N} . \end{equation}

Here, $k$ is horizontal wavenumber magnitude, $e_0$ is the total energy, $m_{\ast }$ is a bandwidth parameter and $N$ is buoyancy frequency.

3.1.1. Reynolds decomposition and Hamiltonian structure

To study interactions between long and short waves we start at (2.8a,b) and make a Reynolds decomposition in wave amplitude

(3.1a–c)\begin{equation} \varPi \rightarrow \varPi_{0}+ \varPi+ {\rm \pi}', \quad \phi=\varPhi+\phi' , \quad \psi = \varPsi +\psi' . \end{equation}

Here, the large-amplitude waves are represented with $\varPi, \varPhi, \varPsi$ and small-amplitude waves are given by $\phi '$, ${\rm \pi} '$ and $\psi '$. Given the potentials $\varPhi$ and $\phi$, the corresponding velocities are

(3.2a,b)\begin{equation} {\mathcal{U}} = \boldsymbol{\nabla} \varPhi + \boldsymbol{\nabla}^{{\perp}} \varPsi, \quad u'=\boldsymbol{\nabla} \phi' + \boldsymbol{\nabla}^{{\perp}} \psi'. \end{equation}

To simplify the presentation, we will utilize the non-rotating approximation ($\kern0.7pt f=0$) in which $(\boldsymbol {\nabla }^{\perp } \varPsi, \boldsymbol {\nabla }^{\perp } \psi ) \rightarrow 0$. The case of rotating ocean $f\ne 0$ is presented in Appendix.

We substitute the Reynolds decomposition (3.1a–c) into the equations of motion (2.1). We then make an assumption that the large scales $\varPi$ and $\varPhi$ exactly satisfy the primitive equations of motion. In other words, gradient operator applied to the inertial waves with no horizontal structure will return zero, thus these inertial waves are exact analytical solutions of the primitive Boussinesq equations (2.1). For the application of oceanic internal waves, this assumption is realized if the large-scale waves are a collection of inertial waves having frequency $\sigma =f$: inertial waves do not have horizontal structure and thus a superposition of inertial waves exactly satisfies the primitive equations of motion. We then subtract equations for the large-amplitude waves. The result is given by

(3.3) \begin{equation} \left.\begin{gathered} \dot{\rm \pi}' + \boldsymbol{\nabla}\,\boldsymbol{\cdot}\,((\varPi_{0}+ \varPi+{\rm \pi}') \boldsymbol{\nabla} \varPhi)+ \boldsymbol{\nabla}({\rm \pi}'\varPhi)=0, \\ \dot{\phi}'+\frac{|\boldsymbol{\nabla}\phi'|^2}{2}+\boldsymbol{\nabla}\phi' \,\boldsymbol{\cdot}\, \boldsymbol{\nabla}\varPhi + \frac{g}{\rho_0^2} \int\int {\rm d}\rho' \,{\rm d}\rho^{\prime\prime} \,{\rm \pi}'=0. \end{gathered}\right\}\end{equation}

In these equations, $\varPi$ and $\varPhi$ are given time–space-dependent functions representing the large-amplitude waves.

These equations are also Hamilton's equations,

(3.4a,b)\begin{equation} \frac{\partial {\rm \pi}'}{\partial t} = \frac{\delta {\mathcal{H}}}{\delta \phi'} , \quad \frac{\partial \phi'}{\partial t} ={-}\frac{\delta {\mathcal{H}}}{\delta {\rm \pi}'} , \end{equation}

with the time-dependent Hamiltonian given by (this Hamiltonian may be obtained by substituting (3.1a–c) to (2.9))

(3.5) \begin{equation} {\mathcal{H}} =\frac{1}{2} \int {\rm d} {\boldsymbol{r}} \Big( \Big( \varPi_{0}+\varPi+{\rm \pi}' \Big) \Big|\boldsymbol{\nabla} \phi' \Big|^2 + 2{\rm \pi}'\boldsymbol{\nabla}\phi'\,\boldsymbol{\cdot}\,\boldsymbol{\nabla}\varPhi - g \Big|\int^{\rho} {\rm d}\hat{\rho} \,\frac{{\rm \pi}'}{\hat{\rho}} \Big|^2 \Big) . \end{equation}

There are two types of terms here – those that will ultimately describe a sea of interacting small-scale waves and those that will describe the influence of large-amplitude large-scale waves on the small-scale waves. In our previous efforts (Lvov & Tabak Reference Lvov and Tabak2001, Reference Lvov and Tabak2004; Lvov et al. Reference Lvov, Polzin, Tabak and Yokoyama2010; Polzin & Lvov Reference Polzin and Lvov2011, Reference Polzin and Lvov2017) these terms are comingled. Comparing (2.9) and (3.5), we see that (3.5) contains additional terms $\varPi |\boldsymbol {\nabla } \phi ' |^2$ and ${\rm \pi} '\boldsymbol {\nabla }\phi '\boldsymbol {\cdot }\boldsymbol {\nabla }\varPhi$ that are explicit representations of what will be scale-separated interactions. The term ${\rm \pi} '\boldsymbol {\nabla }\phi '\boldsymbol {\cdot }\boldsymbol {\nabla }\varPhi$ describes the advection of the small-scale internal field by the given large-scale large-amplitude field. The term $\varPi |\boldsymbol {\nabla } \phi ' |^2$ represents a coupling of small scales to large through changes in the stratification by the large-scale wave. The term ${\rm \pi} ' |\boldsymbol {\nabla } \phi ' |^2$ will represent interactions local in wavenumber.

We now express the space-dependent variables ${\rm \pi} ^{\prime }({\boldsymbol {r}}),\varPi ({\boldsymbol {r}}),\phi ^{\prime }({\boldsymbol {r}})$ and $\varPhi ({\boldsymbol {r}})$ in terms of their Fourier images ${\rm \pi} ^{\prime }({\boldsymbol {p}}),\varPi ({\boldsymbol {p}}),\phi ^{\prime }({\boldsymbol {p}})$ via (2.10), make the Boussinesq approximation ${\varPi }/{\rho }\simeq {\varPi }/{\rho _0}$ and use $\int \textrm {d} {\boldsymbol {p}} \,\textrm {e}^{\textrm {i} {\boldsymbol {p}}\boldsymbol {\cdot } {\boldsymbol {r}}} =(2{\rm \pi} )^3 \delta ({\boldsymbol {p}})$ to obtain

(3.6) \begin{gather} \left.\begin{gathered} {\mathcal{H}} = {\mathcal{H}}_{linear} +{\mathcal{H}}_{nonlinear},\\ {\mathcal{H}}_{linear} = \frac{1}{2}\int {\rm d} {\boldsymbol{p}} \Big( \varPi_{0} |{\boldsymbol{k}}|^2 |\phi'_{{\boldsymbol{p}}}|^2- \frac{g}{\rho_0^2}\frac{| \phi'_{{\boldsymbol{p}}} |^2}{m^2}\Big), \\ {\mathcal{H}}_{nonlinear} = {\mathcal{H}}_{local}+ {\mathcal{H}}_{sweeping} + {\mathcal{H}}_{density}, \end{gathered}\right\} \end{gather}

(3.7)\begin{gather} \left.\begin{gathered} {\mathcal{H}}_{\textit{local}} ={-} {\frac{\textbf{1}}{\textbf{2}(\textbf{2}\boldsymbol{\rm \pi})^{{\textbf{3}}/{\textbf{2}}}}} \int {\rm d} {\boldsymbol{p}}_1 \,{\rm d} {\boldsymbol{p}}_2 \,{\rm d} {\boldsymbol{p}}_3 \delta(\kern1pt {\boldsymbol{p}}_1+{\boldsymbol{p}}_2+{\boldsymbol{p}}_3) {\boldsymbol{k}}_2 \,\boldsymbol{\cdot}\, {\boldsymbol{k}}_3 {\rm \pi}'_{{\boldsymbol{p}}_1} \phi'_{{\boldsymbol{p}}_2} \phi'_{{\boldsymbol{p}}_3}, \\ {\mathcal{H}}_{sweeping} ={-} {\frac{\textbf{1}}{\textbf{2}(\textbf{2}\boldsymbol{\rm \pi})^{{\textbf{3}}/{\textbf{2}}}}} \int {\rm d} {\boldsymbol{p}}_1\,{\rm d} {\boldsymbol{p}}_2 \,{\rm d} {\boldsymbol{p}}_3 \delta(\kern1pt {\boldsymbol{p}}_1+{\boldsymbol{p}}_2+{\boldsymbol{p}}_3) 2 {\boldsymbol{k}}_1\,\boldsymbol{\cdot}\,{\boldsymbol{k}}_2 \varPhi_{{\boldsymbol{p}}_1} \phi'_{{\boldsymbol{p}}_2} {\rm \pi}'_{{\boldsymbol{p}}_3}, \\ {\mathcal{H}}_{density} ={-}{\frac{\textbf{1}}{\textbf{2}(\textbf{2}\boldsymbol{\rm \pi})^{{\textbf{3}}/{\textbf{2}}}}} \int {\rm d} {\boldsymbol{p}}_1 \,{\rm d} {\boldsymbol{p}}_2 \,{\rm d} {\boldsymbol{p}}_3 \delta(\kern1pt {\boldsymbol{p}}_1+{\boldsymbol{p}}_2+{\boldsymbol{p}}_3) {\boldsymbol{k}}_2\,\boldsymbol{\cdot}\,{\boldsymbol{k}}_3 \varPi_{{\boldsymbol{p}}_1} \phi'_{{\boldsymbol{p}}_2} \phi'_{{\boldsymbol{p}}_3}. \end{gathered}\right\}\end{gather}

The Hamiltonian ${\mathcal {H}}_{nonlinear}$ is the sum of three terms:

1. (i) ${\mathcal {H}}_{local}$ represents small amplitudes interacting with small amplitudes. We will call this term ‘local’ interactions in anticipation of making a scale separation between large-amplitude large-scale and small-amplitude small-scale waves in §§ 3.1.2 and 3.1.3;

2. (ii) ${\mathcal {H}}_{density}$ is the term that describes the variations of stratification that small-amplitude waves experience due to the compression and rarefication of isopycnals associated with the large-amplitude waves. We will refer to this term as a density term;

3. (iii) ${\mathcal {H}}_{sweeping}$ is the term that describes the advection (sweeping) of small-amplitude waves by large-amplitude waves. In the future, we refer to this term as a sweeping term.

The main focus of this manuscript is to investigate how the density and sweeping terms affect the overall spectral energy density.

3.1.2. Sweeping Hamiltonian

Following the traditional wave turbulence approach, we make a transformation to the wave-action variables that represent wave amplitude and phase

(3.8a,b)\begin{equation} \phi'_{\boldsymbol{p}}=\frac{{\rm i}N \sqrt{\sigma_{\boldsymbol{p}}}}{\sqrt{2 g} |{\boldsymbol{k}}|} (a_{\boldsymbol{p}}-a^*_{-{\boldsymbol{p}}}) , \quad {\rm \pi}'_{\boldsymbol{p}} =\frac{\sqrt{g} |{\boldsymbol{k}}|}{\sqrt{2\sigma_{\boldsymbol{p}}}N} (a _{\boldsymbol{p}}+ a^*_{-{\boldsymbol{p}}}) . \end{equation}

We substitute (3.8a,b) into ${\mathcal {H}}_{sweeping}$ of (3.7), and obtain

(3.9) \begin{align} {\mathcal{H}}_{sweeping}&={-} {\frac{\textbf{1}}{\textbf{2}(\textbf{2}\boldsymbol{\rm \pi})^{{\textbf{3}}/{\textbf{2}}}}} \int {\rm d} {\boldsymbol{p}}_1 \,{\rm d} {\boldsymbol{p}}_2 \,{\rm d} {\boldsymbol{p}}_3 \delta(\kern1pt {\boldsymbol{p}}_1+{\boldsymbol{p}}_2+{\boldsymbol{p}}_3) 2 {\boldsymbol{k}}_1\,\boldsymbol{\cdot}\,{\boldsymbol{k}}_2 \varPhi_{{\boldsymbol{p}}_1} \frac{{\rm i}N \sqrt{\sigma_{{\boldsymbol{p}}_2}}}{\sqrt{2 g}|{\boldsymbol{k}}_2|} \frac{\sqrt{g} |{\boldsymbol{k}}_3|}{\sqrt{2\sigma_{{\boldsymbol{p}}_3}}N}\nonumber\\ &\quad \times\big(a_{{\boldsymbol{p}}_2}-a^*_{-{\boldsymbol{p}}_2}\big) \big(a _{{\boldsymbol{p}}_3}+ a^*_{-{\boldsymbol{p}}_3}\big) . \end{align}

The next step is algebraically trivial but conceptually fundamental. We invoke an extreme scale-separated limit in which the two small-amplitude waves have similar frequency and horizontal wavenumber magnitude. No condition is required on the vertical wavenumber. This conditioning retains both the ID and Bragg scattering (ES) branches of the resonant manifold; see figure 1.

In this scale-separated limit of the internal wave problem, $\sigma _{{\boldsymbol {p}}_2} \cong \sigma _{{\boldsymbol {p}}_3}$ and $|{\boldsymbol {k}}_2| \cong |{\boldsymbol {k}}_3|$. Beyond the obvious algebraic simplifications, upon expanding the brackets, we find terms of the type $a_{{\boldsymbol {p}}_2}a_{{\boldsymbol {p}}_3}$, $a^*_{-{\boldsymbol {p}}_2}a^*_{-{\boldsymbol {p}}_3}$ and $a_{{\boldsymbol {p}}_2}a^*_{-{\boldsymbol {p}}_3}$, $a^*_{-{\boldsymbol {p}}_2}a_{-{\boldsymbol {p}}_3}$. In what follows, we neglect the $a_{{\boldsymbol {p}}_2}a_{{\boldsymbol {p}}_3}$ and $a^*_{-{\boldsymbol {p}}_2}a^*_{-{\boldsymbol {p}}_3}$ terms and retain the $a_{{\boldsymbol {p}}_2}a^*_{-{\boldsymbol {p}}_3}$, $a^*_{-{\boldsymbol {p}}_2}a_{-{\boldsymbol {p}}_3}$ terms, since the former terms are non-resonant, while the latter may be in the resonance for some wavenumbers. The discarded terms lead to a process when one lower-frequency wave decays into two high-frequency waves and thus the frequencies do not sum to zero. Such decay is a non-resonant process, so we can remove these terms at the onset.

After relabelling subscripts, in which $2\to 1$ and $3\to 2$, we obtain

(3.10)\begin{equation} \left.\begin{gathered} {\mathcal{H}}_{sweeping} =\int {\rm d} {\boldsymbol{p}}_1 \,{\rm d} {\boldsymbol{p}}_2 \,A_{sweeping}(\kern1pt {\boldsymbol{p}}_1, {\boldsymbol{p}}_2) a_{{\boldsymbol{p}}_1}a^*_{{\boldsymbol{p}}_2},\\ {\rm with}\\ A_{sweeping}(\kern1pt {\boldsymbol{p}}_1, {\boldsymbol{p}}_2) ={-}\frac{1}{2}{\rm i} \left({{\boldsymbol{k}}_1-{\boldsymbol{k}}_2}\right){\cdot}\left( {\boldsymbol{k}}_1+{\boldsymbol{k}}_2\right) \varPhi_{{\boldsymbol{p}}_1-{\boldsymbol{p}}_2} . \end{gathered}\right\}\end{equation}

In the limit of large vertical background scales, (3.10) describes a quasi-coherent translation of small-scale small-amplitude waves by the large-scale background. At $1/2$ the vertical scale of ${\boldsymbol {p}}_1$ and ${\boldsymbol {p}}_2$, (3.10) describes a Bragg scattering process. This is distinct from ‘local’ interactions as the large-amplitude wave has a much larger horizontal scale.

3.1.3. Density Hamiltonian

We now can repeat the same steps for the density Hamiltonian. We substitute (3.8a,b) into the ${\mathcal {H}}_{density}$ of (3.7), use $\sigma _{{\boldsymbol {p}}_2} \cong \sigma _{{\boldsymbol {p}}_3}$ and $|{\boldsymbol {k}}_2| \cong |{\boldsymbol {k}}_3|$, relabel subscripts and obtain

(3.11)\begin{equation} \left.\begin{gathered} {\mathcal{H}}_{density} =\int {\rm d} {\boldsymbol{p}}_1 \,{\rm d} {\boldsymbol{p}}_2 \,A_{density}(\kern1pt {\boldsymbol{p}}_1, {\boldsymbol{p}}_2) a_{{\boldsymbol{p}}_1}a^*_{{\boldsymbol{p}}_2}, \\ {\rm with}\\ A_{density}(\kern1pt {\boldsymbol{p}}_1, {\boldsymbol{p}}_2) = \frac{1}{2 \varPi_{0}} \varPi_{{\boldsymbol{p}}_2-{\boldsymbol{p}}_1} \sqrt{ \sigma_{{\boldsymbol{p}}_1}\sigma_{{\boldsymbol{p}}_2}} . \end{gathered}\right\} \end{equation}

In the limit of large vertical background scales, (3.11) describes the modulation of the background stratification. At $1/2$ the vertical scale of ${\boldsymbol {p}}_1$ and ${\boldsymbol {p}}_2$, (3.11) describes a Bragg scattering process. This is distinct from ‘local’ interactions as the large-amplitude wave has a much larger horizontal scale.

3.1.4. Quadratic Hamiltonian for inhomogeneous wave turbulence

Let us neglect the local interaction term ${\mathcal {H}}_{local}$ in Hamiltonian (3.7). Then the Hamiltonian (3.7) of small-amplitude small horizontal scale internal waves $a_{\boldsymbol {p}}$ superimposed into a field of large-amplitude large horizontal scale internal waves given by space–time-dependent $\varPi$ and $\varPhi$ are given by a form that is quadratic in the small-scale variables $a_{{\boldsymbol {p}}}$

(3.12)\begin{equation} \left.\begin{gathered} {\mathcal{H}} =\int {\rm d} {\boldsymbol{p}}_1 \,{\rm d} {\boldsymbol{p}}_2 \,A(\kern1pt {\boldsymbol{p}}_1, {\boldsymbol{p}}_2) a_{{\boldsymbol{p}}_1}a^*_{{\boldsymbol{p}}_2}, \\ {\rm with} \\ A(\kern1pt {\boldsymbol{p}}_1,{\boldsymbol{p}}_2)= \sigma_{{\boldsymbol{p}}} \delta(\kern1pt {\boldsymbol{p}}_1-{\boldsymbol{p}}_2) + A_{sweeping}(\kern1pt {\boldsymbol{p}}_1, {\boldsymbol{p}}_2)+ A_{density }(\kern1pt {\boldsymbol{p}}_1, {\boldsymbol{p}}_2), \end{gathered}\right\} \end{equation}

where $A_{sweeping}({\boldsymbol {p}}_1, {\boldsymbol {p}}_2)$ and $A_{density }({\boldsymbol {p}}_1, {\boldsymbol {p}}_2)$ are given in (3.10) and (3.11). The $A({\boldsymbol {p}}_1,{\boldsymbol {p}}_2)$ are time dependent and depend upon the phases of the external field.

3.2.1. The wave packet transport equation

In the spatially homogeneous wave turbulence of § 2.2 we used wave-action density $n_{\boldsymbol {p}}$ and a linear dispersion relation $\sigma _{\boldsymbol {p}}$, both being a function of a wavenumber ${\boldsymbol {p}}$. In the case when there is a slowly varying large-scale background, i.e. a system where spatial inhomogeneity is present, the properties of wave action and the dispersion relation will depend on the position in space. Then it makes sense to introduce an additional parameter, a position vector ${\boldsymbol {r}}$ in wave action and linear dispersion relation. The theory for this spatial dependence is developed in Gershgorin, Lvov & Nazarenko (Reference Gershgorin, Lvov and Nazarenko2009) using a Gabor transform to represent the envelope structure describing the spatial localization and carrier frequency. The leading-order balance in Gershgorin et al. (Reference Gershgorin, Lvov and Nazarenko2009) leads to action and phase conservation along ray paths. The balance on the envelope scale implies phase modulation, for which we direct the reader to the appendix of Cohen & Lee (Reference Cohen and Lee1990) for clarity. Associated with the envelope structure is a residual circulation (Bühler & McIntyre Reference Bühler and McIntyre2005). The potential for the wave packet to interact with its envelope structure is possible (Bühler & McIntyre Reference Bühler and McIntyre2005; Dosser & Sutherland Reference Dosser and Sutherland2011) but would require a modification of the uniform potential vorticity statement (2.7). It is at this stage that one might also want to consider the potential for nonlinear wave steepening effects associated with ${\mathcal {H}}_{local}$ to counter the dispersion of ray characteristics as a precursor to the description of the solitary wave dynamics.

The familiar statement of action conservation that we are after is obtained in Gershgorin et al. (Reference Gershgorin, Lvov and Nazarenko2009) by assuming the scale of the envelope structure is large in comparison with the inverse wavenumber of the small-scale wave, i.e. the wave packet contains many oscillations. The result required here can be obtained more simply by using a Wigner transform alone to define the space–time-dependent wave-action spectral density

(3.13)\begin{equation} n_{{\boldsymbol{p}},{\boldsymbol{r}}}\equiv \int {\rm e}^{{\rm i}{\boldsymbol{q}}\,\boldsymbol{\cdot}\,{\boldsymbol{r}}} \langle a_{{\boldsymbol{p}}+{\boldsymbol{q}}/2} a_{ {\boldsymbol{p}}-{\boldsymbol{q}}/2}^*\rangle \,{\rm d} {\boldsymbol{q}} , \end{equation}

in which the transform variable ${\boldsymbol {q}}$ is the difference between two large wavenumbers ${\boldsymbol {p}}_1$ and ${\boldsymbol {p}}_2$, ${\boldsymbol {q}}={\boldsymbol {p}}_1-{\boldsymbol {p}}_2$ and the field variables $a$ are evaluated at ${\boldsymbol {p}} \pm {\boldsymbol {q}}/2$. We similarly introduce the space-dependent intrinsic frequency $\omega _{{\boldsymbol {p}},{\boldsymbol {r}}}$

(3.14)\begin{equation} \omega_{{\boldsymbol{p}},{\boldsymbol{r}}}=\int {\rm e}^{{\rm i} {\boldsymbol{r}}\,\boldsymbol{\cdot}\,{\boldsymbol{q}}}A({{\boldsymbol{p}}+{\boldsymbol{q}}/2,{\boldsymbol{p}}-{\boldsymbol{q}}/2}) \,{\rm d}{\boldsymbol{q}}.\end{equation}

This is a generalization of our ‘traditional’ wave action (2.15) which allows for slow variations of the background. It will incorporate the large vertical scale contributions of ${\mathcal {H}}_{sweeping}$ and ${\mathcal {H}}_{density}$.

To derive the transport equation, we take definition (3.13), differentiate it with respect to time, use Hamilton's equation of motion (2.13), with the Hamiltonian (3.12). The spatial dependence is then made explicit by representing $A({\boldsymbol {p}}_1,{\boldsymbol {p}}_2)$ and $a({\boldsymbol {p}}_1)a^{\ast }({\boldsymbol {p}}_2)$ with their corresponding Fourier transforms. After an inspired change of variables (Lvov & Rubenchik Reference Lvov and Rubenchik1977; Gershgorin et al. Reference Gershgorin, Lvov and Nazarenko2009) that maps the ID portion of the resonant manifold, one obtains an intermediary expression

(3.15) \begin{align} {\rm i} \frac{\partial n_{{\boldsymbol{p}},{\boldsymbol{r}}}}{\partial t} &= \int \frac{{\rm d}{\boldsymbol{r}}^{\prime}\,{\rm d}{\boldsymbol{r}}^{\prime\prime}}{(2{\rm \pi})^6} {\rm d}q^{\prime}\,{\rm d}q^{\prime\prime} \exp({{\rm i}q^{\prime}\,\boldsymbol{\cdot}\,({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime}) + {\rm i}q^{\prime\prime}\,\boldsymbol{\cdot}\, ({\boldsymbol{r}}-{\boldsymbol{r}}^{\prime\prime}) }) \nonumber\\ &\quad \times[\omega_{{\boldsymbol{p}}+{q^{\prime\prime}}/{2}}({\boldsymbol{r}}^{\prime}) n_{{\boldsymbol{p}}-{q^{\prime}}/{2}} ({\boldsymbol{r}}^{\prime\prime}) - \omega_{{\boldsymbol{p}}-{q^{\prime\prime}}/{2}}({\boldsymbol{r}}^{\prime}) n_{{\boldsymbol{p}}+{q^{\prime}}/{2}}({\boldsymbol{r}}^{\prime\prime}) ]. \end{align}

After expanding $\omega$ and $n$ in Taylor series with respect to ${\boldsymbol {p}}$ and truncating higher-order terms, one ultimately arrives at the action balance

(3.16)\begin{equation} \frac{\partial n_{{\boldsymbol{p}},{\boldsymbol{r}}}}{\partial t}+ \boldsymbol{\nabla}_{{\boldsymbol{p}}} \sigma_{{\boldsymbol{p}},{\boldsymbol{r}}}\,\boldsymbol{\cdot}\, \boldsymbol{\nabla}_{{\boldsymbol{r}}}n_{{\boldsymbol{p}},{\boldsymbol{r}}}-\boldsymbol{\nabla}_{{\boldsymbol{r}}}\sigma_{{\boldsymbol{p}},{\boldsymbol{r}}}\,\boldsymbol{\cdot}\,\boldsymbol{\nabla}_{{\boldsymbol{p}}} n_{{\boldsymbol{p}},{\boldsymbol{r}}}=0, \end{equation}

or alternately

(3.17) \begin{equation} \frac{\partial n_{{\boldsymbol{p}},{\boldsymbol{r}}}}{\partial t}+ \boldsymbol{\nabla}_{{\boldsymbol{r}}} \,\boldsymbol{\cdot}\, [n_{{\boldsymbol{p}},{\boldsymbol{r}}} \boldsymbol{\nabla}_{{\boldsymbol{p}}} \sigma_{{\boldsymbol{p}},{\boldsymbol{r}}} ] - \boldsymbol{\nabla}_{{\boldsymbol{p}}} \,\boldsymbol{\cdot}\, [n_{{\boldsymbol{p}},{\boldsymbol{r}}} \boldsymbol{\nabla}_{{\boldsymbol{r}}} \sigma_{{\boldsymbol{p}},{\boldsymbol{r}}} ] = 0. \end{equation}

Integration of (3.17) over wavenumber provides a connection to the space–time variational formulations found in Witham (Reference Witham1974, Chapters 11.7 and 14). The Bragg scattering process residing in the scale-separated Hamiltonian (3.10) has been eliminated by using the Wigner transform and Taylor series expansion that serendipitously exploits the symmetries associated with the ID resonance.

3.2.2. Application for internal waves

We now take the expression for $A_{{\boldsymbol {p}}_1, {\boldsymbol {p}}_2}$ from (3.12) and substitute it into (3.14) for the space-dependent linear dispersion relation. The result is given by

(3.18)\begin{equation} \omega_{{\boldsymbol{p}},{\boldsymbol{r}}}= \sigma_{{\boldsymbol{p}}} - {\boldsymbol{k}}\,\boldsymbol{\cdot}\, {\mathcal{U}}({\boldsymbol{r}}) + \sigma_{{\boldsymbol{p}}}\frac{\varPi({\boldsymbol{r}})}{2\varPi_0}, \end{equation}

where ${\boldsymbol {p}}=({\boldsymbol {k}},m)$ is the wavevector, ${\mathcal {U}}({\boldsymbol {r}})$ is the time-dependent horizontal velocity of the external large-scale wave field and $\varPi ({\boldsymbol {r}})$ is the time-dependent stratification of the external wave field. Here, the $\sigma _{{\boldsymbol {p}}}$ term was produced by the term proportional to $\delta ({\boldsymbol {p}}_1-{\boldsymbol {p}}_2)$, ${\boldsymbol {k}}\boldsymbol {\cdot }{\mathcal {U}}$ comes from the sweeping term in the Hamiltonian and $\sigma _{{\boldsymbol {p}}} \varPi / 2 \varPi _0$ comes from density term in the Hamiltonian. The frequency $\sigma _{{\boldsymbol {p}}}$ is given by (2.12) using $f=0$. For a high-frequency wave in a rotating ocean, $\mathcal {U}$ is replaced by $\boldsymbol {\nabla } \varPhi + \boldsymbol {\nabla }^{\perp } \varPsi$. The derivation for $f \ne 0$ appears in the Appendix.

3.2.3. The ray-path wave packet transport equation

The transport equation (3.16) or its alternative formulation (3.17) with the space–time-dependent linear dispersion relationship (3.18) is a fundamental result expressing wave-action conservation that provides the basis for the analyses to follow.

The action balance (3.16) can be simply solved by the method of characteristics. The characteristics of (3.16), also called rays, are defined by

(3.19a,b)\begin{equation} \dot{{\boldsymbol{r}}}(t) \equiv \boldsymbol{\nabla}_{{\boldsymbol{p}}} \sigma_{{\boldsymbol{p}},{\boldsymbol{r}}} , \quad \dot{{\boldsymbol{p}}}(t) \equiv{-} \boldsymbol{\nabla}_{\boldsymbol{r}}\sigma_{{\boldsymbol{p}},{\boldsymbol{r}}} . \end{equation}

Equations (3.19a,b) imply that wave-action spectral density is conserved along these characteristics

(3.20)\begin{equation} \frac{\partial n_{{\boldsymbol{p}},{\boldsymbol{r}}}}{\partial t}+ \dot{{\boldsymbol{r}}}\,\boldsymbol{\cdot}\,\boldsymbol{\nabla}_{{\boldsymbol{r}}}n_{{\boldsymbol{p}},{\boldsymbol{r}}}+\dot{{\boldsymbol{p}}} \,\boldsymbol{\cdot}\,\boldsymbol{\nabla}_{{\boldsymbol{p}}}n_{{\boldsymbol{p}},{\boldsymbol{r}}}=0. \end{equation}

This is the classical representation for the conservation of action spectral density along ray trajectories. A derivation for a general Hamiltonian set is presented in Gershgorin et al. (Reference Gershgorin, Lvov and Nazarenko2009), here, the derivation is specifically for internal waves in the scale-separated limit. Integration over wavenumber provides the space–time result presented in Witham (Reference Witham1974). This result often appears as an analogy to Liouville's theorem for the conservation of the phase volume of particles without justification. Note that it is the balance of first-order terms in a Taylor series expansion of (3.15).

3.2.4. An ensemble path transport equation

In what follows we derive a combined advective–diffusive transport equation for the action balance (3.16) by generalizing an approach found in Nazarenko et al. (Reference Nazarenko, Newell and Galtier2001). There are strong parallels here to the discussion of Taylor (Reference Taylor1921) that appears in the appendix of McComas & Bretherton (Reference McComas and Bretherton1977): the illuminating analogy with particle dispersion in Taylor (Reference Taylor1921) is to substitute wavenumber ${\boldsymbol {p}}$ for the Lagrangian particle position ${\boldsymbol {r}}$ and obtain a quantitative approach for discussing the migration and dispersion of wave packets in the spectral domain following ray trajectories. Having said this, it should be intuitively obvious that, since particle dispersion in ${\boldsymbol {r}}$ has little to do with resonance, neither should the issue of dispersion of ${\boldsymbol {p}}$ in phase space be intrinsically tied to resonance! Yet, an emphasis on the resonant paradigm is the interpretive context pursued in McComas & Bretherton (Reference McComas and Bretherton1977) and Nazarenko et al. (Reference Nazarenko, Newell and Galtier2001). A second key departure is that our motivation stems from the fact that the GM76 spectrum is a no-flux solution to the Fokker–Planck equation (2.25). We acknowledge the oceanographic literature in this regard and so are focused upon the issue of a mean drift of wave packets to high wavenumber that is accessible in ray-tracing simulations (e.g. Henyey et al. Reference Henyey, Wright and Flatté1986) but not explicitly represented in (2.25). To underscore the distinction, the existence of a mean drift implies relative dispersion (e.g. Bennett Reference Bennett1984) rather than the issue of absolute dispersion addressed in Taylor (Reference Taylor1921).

The very first step in our derivation is to specify a wave packet ensemble mean drift and departures from the mean drift. This decomposition is an improvement on the arguments presented in McComas & Bretherton (Reference McComas and Bretherton1977), Nazarenko et al. (Reference Nazarenko, Newell and Galtier2001) and Lanchon & Cortet (Reference Lanchon and Cortet2023). Our limited knowledge of the ray literature does not permit us to comment on the originality of our interpretation. It is, however, crucial in understanding intrinsic differences between the kinetic equation and ray theory. Given the nature of our understanding, the issue assuredly carries over to other physical problems such as the interaction of near-inertial waves with lower-frequency flows, (e.g. Young & Jelloul Reference Young and Jelloul1997; Kafiabad et al. Reference Kafiabad, Savva, Miles and Vanneste2019; Dong, Bühler & Smith Reference Dong, Bühler and Smith2020) surface gravity wave interaction with lower-frequency flows (Villas Bôas & Young Reference Villas Bôas and Young2020) and Rossby wave–Rossby wave interactions (Nazarenko Reference Nazarenko2011).

We represent the wave action of a single wave packet as a sum of a spatially homogeneous part $\bar {n}$ and small ‘wiggles’ $\tilde {n}$

(3.21a–c)\begin{equation} n_{{\boldsymbol{p}},{\boldsymbol{r}}} = \bar{n}_{{\boldsymbol{p}}} + \tilde{n}_{{\boldsymbol{p}},{\boldsymbol{r}}};\quad n_{{\boldsymbol{p}}}=\int {\rm d}{\boldsymbol{r}} \,n_{{\boldsymbol{p}},{\boldsymbol{r}}}; \quad \int {\rm d}{\boldsymbol{r}} \,\tilde{n}_{{\boldsymbol{p}},{\boldsymbol{r}}} =0 , \end{equation}

and acknowledge the presence of a mean drift in the spectral domain by adding zero

(3.22)\begin{equation} \dot{{\boldsymbol{p}}} = \dot{{\boldsymbol{p}}} - \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle + \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle , \end{equation}

in which $\langle \dots \rangle$ is an ensemble average for a system with spatially homogenous statistics, so that $\langle \dots \rangle$ is independent of ${\boldsymbol {r}}$. The mean drift arises due to inhomogeneities of the ray-path statistics in the spectral domain. Please note that the dimensions of $\bar {n}_{{\boldsymbol {p}}}$ and $n_{{\boldsymbol {p}}}$ are different.

Starting from the flux form of the action balance (3.17), we substitute (3.21a–c), (3.22) and invoke an ensemble average and integrate over ${\boldsymbol {r}}$ to obtain

(3.23)\begin{equation} \frac{\partial \langle n_{{\boldsymbol{p}}} \rangle }{\partial t}={-} \int {\rm d}{\boldsymbol{r}} \langle \boldsymbol{\nabla}_{{\boldsymbol{p}}} \,\boldsymbol{\cdot}\,[\kern0.7pt \dot{{\boldsymbol{p}}} - \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle ] n_{{\boldsymbol{p}},{\boldsymbol{r}}} \rangle - \boldsymbol{\nabla}_{\boldsymbol{p}} \,\boldsymbol{\cdot}\, [ \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle \langle n_{{\boldsymbol{p}}} \rangle ] . \end{equation}

Closure of this equation depends upon writing $n_{{\boldsymbol {p}},{\boldsymbol {r}}}$ in terms of $n_{{\boldsymbol {p}}}$. At this juncture, we invoke that property that wave-action spectral density does not change along trajectories

(3.24)\begin{align} n_{{\boldsymbol{p}},{\boldsymbol{r}}} &\equiv n(\kern1pt {\boldsymbol{p}}(t),{\boldsymbol{r}}(t),t) = n(\kern1pt {\boldsymbol{p}}(t-T),{\boldsymbol{r}}(t-T),t-T) \nonumber\\ & = \bar{n}\left[{\boldsymbol{p}}(t) - \int_{t-T}^t \dot{p}(t')\,{\rm d} t'; t-T\right] \nonumber\\ &\quad + \tilde{n}\left[{\boldsymbol{p}}(t) - \int_{t-T}^t \dot{p}(t')\,{\rm d} t'; {\boldsymbol{r}}(t) - \int_{t-\tau}^t \dot{{\boldsymbol{r}}} (t') \,{\rm d} t'; t-T\right]. \end{align}

We then execute a Taylor series expansion of $\bar {n}$

(3.25)\begin{align} n(\kern1pt {\boldsymbol{p}}(t-\tau),{\boldsymbol{r}}(t-\tau),t-\tau) \simeq \bar{n}(\kern1pt {\boldsymbol{p}}(t); t-T) - \boldsymbol{\nabla}_p \bar{n}(\kern1pt {\boldsymbol{p}}(t); t-T) \,\boldsymbol{\cdot}\, \int_{t-T}^t \dot{{\boldsymbol{p}}}(t') \,{\rm d} t' + \tilde{n}(\dots) , \end{align}

substitute, add zero once again and, using the definition of ensemble averaging $\langle \dots \rangle$,

(3.26) \begin{equation} \langle \boldsymbol{\nabla}_{{\boldsymbol{p}}} [\kern0.7pt \dot{{\boldsymbol{p}}} - \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle] \bar{n} (\kern1pt {\boldsymbol{p}}(t), t-\tau) \rangle \cong \boldsymbol{\nabla}_{{\boldsymbol{p}}} \left\langle \kern0.7pt\dot{{\boldsymbol{p}}} - \left\langle \kern0.7pt\dot{{\boldsymbol{p}}} \right\rangle \right\rangle \left\langle \bar{n} (\kern1pt {\boldsymbol{p}}(t), t-\tau) \right\rangle = 0, \end{equation}

and

(3.27)\begin{align} &\left\langle \boldsymbol{\nabla}_{{\boldsymbol{p}}} \left[\kern0.7pt \dot{{\boldsymbol{p}}} - \left\langle \kern0.7pt\dot{{\boldsymbol{p}}} \right\rangle\right] \int_{t-\tau}^t \left\langle \kern0.7pt\dot{{\boldsymbol{p}}} \right\rangle {\rm d}t' \,\boldsymbol{\cdot}\,\boldsymbol{\nabla}_{{\boldsymbol{p}}} \bar{n} (\kern1pt {\boldsymbol{p}}(t), t-\tau) \right\rangle \nonumber\\ &\quad \cong \boldsymbol{\nabla}_{{\boldsymbol{p}}} \left\langle \kern0.7pt\dot{{\boldsymbol{p}}} - \left\langle \kern0.7pt\dot{{\boldsymbol{p}}} \right\rangle \right\rangle \int_{t-\tau}^t\left\langle \kern0.7pt\dot{{\boldsymbol{p}}} \right\rangle {\rm d}t' \,\boldsymbol{\cdot}\,\boldsymbol{\nabla}_{{\boldsymbol{p}}} \left\langle \bar{n} (\kern1pt {\boldsymbol{p}}(t), t-\tau) \right\rangle\nonumber\\ &\quad = 0, \end{align}

as $\langle \dot {{\boldsymbol {p}}} - \langle \dot {{\boldsymbol {p}}} \rangle \rangle \equiv 0$ and neglecting an initial transient term

(3.28)\begin{equation} \langle \boldsymbol{\nabla}_{{\boldsymbol{p}}} [\kern0.7pt \dot{{\boldsymbol{p}}}(t) - \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle] \tilde{n}(\kern1pt {\boldsymbol{p}}(t-\tau),{\boldsymbol{r}}(t-\tau),t-\tau) \rangle, \end{equation}

so that (3.23) becomes

(3.29)\begin{equation} \frac{\partial \langle n_{{\boldsymbol{p}}} \rangle}{\partial t}={-} \boldsymbol{\nabla}_{\boldsymbol{p}} \,\boldsymbol{\cdot}\, \int_{t-\tau}^t \langle [\kern0.7pt \dot{{\boldsymbol{p}}}(t) - \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle] [\kern0.7pt \dot{{\boldsymbol{p}}}(t'-\tau) - \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle ]\rangle\,{\rm d}t' \boldsymbol{\nabla}_{{\boldsymbol{p}}} \langle n_{{\boldsymbol{p}}} \rangle - \boldsymbol{\nabla}_{{\boldsymbol{p}}} \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle \langle n_{{\boldsymbol{p}}} \rangle. \end{equation}

We introduce the auto-lag covariance matrix

(3.30) \begin{equation} {\mathcal{C}}_{ij}(\kern1pt {\boldsymbol{p}},t,t') = \langle [\kern0.7pt \dot{{\boldsymbol{p}}}({\boldsymbol{r}}(t)) - \langle \kern0.7pt\dot{{\boldsymbol{p}}}({\boldsymbol{r}}(t)) \rangle]_i [\dot{{\boldsymbol{p}}}({\boldsymbol{r}}(t')) -\langle \kern0.7pt\dot{{\boldsymbol{p}}}({\boldsymbol{r}}(t'))\rangle ]_j\rangle. \end{equation}

The final result is then

(3.31)\begin{equation} \frac{\partial \langle n_{{\boldsymbol{p}}} \rangle}{\partial t}={-} \boldsymbol{\nabla}_{{\boldsymbol{p}}_i}\,\boldsymbol{\cdot}\, \int^t_{t-\tau} {\mathcal{C}}_{ij}(\kern1pt {\boldsymbol{p}},t,t^{\prime}) \,{\rm d}t^{\prime} \,\boldsymbol{\cdot}\,\boldsymbol{\nabla}_{{\boldsymbol{p}}_j} \langle n_{\boldsymbol{p}} \rangle - \boldsymbol{\nabla}_{{\boldsymbol{p}}_i} \langle \kern0.7pt\dot{{\boldsymbol{p}}}_i \rangle \langle n_{\boldsymbol{p}} \rangle. \end{equation}

Convergence of the time integral, in which one can replace the lower limit of integration by $-\infty$, is the hallmark of a Markov approximation which we investigate further in Polzin & Lvov (Reference Polzin and Lvov2023). This equation encapsulates our fundamental theoretical result: the transport equation changes from diffusion (2.25) to an expression that involves both advection and diffusion in the spectral domain. Instead of being a no-flux stationary state, the GM76 spectrum, for which $\langle n_{\boldsymbol {p}} \rangle \propto m^0$, now supports a downscale action flux.

Our decomposition of wavenumber tendency into mean drift and dispersion about that mean drift provides a concrete mathematical interpretation for an interaction time scale $\tau _i$

(3.32)\begin{equation} \tau_i^{{-}1} = | \langle \kern0.7pt\dot{{\boldsymbol{p}}} \rangle | / | {\boldsymbol{p}} | \end{equation}

and correlation time scale $\tau _c$:

(3.33)\begin{equation} \tau_c = \int^t_{-\infty} {\mathcal{C}}_{ii}(\kern1pt {\boldsymbol{p}},t,t^{\prime}) \,{\rm d}t^{\prime} / \langle | \dot{{\boldsymbol{p}}}({\boldsymbol{r}}(t)) - \langle \kern0.7pt\dot{{\boldsymbol{p}}}({\boldsymbol{r}}(t)) |^2 \rangle \rangle . \end{equation}

Intuitive notions of the interplay between $\tau _i$ and $\tau _c$ are discussed in Müller et al. (Reference Müller, Holloway, Henyey and Pomphrey1986) and in Nazarenko et al. (Reference Nazarenko, Newell and Galtier2001) using expressions for ${\mathcal {C}}_{ij}$ (3.30) in which the ensemble mean drift has not been subtracted. The sentiment in Müller et al. (Reference Müller, Holloway, Henyey and Pomphrey1986) is that a separation between $\tau _c$ and $\tau _i$ is problematic for the oceanic internal wave field. Numerical ray-tracing results (Henyey & Pomphrey Reference Henyey and Pomphrey1983; Henyey, Pomphrey & Meiss Reference Henyey, Pomphrey and Meiss1984) do not elucidate why this might be. Our decomposition of wavenumber tendency into mean drift and dispersion about that mean drift, the revised Fokker–Planck (3.31) and revised covariance matrix (3.30) provide a concrete mathematical interpretation for such judgements about $\tau _i$ and $\tau _c$.

Having summoned the analogy between particle dispersion Taylor (Reference Taylor1921) and dispersion of wave packets in wavenumber, we are led to a degree of scepticism concerning the McComas & Bretherton (Reference McComas and Bretherton1977) and Nazarenko et al. (Reference Nazarenko, Newell and Galtier2001) interpretation that ray tracing should collapse onto the Fokker–Planck equation and diffusivity derived from the resonant kinetic equation (2.25). Convergence of the time-lagged auto-covariance (3.33) relates a finite diffusivity to the product of a covariance and correlation time scale. If one casts this as a resonant process, the covariance will be infinitely small and the correlation time scale infinitely long, thus leading to an inconsistency between interaction and correlation time scales in the resonant limit. Introducing a broadened kinetic equation (2.2.1) with finite bandwidth (2.19) might naively be considered progress. However, that bandwidth represents a spike decay rate rather than a correlation time scale and is not a resolution. As one runs to the finite amplitude of the weakly nonlinear problem, the resonant bandwidth becomes the r.m.s. Doppler shift (Polzin & Lvov Reference Polzin and Lvov2017). This is aphysical.

We find through numerical experimentation in Polzin & Lvov (Reference Polzin and Lvov2023) that the mean drift $\langle \dot {{\boldsymbol {p}}} \rangle$ is a resonant process and dispersion about the mean drift ${\mathcal {C}}_{ij}$ is non-resonant. The latter should not come as a surprise once one appreciates the direct analogy between particle dispersion and ray tracing originally suggested in McComas & Bretherton (Reference McComas and Bretherton1977): particle dispersion in turbulence has nothing to do with the concept of resonance. The moments $\langle \dot {{\boldsymbol {p}}} \rangle$ and ${\mathcal {C}}_{ij}$, and changes in the structure and scaling of the resonant bandwidth, are the signature differences of ray theory vs the kinetic equation, paralleling differences between Eulerian (Kraichnan Reference Kraichnan1959) and Lagrangian (Kraichnan Reference Kraichnan1965) representations of 3-D turbulence. These differences are rooted in the distinctions between amplitude-modulated and frequency-modulated signals.