1 Introduction
The propagation of inertia–gravity waves in the ocean has received a great deal of attention, mainly motivated by the role they play in the large- and mesoscale circulation, through wave–mean-flow interaction, mixing and dissipation. The inertia–gravity wave spectrum is dominated by two types of waves: near-inertial oscillations, with frequencies close to the inertial frequency
$f$
, which are mainly generated by winds, and internal tides (ITs), primarily at the semi-diurnal lunar frequency, which are generated by the interaction of the barotropic tide with topography (e.g. Ferrari & Wunsch Reference Ferrari and Wunsch2009). Near-inertial oscillations have distinctive dynamics, including weak dispersion and weak vertical motion, that stems from their unique place at the low-frequency end of the inertia–gravity wave spectrum (Alford et al.
Reference Alford, McKinnon, Simmons and Nash2016); ITs, in contrast, are generic mid-frequency inertia–gravity waves with their externally imposed frequency as their sole defining property.
The ocean’s highly energetic quasigeostrophic turbulence has a strong impact on the structure of both inertial oscillations and ITs and hence, in the case of ITs, on their signature on the sea-surface height (Rainville & Pinkel Reference Rainville and Pinkel2006; Ray & Zaron Reference Ray and Zaron2016). There is by now an extensive literature devoted to this impact, with a recent impetus provided by upcoming high-resolution satellite-altimetry instruments and the need to disentangle ITs from mesoscale (balanced) motion in the observed sea-surface height. We refer the reader to the recent papers by Dunphy et al. (Reference Dunphy, Ponte, Klein and Le Gentil2017) and Wagner, Ferrando & Young (Reference Wagner, Ferrando and Young2017) for further background.
A key aspect of the interactions between quasigeostrophic turbulence and both near-inertial oscillations and low-mode ITs is that turbulence and waves share similar horizontal scales, of the order of 100 km. A consequence is that, in such a regime, the Wentzel–Kramers–Brillouin (WKB) approximation on which much of our understanding of inertia–gravity wave propagation is built is not valid. This has prompted the development of simplified, wave-averaged models that rely only on time scale separation to represent the interactions between waves and flow in a simplified manner. Models of this kind include the Young–Ben Jelloul model of near-inertial oscillations in a quasigeostrophic flow (Young & Ben Jelloul Reference Young and Ben Jelloul1997) and its extensions accounting for the feedback of the waves on the flow (Xie & Vanneste Reference Xie and Vanneste2015; Wagner & Young Reference Wagner and Young2016; Thomas, Smith & Bühler Reference Thomas, Smith and Bühler2017). Wagner et al. (Reference Wagner, Ferrando and Young2017) recently derived an analogue of the Young–Ben Jelloul model equation for ITs. This model is formulated in physical space and retains a stiff term which enforces the constraint that the ITs’ fixed frequency imposes on their spatial structure and which cannot be eliminated without resorting to a Fourier-space formulation.
The present paper focuses similarly on the impact of a quasigeostrophic flow on ITs, making no asymptotic assumption about the relative horizontal scales of ITs and flow. Our focus is on quantifying the scattering induced by a barotropic (i.e.
$z$
-independent) geostrophic turbulent flow which we model as a spatially homogeneous random field. We take advantage of the assumption of barotropic flow to use a vertical-mode expansion and thus reduce the problem to the study of an (equivalent) shallow-water model. By applying the theory of waves in random media developed by Ryzhik, Papanicolaou & Keller (Reference Ryzhik, Papanicolaou and Keller1996), we derive a kinetic equation describing, in a statistically averaged sense, the energy exchanges between ITs with different wavevectors. The theory is formulated in terms of a wavevector-resolving energy density,
$a(\boldsymbol{x},\boldsymbol{k},t)$
, which makes it possible to capture spatial variations of the wave energy. The form of the scattering term in the kinetic equation for
$a(\boldsymbol{x},\boldsymbol{k},t)$
shows that energy transfers are restricted to waves with the same frequency or, equivalently, the same wavenumber
$|\boldsymbol{k}|$
. These transfers result from interactions within resonant triads consisting of two ITs of equal frequencies with a zero-frequency flow (vortical) mode – the so-called catalytic interactions of Lelong & Riley (Reference Lelong and Riley1991) and Bartello (Reference Bartello1995). The rate of these transfers is proportional to the energy spectrum of the geostrophic flow. In the case of an isotropic flow, the scattering leads to the relaxation of the energy density towards a locally isotropic density
$a(\boldsymbol{x},|\boldsymbol{k}|)$
. Our theory complements that developed by Ward & Dewar (Reference Ward and Dewar2010), shifting from a deterministic to a statistical treatment that can be regarded as a version of wave turbulence (e.g. Nazarenko Reference Nazarenko2011) in which the statistics of the flow are prescribed. It generalises the theory developed by Danioux & Vanneste (Reference Danioux and Vanneste2016) for near-inertial oscillations to inertia–gravity waves of arbitrary frequencies. Note that our statistical approach focuses on a homogeneous field of scatterers (resulting from the turbulent flow) and that different analysis techniques apply to waves incident on isolated scatterers (e.g. Olbers Reference Olbers1981).
We analyse the predictions of the kinetic equation, focusing our attention on parameters representative of the first baroclinic mode of the semidiurnal lunar tide
$M_{2}$
. These predictions include a time scale for wave isotropisation applicable to statistically homogeneous wave fields (i.e. such that
$\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a=0$
) and in particular to the isotropisation of an initially plane wave examined numerically by Ward & Dewar (Reference Ward and Dewar2010) in a shallow-water set-up. The kinetic equation applies to more general, non-homogeneous situations in which the energy density is modulated spatially (
$\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a\neq 0$
). This makes it possible to study the scattering of ITs generated by a localised source such as a topographic ridge. Ponte & Klein (Reference Ponte and Klein2015) and Dunphy et al. (Reference Dunphy, Ponte, Klein and Le Gentil2017) recently used three-dimensional Boussinesq simulations to study this problem and quantify the temporal incoherence of the ITs that results from the presence of a time-dependent turbulent flow. (See also Kelly et al. (Reference Kelly, Lermusiaux, Duda and Haley2016) and Kelly & Lermusiaux (Reference Kelly and Lermusiaux2016) for simulations of ITs in realistic configurations.) We carry out shallow-water simulations in a set-up analogous to theirs, and compare the results with direct solutions of the kinetic equation. This provides an estimate for the length scale over which the wave field becomes isotropic and, more broadly, sheds light on the interplay between transport of the wave energy by the group velocity and scattering. We emphasise that our theory concentrates on the statistical properties of the IT energy and makes no predictions for their phase. In the regime we consider, with a flow assumed to vary on a time scale much larger than the tidal period, the stationarity of the turbulent-energy spectrum ensures that the tidal energy remains concentrated at the single wavenumber dictated by the fixed tide frequency.
The plan of the paper is as follows. We describe the equations satisfied by linear internal waves propagating on a barotropic quasigeostrophic flow in § 2, expanding them in vertical modes to obtain an equivalent shallow-water system for each mode. We then sketch the derivation of the kinetic equation using the method of Ryzhik et al. (Reference Ryzhik, Papanicolaou and Keller1996), relegating the technical computations to appendix B. We focus on the application of the kinetic equation to the case of an isotropic flow in § 3 where we derive explicit estimates for the time- and length scales over which the IT field becomes isotropic. In § 4 we compare our theoretical predictions with direct simulations of the linearised shallow-water equations and with numerical solutions of the kinetic equation itself in a configuration where a wavemaker forces a plane IT in a turbulent flow. We conclude in § 5 with a discussion.
2 Scattering theory for internal tides
2.1 Model
We model the propagation of ITs through a turbulent quasigeostrophic eddy field using the hydrostatic Boussinesq equations linearised about a slowly evolving barotropic flow. The background flow is time dependent and geostrophically and hydrostatically balanced, given by
$\boldsymbol{U}=(U,V,0)=(-\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D713},\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D713},0)$
in terms of a
$z$
-independent streamfunction
$\unicode[STIX]{x1D713}$
. With these assumptions, the linearised hydrostatic Boussinesq equations read
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\boldsymbol{u}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U}+f\hat{\boldsymbol{z}}\times \boldsymbol{u}=-\unicode[STIX]{x1D735}p,\\ \unicode[STIX]{x2202}_{z}\,p=b,\\ \unicode[STIX]{x2202}_{t}b+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}b+N^{2}w=0,\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}+\unicode[STIX]{x2202}_{z}w=0,\end{array}\right\}\end{eqnarray}$$
where
$(\boldsymbol{u},w)$
denotes the IT velocity,
$\hat{\boldsymbol{z}}$
is the vertical unit vector,
$p$
is the pressure normalised by a reference density,
$b$
the buoyancy,
$f$
the Coriolis parameter and
$N(z)$
the buoyancy frequency. We use the notation
$\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}_{x},\unicode[STIX]{x2202}_{y},0)$
for the horizontal gradient throughout.
Assuming a flat bottom boundary and rigid lid, we project (2.1) onto baroclinic modes to obtain a set of rotating shallow-water equations governing their amplitudes:
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{t}\boldsymbol{u}_{m}+\boldsymbol{u}_{m}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{U}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}_{m}+f\hat{\boldsymbol{z}}\times \boldsymbol{u}_{m}=-g\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}_{m},\\ \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D702}_{m}+\boldsymbol{U}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D702}_{m}+h_{m}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}_{m}=0,\end{array}\right\}\end{eqnarray}$$
where
$\unicode[STIX]{x1D702}_{m}$
is the equivalent surface height and
$h_{m}$
is the equivalent depth (see appendix A for details). Note that this system differs from the one obtained by linearising the shallow-water equations about a background flow in geostrophic balance since the latter system includes a contribution from the (sloping) background free surface.
For physical applications in later sections, we take parameters corresponding to the first baroclinic mode only, since this contains the majority of the IT energy. We drop the subscript
$m$
in (2.2) from this point on. In the ocean, energy is transferred between vertical modes as a result of vertical shear. However, as discussed by Dunphy & Lamb (Reference Dunphy and Lamb2014) and Ponte & Klein (Reference Ponte and Klein2015) the effect is small. Simulations in Dunphy et al. (Reference Dunphy, Ponte, Klein and Le Gentil2017) put the transfer of energy from the first mode to higher modes at 3 % in their most extreme cases, with a highly energetic background flow, and less for typical ocean conditions.
2.2 Derivation of the kinetic equation
We study IT scattering in the distinguished regime where the spatial scale of the flow,
$L_{\ast }$
say, is of the same order as the wavelength, that is,
$|\boldsymbol{k}|L_{\ast }=O(1)$
, where
$\boldsymbol{k}=(k,l)$
is the IT horizontal wavevector. The assumption of a geostrophic flow requires a small Rossby number
$Ro=U_{\ast }/(fL_{\ast })\ll 1$
, where
$U_{\ast }$
is a typical flow velocity; in turn, this implies that the background flow velocities are small compared with the wave phase speed
$\unicode[STIX]{x1D714}/|\boldsymbol{k}|$
, where
$$\begin{eqnarray}\unicode[STIX]{x1D714}=\sqrt{f^{2}+gh|\boldsymbol{k}|^{2}}\end{eqnarray}$$
is the IT frequency, since
$U_{\ast }/(\unicode[STIX]{x1D714}/|\boldsymbol{k}|)=O(U_{\ast }/(\unicode[STIX]{x1D714}L_{\ast }))=O(Ro)$
, given that
$\unicode[STIX]{x1D714}=O(f)$
away from the equator. With the flow time scale
$T_{\ast }$
taken as the natural advective time scale
$L_{\ast }/U_{\ast }$
, this also implies that the flow evolves slowly compared with the IT time scale since
$\unicode[STIX]{x1D714}T_{\ast }=O(Ro)\ll 1$
. We further assume that, while the IT phases vary over the length scale
$|\boldsymbol{k}|^{-1}$
, their amplitudes vary over a much larger scale
$\unicode[STIX]{x1D700}^{-1}|\boldsymbol{k}|$
, where
$\unicode[STIX]{x1D700}\ll 1$
. We adopt the scaling
$\unicode[STIX]{x1D700}=O(Ro^{2})$
. As emerges below, this is the distinguished scaling that ensures that transport and scattering affect the wave field at the same order and are both captured at leading order by our asymptotic model.
Since our focus is on generic, statistical properties of the IT field, we model the turbulent background flow by a random streamfunction with homogeneous and stationary statistics. With our scaling assumptions, it is then possible to derive a single equation that describes the scattering and transport of IT energy following the theory of Ryzhik et al. (Reference Ryzhik, Papanicolaou and Keller1996). This theory is formulated in terms of the Wigner transform
$$\begin{eqnarray}\unicode[STIX]{x1D652}(\boldsymbol{x},\boldsymbol{k},t)=\frac{1}{(2\unicode[STIX]{x03C0})^{2}}\int _{\mathbb{R}^{2}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{y}}\unicode[STIX]{x1D753}(\boldsymbol{x}-\boldsymbol{y}/2,t)\unicode[STIX]{x1D753}^{\ast }(\boldsymbol{x}+\boldsymbol{y}/2,t)\,\text{d}\boldsymbol{y},\end{eqnarray}$$
where
$\unicode[STIX]{x1D753}=(u,v,\unicode[STIX]{x1D702})^{\text{T}}$
groups the dynamical fields and the asterisk denotes conjugate transpose. The Wigner transform is a Hermitian matrix; it is not necessarily positive definite, although its integral over the wavevector
$\boldsymbol{k}$
, simply given by
$\unicode[STIX]{x1D753}(\boldsymbol{x},t)\unicode[STIX]{x1D753}^{\ast }(\boldsymbol{x},t)$
, is.
The equation derived by Ryzhik et al. (Reference Ryzhik, Papanicolaou and Keller1996) governs the evolution of a scalar amplitude
$a(\boldsymbol{x},\boldsymbol{k},t)$
which appears naturally in an eigenvector decomposition of the matrix
$\unicode[STIX]{x1D652}(\boldsymbol{x},\boldsymbol{k},t)$
(see appendix B for details). Physically, this amplitude is interpreted as a wavevector-resolving energy density, related to the (leading-order) energy density of the system by
$$\begin{eqnarray}{\mathcal{E}}_{0}(\boldsymbol{x},t)=\frac{1}{2}\int _{\mathbb{R}^{2}}a(\boldsymbol{x},\boldsymbol{k},t)\,\text{d}\boldsymbol{k}.\end{eqnarray}$$
To avoid any confusion, we emphasise that
$a(\boldsymbol{x},\boldsymbol{k},t)$
itself represents a wave-energy density and not a wave amplitude; as its definition in (B 29) makes clear, it is a quadratic function of the wave fields, like the Wigner transform
$\unicode[STIX]{x1D652}(\boldsymbol{x},\boldsymbol{k},t)$
.
In appendix B, we show that
$a(\boldsymbol{x},\boldsymbol{k},t)$
satisfies the kinetic equation
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}a+\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a={\mathcal{L}}a-\unicode[STIX]{x1D6F4}a,\end{eqnarray}$$
where the notation
$\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}=\unicode[STIX]{x1D735}=(\unicode[STIX]{x2202}_{x},\unicode[STIX]{x2202}_{y})$
emphasises that the spatial gradient applies to functions of both
$\boldsymbol{x}$
and
$\boldsymbol{k}$
, and where
$\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}=(\unicode[STIX]{x2202}_{k},\unicode[STIX]{x2202}_{l})$
. Here
$\unicode[STIX]{x1D714}$
is determined by the IT dispersion relation (2.3), so that
$\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}$
is the group velocity and the left-hand side of (2.6) represents the familiar wave transport. (The term
$-\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}a$
would be added if
$\unicode[STIX]{x1D714}$
depended explicitly on
$\boldsymbol{x}$
.) The right-hand side represents wave scattering by the background flow. The first term, given by
$$\begin{eqnarray}{\mathcal{L}}a(\boldsymbol{x},\boldsymbol{k},t)=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })a(\boldsymbol{x},\boldsymbol{k}^{\prime },t)\,\text{d}\boldsymbol{k}^{\prime },\end{eqnarray}$$
quantifies the transfers of energy from all wavevectors
$\boldsymbol{k}^{\prime }$
into wavevector
$\boldsymbol{k}$
that result from interactions with the background flow; the second term, where
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}=\unicode[STIX]{x1D6F4}(\boldsymbol{k})=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })\,\text{d}\boldsymbol{k}^{\prime },\end{eqnarray}$$
is the total scattering cross-section, quantifies the energy lost by wavevector
$\boldsymbol{k}$
to all other wavevectors.
The function
$\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })$
that appears in (2.7)–(2.8) is the main object of interest. It is known as the differential scattering cross-section and measures the rate at which energy is scattered from
$\boldsymbol{k}^{\prime }$
to
$\boldsymbol{k}$
at a position
$\boldsymbol{x}$
in space. We obtain it in the form
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime }) & = & \displaystyle \frac{2\unicode[STIX]{x03C0}}{gh\unicode[STIX]{x1D714}^{3}|\boldsymbol{k}|^{5}}\{\!|\boldsymbol{k}\times \boldsymbol{k}^{\prime }|^{2}[(\unicode[STIX]{x1D714}^{2}+f^{2})\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}^{\prime }-f^{2}|\boldsymbol{k}|^{2}]^{2}\nonumber\\ \displaystyle & & \displaystyle +\,f^{2}\unicode[STIX]{x1D714}^{2}[|\boldsymbol{k}\times \boldsymbol{k}^{\prime }|^{2}+\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}^{\prime }(|\boldsymbol{k}|^{2}-\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}^{\prime })]^{2}\!\}\frac{\widehat{E}(\boldsymbol{k}-\boldsymbol{k}^{\prime })}{|\boldsymbol{k}-\boldsymbol{k}^{\prime }|^{2}}\unicode[STIX]{x1D6FF}(|\boldsymbol{k}|-|\boldsymbol{k}^{\prime }|),\end{eqnarray}$$
where
$\boldsymbol{k}\times \boldsymbol{k}^{\prime }=kl^{\prime }-lk^{\prime }$
, and
$\widehat{E}$
is the energy spectrum of the flow. We note that
$\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })$
is real, positive, and symmetric with respect to the exchange between
$\boldsymbol{k}$
and
$\boldsymbol{k}^{\prime }$
. In appendix B, we also show that these properties ensure conservation of the leading-order energy density (2.5):
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}{\mathcal{E}}_{0}+\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}\boldsymbol{\cdot }{\mathcal{F}}_{0}=0,\end{eqnarray}$$
where
$$\begin{eqnarray}{\mathcal{F}}_{0}(\boldsymbol{x},t)=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}(\boldsymbol{k})\,a(\boldsymbol{x},\boldsymbol{k},t)\,\text{d}\boldsymbol{k}\end{eqnarray}$$
is the leading-order energy flux (see (B 56)).
The presence of the factor
$\unicode[STIX]{x1D6FF}(|\boldsymbol{k}|-|\boldsymbol{k}^{\prime }|)$
in (2.9) implies that energy is only exchanged between wavevectors of the same magnitude, that is, between waves with the same frequency, as a result of the assumed slow time dependence of the background flow. Thus, in the regime considered, the IT energy is confined to the constant-frequency circle
$|\boldsymbol{k}|=((\unicode[STIX]{x1D714}^{2}-f^{2})/(gh))^{1/2}$
in the wavevector plane. This can be related to the observation that the background flow only enters
$\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })$
through its energy spectrum
$\widehat{E}$
, which, for the statistically stationary flows considered, is time independent. The scattering described by (2.6) results from the resonant interactions of two ITs, with wavevectors
$\boldsymbol{k}$
and
$\boldsymbol{k}^{\prime }$
and identical frequencies, with a vortical flow mode of wavevector
$\boldsymbol{k}-\boldsymbol{k}^{\prime }$
and zero frequency. Because of potential-vorticity conservation, these interactions would leave the vortical mode unaffected even if it were allowed to evolve freely; hence they have been termed catalytic interactions (Lelong & Riley Reference Lelong and Riley1991; Bartello Reference Bartello1995; Ward & Dewar Reference Ward and Dewar2010). We emphasise that (2.6) captures the net effect of multiple triadic interactions acting over long time scales. This is why the time scale of evolution is not linear in the flow amplitude but quadratic, dictated by the energy spectrum of the flow, in a manner familiar from wave turbulence (e.g. Nazarenko Reference Nazarenko2011).
3 Scattering in isotropic turbulence
3.1 Isotropisation
In this section we use the kinetic equation (2.6) to make predictions about the scattering process and quantify the time and length scales over which ITs lose their spatial coherence. For simplicity, we assume that the flow is isotropic,
$\widehat{E}(\boldsymbol{k})=\widehat{E}(|\boldsymbol{k}|)$
. This motivates the use of polar coordinates for the wavevector, such that
$$\begin{eqnarray}\boldsymbol{k}=|\boldsymbol{k}|\left(\begin{array}{@{}c@{}}\cos \unicode[STIX]{x1D703}\\ \sin \unicode[STIX]{x1D703}\end{array}\right)\quad \text{and}\quad \boldsymbol{k}^{\prime }=|\boldsymbol{k}^{\prime }|\left(\begin{array}{@{}c@{}}\cos (\unicode[STIX]{x1D703}+\unicode[STIX]{x1D703}^{\prime })\\ \sin (\unicode[STIX]{x1D703}+\unicode[STIX]{x1D703}^{\prime })\end{array}\right),\end{eqnarray}$$
where
$\unicode[STIX]{x1D703}^{\prime }$
is the angle between
$\boldsymbol{k}$
and
$\boldsymbol{k}^{\prime }$
. The change of coordinates reduces the scattering operator (2.7) to
$$\begin{eqnarray}{\mathcal{L}}a(\boldsymbol{x},|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)=\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703}^{\prime })a(\boldsymbol{x},|\boldsymbol{k}|,\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\prime },t)\,\text{d}\unicode[STIX]{x1D703}^{\prime },\end{eqnarray}$$
where
$$\begin{eqnarray}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703}^{\prime }):=\int _{0}^{\infty }\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })|\boldsymbol{k}^{\prime }|\,\text{d}|\boldsymbol{k}^{\prime }|.\end{eqnarray}$$
Note that we have used the evenness of
$\unicode[STIX]{x1D70E}^{\prime }$
in
$\unicode[STIX]{x1D703}^{\prime }$
to rewrite (3.2) as a convolution. Note also that
$\unicode[STIX]{x1D70E}^{\prime }$
is independent of the direction
$\unicode[STIX]{x1D703}$
of
$\boldsymbol{k}$
because the scattering process is rotationally invariant. The scattering cross-section (2.9) can be written explicitly as
$$\begin{eqnarray}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703})=\frac{\unicode[STIX]{x03C0}|\boldsymbol{k}|^{2}}{gh\unicode[STIX]{x1D714}^{3}}\{(\!1+\cos \unicode[STIX]{x1D703})[(\unicode[STIX]{x1D714}^{2}+f^{2})\cos \unicode[STIX]{x1D703}-f^{2}]^{2}+\unicode[STIX]{x1D714}^{2}f^{2}[\!1-\cos (3\unicode[STIX]{x1D703})]\}\widehat{E}(2|\boldsymbol{k}\sin (\unicode[STIX]{x1D703}/2)|),\end{eqnarray}$$
where we have removed the prime from
$\unicode[STIX]{x1D703}^{\prime }$
for convenience. Similarly, the total scattering cross-section (2.8) reduces to
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}(|\boldsymbol{k}|)=\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$
In the limit
$\unicode[STIX]{x1D714}\rightarrow f$
corresponding to near-inertial waves the cross-section reduces to
$$\begin{eqnarray}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703})=\frac{2\unicode[STIX]{x03C0}f|\boldsymbol{k}|^{2}}{gh}\widehat{E}(2|\boldsymbol{k}\sin (\unicode[STIX]{x1D703}/2)|),\end{eqnarray}$$
which recovers the result obtained by Danioux & Vanneste (Reference Danioux and Vanneste2016) starting from the Young–Ben Jelloul model.
With the scattering cross-section (3.4), the scattering operator (3.2) can be diagonalised using a Fourier series, or more precisely a cosine series since
$a$
is even in
$\unicode[STIX]{x1D703}$
. Denoting the cosine transform by a hat, with
$$\begin{eqnarray}\widehat{a}_{n}(\boldsymbol{x},|\boldsymbol{k}|,t)=\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}\cos (n\unicode[STIX]{x1D703})a(\boldsymbol{x},|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)\,\text{d}\unicode[STIX]{x1D703},\end{eqnarray}$$
we find that
$$\begin{eqnarray}(\widehat{{\mathcal{L}}a})_{n}=\unicode[STIX]{x1D706}_{n}\widehat{a}_{n},\quad n=0,1,\ldots ,\end{eqnarray}$$
with the eigenvalues
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{n}=\unicode[STIX]{x1D706}_{n}(|\boldsymbol{k}|):=2\unicode[STIX]{x03C0}\widehat{\unicode[STIX]{x1D70E}^{\prime }}=\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E}^{\prime }(|\boldsymbol{k}|,\unicode[STIX]{x1D703})\cos (n\unicode[STIX]{x1D703})\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$
Fourier transforming the kinetic equation (2.6) then gives
It follows from (3.9) and the non-negativity of
$\unicode[STIX]{x1D70E}^{\prime }$
in (3.4) that
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{0}=\unicode[STIX]{x1D6F4}(|\boldsymbol{k}|)\quad \text{and}\quad |\unicode[STIX]{x1D706}_{n\geqslant 1}|<\unicode[STIX]{x1D706}_{0}.\end{eqnarray}$$
Thus the scattering term on the right-hand side of (3.10) vanishes for
$n=0$
and represents a damping for
$n\geqslant 1$
.
The implications are clearly seen for a wave field that is spatially homogeneous, that is, with
$\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a=0$
: the solution of (2.6), with initial condition
$$\begin{eqnarray}a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t=0)=A(|\boldsymbol{k}|,\unicode[STIX]{x1D703}),\end{eqnarray}$$
is then simply
$$\begin{eqnarray}a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)=\mathop{\sum }_{n=0}^{\infty }\widehat{A}_{n}(|\boldsymbol{k}|)\text{e}^{(\unicode[STIX]{x1D706}_{n}-\unicode[STIX]{x1D6F4})t}\cos (n\unicode[STIX]{x1D703}).\end{eqnarray}$$
This describes the relaxation of the solution towards a stationary, isotropic (
$\unicode[STIX]{x1D703}$
-independent) solution, since
$$\begin{eqnarray}\lim _{t\rightarrow \infty }a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)=\widehat{A}_{0}(|\boldsymbol{k}|)=\frac{1}{2\unicode[STIX]{x03C0}}\int _{-\unicode[STIX]{x03C0}}^{\unicode[STIX]{x03C0}}a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t=0)\,\text{d}\unicode[STIX]{x1D703}.\end{eqnarray}$$
This is a key feature of the scattering: the main impact of the random isotropic flow is to lead to the isotropisation of the IT field regardless of the initial condition. Note that, with
$a(|\boldsymbol{k}|,\unicode[STIX]{x1D703},t)$
the wave-energy density,
$\widehat{A}_{0}$
represents the total,
$\unicode[STIX]{x1D703}$
-integrated energy at wavevector
$|\boldsymbol{k}|$
, while the
$\widehat{A}_{n}$
for
$n\neq 0$
capture the energy’s dependence on
$\unicode[STIX]{x1D703}$
.
We can identify two time scales for the scattering process. First, the scattering time
$$\begin{eqnarray}T_{\mathit{scat}}=\unicode[STIX]{x1D6F4}^{-1}\end{eqnarray}$$
estimates the time over which energy concentrated at
$\boldsymbol{k}$
in spectral space is reduced by a factor of
$\text{e}^{-1}$
while converted to waves with other wavevectors. In other words, it is the time scale over which scattering effects become significant. Second, the time scale for convergence to an isotropic wave field is given by
$$\begin{eqnarray}T_{\mathit{iso}}=(\unicode[STIX]{x1D6F4}-\unicode[STIX]{x1D706}^{\prime })^{-1},\quad \text{where}\;\unicode[STIX]{x1D706}^{\prime }:=\max _{n\geqslant 1}\unicode[STIX]{x1D706}_{n}.\end{eqnarray}$$
This is the time for the last surviving anisotropic (i.e.
$n\neq 0$
) Fourier mode to decay by a factor of
$\text{e}^{-1}$
. Scattering length scales associated with the time scales (3.15) and (3.16) can be defined as
$$\begin{eqnarray}L_{\mathit{scat}}=c_{g}T_{\mathit{scat}},\quad L_{\mathit{iso}}=c_{g}T_{\mathit{iso}},\end{eqnarray}$$
where
$c_{g}=|\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}|=gh|\boldsymbol{k}|/\unicode[STIX]{x1D714}(|\boldsymbol{k}|)$
is the group speed.
3.2 Predicted behaviour
In this section, we use the time and length scales (3.15)–(3.17) to examine how the scattering depends on the Coriolis parameter
$f$
, and on the strength and horizontal scales of the eddies as encoded in
$\widehat{E}$
. Since we focus on ITs, we regard the frequency
$\unicode[STIX]{x1D714}$
as fixed and deduce
$|\boldsymbol{k}|$
from the dispersion relation (2.3). We test some of our predictions against numerical simulations in § 4.
We assume an isotropic energy spectrum of the form
$$\begin{eqnarray}\widehat{{\mathcal{E}}}(|\boldsymbol{k}|):=2\unicode[STIX]{x03C0}|\boldsymbol{k}|\,\widehat{E}(|\boldsymbol{k}|)=\left\{\begin{array}{@{}ll@{}}c_{1}|\boldsymbol{k}|,\quad & |\boldsymbol{k}|\leqslant \unicode[STIX]{x1D705},\\ c_{2}|\boldsymbol{k}|^{-3.5},\quad & |\boldsymbol{k}|\geqslant \unicode[STIX]{x1D705}.\end{array}\right.\end{eqnarray}$$
This depends on two parameters:
$\unicode[STIX]{x1D705}$
, a peak wavenumber which sets the dominant length scale of the flow, and the root-mean-square velocity defined by
$v_{\mathit{rms}}^{2}=\int _{0}^{\infty }\widehat{{\mathcal{E}}}\,\text{d}|\boldsymbol{k}|$
. The constants
$c_{1}$
and
$c_{2}$
are determined by
$\unicode[STIX]{x1D705}$
and
$v_{\mathit{rms}}$
and the requirement of continuity at
$|\boldsymbol{k}|=\unicode[STIX]{x1D705}$
. In practice, we choose
$\unicode[STIX]{x1D705}$
so that the correlation length
$l_{c}:=\unicode[STIX]{x03C0}/k_{c}$
, where
$k_{c}:=\iint |\boldsymbol{k}|\widehat{E}\,\text{d}\boldsymbol{k}/\iint \widehat{E}\,\text{d}\boldsymbol{k}$
, is similar to the wavelength of the IT; calculation of the integrals using (3.18) gives
$\unicode[STIX]{x1D705}=9\unicode[STIX]{x03C0}/(10l_{c})$
. Although quasigeostrophic theory predicts a kinetic-energy spectrum that decays as
$|\boldsymbol{k}|^{-3}$
for balanced geostrophic turbulence, the slope is often observed to be slightly steeper, a result which is typically attributed to the presence of large-scale coherent structures that emerge in the turbulent flow (McWilliams, Weiss & Yavneh Reference McWilliams, Weiss and Yavneh1994; Bartello Reference Bartello1995; Kafiabad & Bartello Reference Kafiabad and Bartello2016). This motivates the form of (3.18) as representative of balanced geostrophic eddy fields in the ocean. Note that we have chosen an energy spectrum with non-zero energy for all
$|\boldsymbol{k}|$
. This matters because only the range
$[0,2|\boldsymbol{k}|]$
of the energy spectrum contributes to the scattering of ITs with wavenumber
$|\boldsymbol{k}|$
, as the factor
$\widehat{E}(2|\boldsymbol{k}\sin (\unicode[STIX]{x1D703}^{\prime }/2)|)$
in the scattering cross-section (3.4) indicates. A lower cutoff of the spectrum, say at some wavenumber
$k_{cut}$
, would then imply that waves with
$|\boldsymbol{k}|\leqslant k_{cut}/2$
are unaffected by scattering. This is related to the resonant triad view that two ITs with the same wavenumber
$|\boldsymbol{k}|$
and hence the same frequency can only form a resonant triad with a flow mode if this has a wavenumber in
$[0,2|\boldsymbol{k}|]$
.
Figure 1. Scattering cross-section
$\unicode[STIX]{x1D70E}^{\prime }$
in (3.4) for the energy spectrum (3.18) as a function of the peak wavenumber
$\unicode[STIX]{x1D705}$
and angle
$\unicode[STIX]{x1D703}$
treated as polar coordinates. The IT wavenumber is fixed as
$|\boldsymbol{k}|=3\times 10^{-5}~\text{m}^{-1}$
corresponding to the mode-1
$M_{2}$
tide at
$45^{\circ }$
latitude, and the flow’s root-mean-square velocity as
$v_{\mathit{rms}}=0.25~\text{m}~\text{s}^{-1}$
. (a) Peak wavenumber in the range
$0.2|\boldsymbol{k}|\leqslant \unicode[STIX]{x1D705}\leqslant |\boldsymbol{k}|$
, (b) in the range
$0.5|\boldsymbol{k}|\leqslant \unicode[STIX]{x1D705}\leqslant 2.5|\boldsymbol{k}|$
. Three (circular) contours of
$\unicode[STIX]{x1D705}$
are labelled in each panel in units of
$10^{-5}~\text{m}^{-1}$
.
Figure 1 shows the scattering cross-section
$\unicode[STIX]{x1D70E}^{\prime }$
in (3.4) as a function of
$\unicode[STIX]{x1D703}$
and
$\unicode[STIX]{x1D705}$
for
$v_{\mathit{rms}}=0.25~\text{m}~\text{s}^{-1}$
and for
$f=1.028\times 10^{-4}~\text{s}^{-1}$
corresponding to
$45^{\circ }$
latitude. The equivalent depth is set to
$h=1.2~\text{m}$
, as appropriate for the first baroclinic mode (Olbers, Willebrand & Eden Reference Olbers, Willebrand and Eden2012), and the frequency to
$\unicode[STIX]{x1D714}=2\unicode[STIX]{x03C0}/12.42~\text{h}^{-1}$
, corresponding to the
$M_{2}$
tide. The horizontal wavenumber is then
$|\boldsymbol{k}|=3\times 10^{-5}~\text{m}^{-1}$
corresponding to a wavelength of approximately
$200~\text{km}$
. The figure indicates that scattering is local in the angular coordinate, that is, ITs are preferentially scattered into waves with nearby directions. This is especially the case for small values of
$\unicode[STIX]{x1D705}$
, corresponding to flows with typical scales much larger than the IT wavelength (figure 1
a), when the values of
$\unicode[STIX]{x1D70E}^{\prime }$
are also the largest. For larger values of
$\unicode[STIX]{x1D705}$
, that is, for flows with smaller scales, the energy transfers are slower, but less localised in the angular direction (figure 1
b).
The net effect of the scattering depends on both the value of
$\unicode[STIX]{x1D70E}^{\prime }$
at fixed
$\unicode[STIX]{x1D703}$
and the range of
$\unicode[STIX]{x1D703}$
where
$\unicode[STIX]{x1D70E}^{\prime }$
is substantial; it is best measured by the scattering and isotropisation time and length scales introduced in (3.15)–(3.17). These scales are deduced from the cosine transform of
$\unicode[STIX]{x1D70E}^{\prime }$
which give the eigenvalues of the scattering operator. The eigenvalues are shown in figure 2(a) for
$\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$
, corresponding to a flow correlation length
$l_{c}\approx 180~\text{km}$
, with all other parameters as in figure 1. The most important eigenvalues are the two largest,
$n=0$
and here
$n=1$
, since they control the scattering and isotropisation time and length scales. These scales are displayed as functions of
$\unicode[STIX]{x1D705}$
in figure 2(b). The figure shows that large-scale flows lead to rapid scattering but slow isotropisation. This can be easily understood: large-scale flows cause rapid energy transfers but, because of the localised nature of the scattering, these transfers are limited to waves of similar directions and a long time is needed for energy to be distributed near uniformly in the angular direction. (The large-scale flow regime can be tackled using ray tracing as has frequently been applied for deterministic flows (e.g. Rainville & Pinkel Reference Rainville and Pinkel2006; Chavanne et al.
Reference Chavanne, Flament, Carter, Merrifield, Luther, Zaron and Gurgel2010). For weak random flows as assumed here, the ray equations can be analysed asymptotically using methods developed for noisy Hamiltonian systems (e.g. Bal, Komorowski & Ryzhik Reference Bal, Komorowski and Ryzhik2010) to show that the IT wavevector diffuses along the constant-frequency circle
$|\boldsymbol{k}|=\text{const.}$
, consistent with the kinetic-equation description; see Müller (Reference Müller1976, Reference Müller1977) for early treatments in this spirit.) Isotropisation is most effective when
$\unicode[STIX]{x1D705}$
has an order of magnitude similar to
$|\boldsymbol{k}|$
: for the chosen energy spectrum, isotropisation is fastest for
$\unicode[STIX]{x1D705}\approx 6\times 10^{-5}~\text{m}^{-1}$
corresponding to a flow correlation length of approximately 50 km. Isotropisation slows down for larger values of
$\unicode[STIX]{x1D705}$
simply because the total flow energy in the useful range
$[0,2|\boldsymbol{k}|]$
decreases with
$\unicode[STIX]{x1D705}$
.
Figure 2. (a) Eigenvalues
$\unicode[STIX]{x1D706}_{n}$
of the scattering operator, given by (3.9), for the energy spectrum in (3.18) with
$\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$
and
$v_{\mathit{rms}}=0.25~\text{m}~\text{s}^{-1}$
, and for an IT with
$|\boldsymbol{k}|=3\times 10^{-5}~\text{m}^{-1}$
and
$f=1.028\times 10^{-4}~\text{s}^{-1}$
. (b) Scattering and isotropisation length and time scales
$L_{scat}$
,
$L_{iso}$
,
$T_{scat}$
and
$T_{iso}$
as functions of the peak wavenumber
$\unicode[STIX]{x1D705}$
, with all the other parameters as in (a). (c) As in (b) but as functions of
$v_{\mathit{rms}}$
and for
$\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$
. (d) As in (b) but as functions of latitude.
Figure 2(c) shows the scattering and isotropisation times and lengths as functions of
$v_{\mathit{rms}}$
and for
$\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$
. The dependence is simply in
$v_{\mathit{rms}}^{-2}$
. The figure suggests that full isotropisation of ITs generated at localised topographical features is rare in the ocean since the length scales required exceed the basin scales even for strong flows. On the other hand, scattering is effective over much shorter spatial scales, of the order of a few hundreds of kilometres, and over time scales of a week or so, comparable to other dynamical time scales in the ocean. The conclusion, then, is that typical ITs are strongly influenced by the quasigeostrophic flow, though not to the extent that they become completely isotropic. As highlighted by Ward & Dewar (Reference Ward and Dewar2010), the time scale of a week or so is shorter than the characteristic time scales of nonlinear wave–wave interactions except, perhaps, for the special case of parametric subharmonic instability at the critical latitude of
$29^{\circ }$
(MacKinnon & Winters Reference MacKinnon and Winters2005). It is likely, then, that scattering by the geostrophic flow plays a more important role than wave–wave interactions in determining the characteristics of oceanic inertia–gravity waves. We should note, however, that the most energetic regions of the ocean, such as western boundary currents, where scattering is most effective, are also strongly inhomogeneous so that out theory does not apply in a strict sense.
Figure 2(d) explores the dependence of scattering on latitude. Latitude affects the cross-section
$\unicode[STIX]{x1D70E}^{\prime }$
through
$f$
and also through
$|\boldsymbol{k}|$
if we consider a fixed frequency as is done here. In the ocean, different latitudes may also lead to different energy spectra. For simplicity, in plotting figure 2(d) we have taken the same spectrum for each latitude, keeping
$\unicode[STIX]{x1D705}$
fixed. The figure shows that the scattering time increases with latitude. The isotropisation time, however, decreases with latitude with, as far as we can tell, no obvious interpretation; the scattering is determined as the difference between two eigenvalues and is hence difficult to intuit. The scattering and isotropisation lengths both decrease with latitude, partly as a result of a decrease of the group velocity. We note that Ward & Dewar (Reference Ward and Dewar2010) conclude from simulations that scattering and isotropisation weaken with latitude, leading to longer propagation distances (see their figure 11). This apparent contradiction is likely resolved by the fact that their non-dimensional formulation implies that their energy spectrum also changes with latitude, keeping the energy in the range
$[0,2|\boldsymbol{k}|]$
constant as
$f$
changes. A general conclusion we can draw from the form of
$\unicode[STIX]{x1D70E}^{\prime }$
and our parameter dependence study is the fact that the quasigeostrophic energy spectrum is the key factor determining the strength of the scattering.
4 Simulations
In this section we analyse numerical simulations of the linearised equivalent shallow-water system (2.2) and compare them with the theoretical predictions of the previous section and with direct simulations of the kinetic equation (2.6).
4.1 Shallow-water simulations
We solve (2.2) numerically, adding a harmonic forcing term to generate a coherent plane wave. The numerical scheme relies on pseudospectral and splitting methods: the terms independent of the background flow are integrated exactly in Fourier space, while the terms that depend on the background flow are integrated using an Euler scheme in physical space. The domain is a
$7168~\text{km}\times 1024~\text{km}$
channel on an
$f$
-plane centred at
$45\,^{\circ }\text{N}$
. We use a spatial resolution of
$1792\times 256$
, with that
$\unicode[STIX]{x0394}x=\unicode[STIX]{x0394}y=4~\text{km}$
, with periodic boundary conditions in the
$y$
-direction, and absorbing layers 30 grid points wide at each end of the domain in the
$x$
-direction, and take time steps of
$\unicode[STIX]{x0394}t=4000~\text{s}$
. We run an ensemble of 100 simulations with random realisations of the background flow in order to study statistics. Each simulation corresponds to 80 days, which is long enough to study isotropisation in a moderately energetic flow. A wavemaker forces an IT through a term of the form
$$\begin{eqnarray}{\mathcal{F}}=A\sin (\unicode[STIX]{x1D6FA}t)\text{e}^{-(x-x_{0})^{2}/\unicode[STIX]{x1D6E5}^{2}},\end{eqnarray}$$
added to the continuity equation (cf. Ponte & Klein Reference Ponte and Klein2015). Here,
$x_{0}=400~\text{km}$
is the position of the wavemaker in the
$x$
-direction and
$\unicode[STIX]{x1D6E5}=10~\text{km}$
its width;
$\unicode[STIX]{x1D6FA}$
is the tidal frequency. The amplitude
$A$
is arbitrary since we solve a linear system. The forcing is ramped up slowly to reach its maximum amplitude over approximately 1 week. The resulting plane waves that are generated have a wavelength of approximately 150 km, as expected for a first baroclinic mode wave at
$45^{\circ }$
latitude.
Figure 3. (a) Sea-surface elevation
$\unicode[STIX]{x1D702}$
in an equivalent shallow-water simulation of the mode-1
$M_{2}$
tide in a turbulent flow with
$v_{rms}=0.25~\text{m}~\text{s}^{-1}$
at
$45^{\circ }$
latitude. (b) Vorticity field of the turbulent flow.
For the background flow we take a homogeneous isotropic Gaussian random field, tapered in the region
$0<x<1000~\text{km}$
so as not to interfere with the wavemaker. This field is generated numerically as a Fourier series with random coefficients. The energy spectrum is that in (3.18), with
$v_{\mathit{rms}}=0.25~\text{m}~\text{s}^{-1}$
and
$\unicode[STIX]{x1D705}=1.45\times 10^{-5}~\text{m}^{-1}$
leading to flows with a correlation length of approximately 200 km. The other parameters are those of the mode-1
$M_{2}$
tide at
$45^{\circ }$
, as in § 3.2. The results presented below use a time-independent flow. We carried out additional simulations with slowly time-varying flows to confirm the theoretical prediction of appendix B that IT scattering is essentially unaffected by the time dependence of stationary random flows with time scales
$O(Ro^{-1})$
longer than the IT period. Note that the modelling of the background flow by a Gaussian random field is a choice motivated by practicality and the fact that, according to our theory, the only statistical property of the flow that influences the scattering is its energy spectrum. An alternative would be to carry out a large number of quasigeostrophic simulations to generate an ensemble of flows with more realistic statistics. This would however be computationally expensive; it would also require great care to control the flow parameters and to ensure stationary statistics.
Figure 4. Energy density in the
$\boldsymbol{k}$
plane in different regions along the channel, from an ensemble of 100 shallow-water simulations. The regions are five
$1024~\text{km}\times 1024~\text{km}$
boxes centred about the midpoints
$x=\{1000,2250,3500,4750,6000\}~\text{km}$
, as indicated in figure 3.
Figure 3(a) shows one realisation of the IT height field
$\unicode[STIX]{x1D702}$
at
$t=80$
days, when plane waves generated by the wavemaker (indicated by a dashed line) have propagated across the eddy field shown in the bottom panel. Close to the wavemaker the wave field has a plane-wave structure which becomes scrambled in appearance as the phases randomise due to flow scattering. In agreement with our scattering theory, the wave field retains a single length scale – the wavelength set by the tidal frequency – throughout its evolution: scattering does not lead to a scale cascade. This is consistent with the earlier simulation results of Ward & Dewar (Reference Ward and Dewar2010), Ponte & Klein (Reference Ponte and Klein2015), Dunphy et al. (Reference Dunphy, Ponte, Klein and Le Gentil2017) and Wagner et al. (Reference Wagner, Ferrando and Young2017), the latter two in a three-dimensional set-up. Note that, since we solve the linearised system, there is no harmonic generation and consequent formation of smaller wave scales as described by Ward & Dewar (Reference Ward and Dewar2010). The eigenvalues
$\unicode[STIX]{x1D706}_{n}$
for the IT and flow parameters chosen are those shown in figure 2(a) and correspond to
$L_{scat}=420~\text{km}$
and
$L_{iso}=5600~\text{km}$
. This is qualitatively consistent with the wave field in figure 3.
To assess our theoretical results in more detail, we need to estimate the energy density
$a(\boldsymbol{x},\boldsymbol{k},t)$
from the simulations. To this end, we take Fourier transforms of the wave fields in five
$1024~\text{km}\times 1024~\text{km}$
square boxes spanning the length of the domain. For each realisation of the flow, we compute the Fourier transform of
$u,v$
and
$\unicode[STIX]{x1D702}$
in each box at the end of the simulation, project onto the IT eigenmode then average over the ensemble to obtain an approximation
$\tilde{a}(\boldsymbol{k},t)$
, say, of
$a(\boldsymbol{x},\boldsymbol{k},t)$
in the box. The details of this are given in appendix C. The results are shown in figure 4. For the first box, located immediately to the right of the wavemaker, most of the energy is concentrated at the single point
$\boldsymbol{k}=(k_{0},0)$
, where
$k_{0}=\sqrt{(\unicode[STIX]{x1D714}^{2}-f^{2})/gh}$
, indicating a pure plane wave propagating to the right. (There is also a faint signal of left-propagating waves resulting from scattering at larger
$x$
.) In the next boxes, energy spreads around the circle of constant radius
$|\boldsymbol{k}|=k_{0}$
. The spectrum is in fact distributed over a finite-width annulus rather than a circle, as a result of off-resonant interactions between waves and flow.
We obtain a clearer view of the distribution of energy as a function of
$\unicode[STIX]{x1D703}$
by integrating
$a(\boldsymbol{x},\boldsymbol{k},t)$
in the 90 angular sectors
$2(n-1)\unicode[STIX]{x03C0}/90\leqslant \unicode[STIX]{x1D703}\leqslant 2n\unicode[STIX]{x03C0}/90$
,
$n=1,\ldots ,90$
. The results are shown in figure 5(a). They enable a better assessment of the validity of the length scale estimates
$L_{\mathit{scat}}=420~\text{km}$
and
$L_{\mathit{iso}}=5600~\text{km}$
. Note that these lengths should be measured from the point where the background flow starts, which is at approximately
$x=1000~\text{km}$
, so that we would expect to see the fields isotropise at
$x>6000~\text{km}$
along the channel, corresponding to the rightmost box of figure 4. We see however that the field is not fully isotropic in that region. An explanation is that the numerical simulation includes an absorbing layer near the right boundary of the domain, which allows energy to exit the channel but not to re-enter it. As a result, there are no left-propagating waves at the end of the channel. In addition, we emphasise that
$L_{\mathit{iso}}$
is only an order-of-magnitude estimate which, by converting time scale into length scale using the group speed, ignores the directional properties of the transport of wave energy with the group velocity. We next go beyond this order-of-magnitude estimate and make direct predictions for the scattering by solving the kinetic equation numerically.
Figure 5. Estimated energy density
$\tilde{a}$
as a function of
$\unicode[STIX]{x1D703}$
from (a) the ensemble of shallow-water simulations, with the
$x_{i}$
taken as the box midpoints, and (b) the kinetic-equation simulation. The parameters are those of figure 3 and the time corresponds to the end of the simulation.
4.2 Kinetic-equation simulations
We simulate the kinetic equation (2.6) under the assumption of homogeneity in the
$y$
-direction, consistent with the periodic boundary conditions, and using the angle
$\unicode[STIX]{x1D703}$
, with
$\boldsymbol{k}=k_{0}(\cos \unicode[STIX]{x1D703},\sin \unicode[STIX]{x1D703})$
, as an independent variable. This reduces the number of independent variables to 3, with
$a(x,\unicode[STIX]{x1D703},t)$
, so that the kinetic equation becomes
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}a+c_{g}\cos \unicode[STIX]{x1D703}\unicode[STIX]{x2202}_{x}a=({\mathcal{L}}-\unicode[STIX]{x1D6F4})a+{\mathcal{F}}(x,\unicode[STIX]{x1D703}),\end{eqnarray}$$
where
${\mathcal{L}}$
and
$\unicode[STIX]{x1D6F4}$
are given by (3.2) and (3.5), and
${\mathcal{F}}(x,\unicode[STIX]{x1D703})$
is a forcing term mimicking the wavemaker of the shallow-water simulations. We take
${\mathcal{F}}$
to be a Gaussian centred about
$x=400~\text{km}$
,
$\unicode[STIX]{x1D703}=0$
, with width parameters
$\unicode[STIX]{x1D6E5}_{x}=40~\text{km}$
and
$\unicode[STIX]{x1D6E5}_{\unicode[STIX]{x1D703}}=0.1$
, and an amplitude that is scaled to match the initial energy peak from the shallow-water simulation data.
We simulate (4.2) using a pseudospectral splitting method. The advection term is integrated using a semi-Lagrangian finite-difference scheme, and the scattering and forcing terms on the right-hand side are integrated exactly in Fourier space. The domain is
$7168~\text{km}\times 2\unicode[STIX]{x03C0}$
, with a resolution of
$1792\times 256$
and time step
$\unicode[STIX]{x0394}t=900~\text{s}$
up to a final time of 80 days. We apply periodic boundary conditions in the
$\unicode[STIX]{x1D703}$
-direction, and place absorbing layers 30 grid points wide at each end of the domain in the
$x$
-direction, as in the shallow-water simulation. The evolution of
$a(x,\unicode[STIX]{x1D703},t)$
is illustrated in figure 6. The wave energy, initially concentrated at
$(x,\unicode[STIX]{x1D703})=(400,0)$
, gets advected by the group velocity in the
$x$
-direction and spreads in the angular direction. Once energy reaches
$|\unicode[STIX]{x1D703}|>\unicode[STIX]{x03C0}/2$
, it propagates to the left, leading to the weak signal for
$k<0$
observed in figure 4(a).
Figure 6. Wave-energy density
$a(x,\unicode[STIX]{x1D703},t)$
at
$t=8$
, 16, 32 and 80 days obtained by solving the kinetic equations numerically with parameters matching those of figure 3. The
$x$
-axis is in units of 1000 km.
At the end of the simulation, we evaluate
$a(x,\unicode[STIX]{x1D703},t)$
at different points along the channel to get a set of curves
$a(x_{i},\unicode[STIX]{x1D703},t=80~\text{days})$
,
$1\leqslant i\leqslant 5$
. The points
$x_{i}$
are taken to be spaced along the channel in the same way as the windows shown in figure 4. Note that we have to take into account the fact that there is no flow in the shallow-water simulations from the point of generation at
$x=400~\text{km}$
to approximately
$x=1000~\text{km}$
(see figure 3), whereas the kinetic equation assumes a flow is present throughout the domain. To resolve the discrepancy, the values of
$x_{i}$
are taken 600 km less than the midpoints of the boxes used for the shallow-water simulations. The functions
$a(x_{i},\unicode[STIX]{x1D703})$
are shown next to the equivalent shallow-water estimates in figure 5. There is a remarkably good agreement between the solution of the kinetic equation and the shallow-water simulation results. This demonstrates the value of the kinetic equation in predicting the generic properties of the scattering and their dependence on the various parameters in the problem.
Figure 7. Energy density from 10 realisations of shallow-water simulations evaluated (a) in the second box of figure 4, centred about
$x_{2}=2250~\text{km}$
, and (b) in the fourth box, centred about
$x_{4}=4750~\text{km}$
.
We emphasise that the kinetic equation yields only ensemble average predictions and cannot describe the details of the effect of a single flow realisation on the IT. To illustrate how the scattering fluctuates between realisations, we show in figure 7 the energy density
$\tilde{a}$
estimated from single shallow-water simulations. While fluctuations can be large, the typical behaviour is a redistribution of energy in the angular direction that is well captured by the ensemble-averaged predictions. Furthermore, ergodicity implies that the ensemble average deductions of the kinetic equation apply accurately to quantities that are spatial averages over many eddy scales.
5 Discussion
This paper examines the scattering of oceanic ITs – or indeed of any inertia–gravity wave – caused by the turbulent mesoscale flow in which they propagate. Assuming that the flow is barotropic, weak (small Rossby number) and random with stationary and homogeneous statistics, we derive the kinetic equation (2.6) governing energy exchanges between waves travelling in different directions. A key outcome is the scattering cross-section (2.9), or (3.4) for an isotropic flow, which measures the energy transfer rate as a function of the wavevectors involved and other parameters. The scattering cross-section depends linearly on the energy spectrum of the flow.
The form of the scattering cross-section shows that the energy exchanges between waves are restricted to waves with the same frequency and hence the same horizontal wavenumber. Therefore, while scattering results in a complex random wave field, this field has a single spatial scale determined by the forcing frequency. This is obvious for a time-independent flow but perhaps less so for the time-dependent flows we consider for which it is a consequence of the statistical stationarity of the flows. Note that this does not imply that the wave field is completely phase locked in time: slow phase variations result from the interaction with a time-dependent flow (Ponte & Klein Reference Ponte and Klein2015; Dunphy et al.
Reference Dunphy, Ponte, Klein and Le Gentil2017), but these are not described by our analysis, which focuses on the wave amplitude as measured by the energy density
$a(\boldsymbol{x},\boldsymbol{k},t)$
. It would be of interest to study the phase variations from the statistical viewpoint taken here.
For an isotropic flow, scattering leads to an equilibrium isotropic wave field over a time scale that we can estimate from the Fourier transform of the scattering cross-section (3.4). At equilibrium, and in the absence of spatial modulations, the wave energy density is
$a(\boldsymbol{k})=a(|\boldsymbol{k}|)\propto \unicode[STIX]{x1D6FF}(|\boldsymbol{k}|-k_{0})$
with
$k_{0}=\sqrt{(\unicode[STIX]{x1D714}^{2}-f^{2})/gh}$
, corresponding to a correlation function
$\int a(|\boldsymbol{k}|)\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\,\text{d}\boldsymbol{k}\propto \int \unicode[STIX]{x1D6FF}(|\boldsymbol{k}|-k_{0})\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\,\text{d}\boldsymbol{k}$
that is proportional to the Bessel function
$\text{J}_{0}(k_{0}|\boldsymbol{x}|)$
. Thus, while the flow controls the speed of convergence towards the isotropic wave-energy distribution, it has no effect on the form of this distribution.
The time scale necessary for scattering to significantly alter the wave field is deduced from the scattering cross-section and found to be of the order of a few days to a week. This is short compared with time scales associated with nonlinear wave–wave interactions (see Ward & Dewar Reference Ward and Dewar2010), which raises the possibility that scattering is as crucial as the more widely considered wave–wave interactions in shaping the inertia–gravity wave spectrum in the ocean. Note however that, as our asymptotic treatment makes clear, the large time scale separation between inertia–gravity waves and the quasigeostrophic flow implies that scattering causes little frequency broadening and so cannot by itself explain the continuum of observed frequencies.
The conclusion that scattering simply relaxes wave energy towards an isotropic equilibrium depends crucially on the two-dimensional set-up implied by our assumption of barotropic flow. It holds because scattering redistributes wave energy in Fourier space over constant-frequency curves which, in this case, are just circles, hence compact. This picture changes radically in the presence of vertical shear since this causes energy exchanges between different vertical wavenumbers. The relevant constant-frequency set is then a cone in wavevector space. Because this cone is unbounded, no finite-energy equilibrium exists, and the energy of an initially plane wave can be expected to cascade to small scales, both horizontally and vertically and with a fixed aspect ratio, as it spreads on the cone. This scenario, already envisioned by Lelong & Riley (Reference Lelong and Riley1991) and Bartello (Reference Bartello1995), is potentially important for the dissipation of oceanic inertia–gravity waves, with implications for their impact on mixing and mesoscale dynamics and for the maintenance of balance. The framework we have adopted, which regards the flow as a prescribed random field and examines the wave statistics on the basis of a kinetic equation, generalises to the case of vertically sheared flow. It is well suited to describe and quantify the scale cascade that results from what is then a fully three-dimensional scattering. This is the subject of future work.
Acknowledgements
This research is funded by the UK Natural Environment Research Council under the NSFGEO-NERC programme (grant NE/R006652/1). M.A.C.S. was supported by The Maxwell Institute Graduate School in Analysis and its Applications, a Centre for Doctoral Training funded by the UK Engineering and Physical Sciences Research Council (grant EP/L016508/01), the Scottish Funding Council, Heriot-Watt University and the University of Edinburgh.
Appendix A. Vertical-mode decomposition
Since the background flow in (2.1) is barotropic, we project onto a vertical-mode basis according to
$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}u\\ v\\ w\\ p\\ b\end{array}\right)=\mathop{\sum }_{m=0}^{\infty }\left(\begin{array}{@{}cc@{}}u_{m}(x,y,t) & F_{m}(z)\\ v_{m}(x,y,t) & F_{m}(z)\\ w_{m}(x,y,t) & N^{-2}(z)F_{m}^{\prime }\\ p_{m}(x,y,t) & F_{m}(z)\\ b_{m}(x,y,t) & F_{m}^{\prime }\end{array}\right),\end{eqnarray}$$
where the
$F_{m}$
are eigenfunctions of the Sturm–Liouville problem
$$\begin{eqnarray}\mathscr{L}F_{m}=-\frac{1}{r_{m}^{2}}F_{m},\quad F_{m}^{\prime }(0)+N^{2}F_{m}(0)/g=0,\quad F_{m}^{\prime }(-H)=0,\end{eqnarray}$$
where
$\mathscr{L}(\cdot )=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}z[f^{2}/N^{2}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}z(\cdot )]$
and
$H$
is the ocean depth (Olbers et al.
Reference Olbers, Willebrand and Eden2012). The eigenvalues
$r_{m}$
are the Rossby radii of deformation. Orthonormality implies
$$\begin{eqnarray}(F_{n},F_{m})=\int _{-H}^{0}F_{m}F_{n}\,\text{d}z=\unicode[STIX]{x1D6FF}_{mn}.\end{eqnarray}$$
Substituting (A 1) into (2.1) leads to a system for the amplitudes
$u_{m},v_{m}$
, etc. of each baroclinic mode
$m$
. Defining the equivalent depth
$h_{m}=f^{2}r_{m}^{2}/g$
and the equivalent surface height
$\unicode[STIX]{x1D702}_{m}=b_{m}/g$
, we rewrite this system in the shallow-water-like form given by (2.2).
Appendix B. Kinetic-equation derivation
In this appendix we apply the formalism of Ryzhik et al. (Reference Ryzhik, Papanicolaou and Keller1996) to derive a kinetic equation for the equivalent shallow-water system (2.2) (see also Powell & Vanneste Reference Powell and Vanneste2005; Bal et al.
Reference Bal, Komorowski and Ryzhik2010). This formalism exploits the weakness of the flow as measured by the small Rossby number, as well as the scale separation between wavelength and flow scale on the one hand, and the scale of variation of the wave amplitude on the other hand. The associated small parameter,
$\unicode[STIX]{x1D700}$
, is assumed to satisfy
$\unicode[STIX]{x1D700}=O(Ro^{2})$
.
B.1 Shallow-water equations and scaling regime
We start by rewriting (2.2) in the more abstract notation
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D753}+[\unicode[STIX]{x1D647}(\unicode[STIX]{x2202}_{\boldsymbol{x}})+\sqrt{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D649}(\boldsymbol{x},\unicode[STIX]{x2202}_{\boldsymbol{x}},\sqrt{\unicode[STIX]{x1D700}}t)]\unicode[STIX]{x1D753}=0,\end{eqnarray}$$
where the vector
$$\begin{eqnarray}\unicode[STIX]{x1D753}(\boldsymbol{x},t)=(u,\,v,\,\unicode[STIX]{x1D702})^{\text{T}}\end{eqnarray}$$
groups the dynamical variables and we have introduced the matrix operators
$$\begin{eqnarray}\unicode[STIX]{x1D647}(\unicode[STIX]{x2202}_{x})=\left(\begin{array}{@{}ccc@{}}0 & -f & g\unicode[STIX]{x2202}_{x}\\ f & 0 & g\unicode[STIX]{x2202}_{y}\\ h\unicode[STIX]{x2202}_{x} & h\unicode[STIX]{x2202}_{y} & 0\end{array}\right)\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D649}(\boldsymbol{x},\unicode[STIX]{x2202}_{\boldsymbol{x}},\sqrt{\unicode[STIX]{x1D700}}t)=\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D713}_{x}\unicode[STIX]{x2202}_{y}-\unicode[STIX]{x1D713}_{y}\unicode[STIX]{x2202}_{x}-\unicode[STIX]{x1D713}_{xy} & -\unicode[STIX]{x1D713}_{yy} & 0\\ \unicode[STIX]{x1D713}_{xx} & \unicode[STIX]{x1D713}_{x}\unicode[STIX]{x2202}_{y}-\unicode[STIX]{x1D713}_{y}\unicode[STIX]{x2202}_{x}+\unicode[STIX]{x1D713}_{xy} & 0\\ 0 & 0 & \unicode[STIX]{x1D713}_{x}\unicode[STIX]{x2202}_{y}-\unicode[STIX]{x1D713}_{y}\unicode[STIX]{x2202}_{x}\end{array}\right).\end{eqnarray}$$
In (B 1), we have made the relative importance of the various terms explicit by scaling them with the relevant power of the small parameter
$\unicode[STIX]{x1D700}$
. For convenience, we keep the equations in their dimensional form and treat
$\unicode[STIX]{x1D700}$
as a bookkeeping parameter that can be set to
$1$
at the end of the calculation. Note that we could use
$Ro$
as an alternative to
$\unicode[STIX]{x1D700}$
in this bookkeeping role and, in this manner, avoid fractional powers; using
$\unicode[STIX]{x1D700}$
keeps our derivation in closer correspondence to that of Ryzhik et al. (Reference Ryzhik, Papanicolaou and Keller1996). Note also that the dependence of
$\unicode[STIX]{x1D649}$
on
$\sqrt{\unicode[STIX]{x1D700}}t$
arises through the slow time dependence of
$\unicode[STIX]{x1D713}$
.
The depth-averaged energy density for the shallow-water system without background flow is given by
$$\begin{eqnarray}{\mathcal{E}}(\boldsymbol{x},t)={\textstyle \frac{1}{2}}(h(u^{2}+v^{2})+g\unicode[STIX]{x1D702}^{2})={\textstyle \frac{1}{2}}\langle \,\unicode[STIX]{x1D753},\unicode[STIX]{x1D753}\rangle _{\unicode[STIX]{x1D648}},\end{eqnarray}$$
where we have defined the inner product
$$\begin{eqnarray}\langle \,\boldsymbol{f},\boldsymbol{g}\rangle _{\unicode[STIX]{x1D648}}:=\boldsymbol{f}^{\ast }\unicode[STIX]{x1D648}\boldsymbol{g},\quad \text{with }\unicode[STIX]{x1D648}:=\left(\begin{array}{@{}ccc@{}}h & 0 & 0\\ 0 & h & 0\\ 0 & 0 & g\end{array}\right).\end{eqnarray}$$
The matrix
$\unicode[STIX]{x1D647}$
in (B 3) is known as the dispersion matrix, since its eigenvalues give the dispersion relation. Taking the definition in (B 3) and replacing
$\unicode[STIX]{x2202}_{\boldsymbol{x}}$
by
$\text{i}\boldsymbol{k}$
, we find that the solutions to the eigenvalue equation
$\unicode[STIX]{x1D647}(\text{i}\boldsymbol{k})\boldsymbol{b}(\boldsymbol{k})=\text{i}\unicode[STIX]{x1D714}(\boldsymbol{k})\boldsymbol{b}(\boldsymbol{k})$
are given by
$$\begin{eqnarray}\unicode[STIX]{x1D714}_{0}=0,\quad \unicode[STIX]{x1D714}_{\pm }(\boldsymbol{k})=\pm \sqrt{f^{2}+gh|\boldsymbol{k}|^{2}},\end{eqnarray}$$
which is the usual dispersion relation for the rotating shallow-water model. The three eigenvectors are given by
$$\begin{eqnarray}\boldsymbol{b}_{0}=\frac{1}{g\unicode[STIX]{x1D714}_{\pm }^{2}}\left(\begin{array}{@{}c@{}}-\text{i}gl\\ \text{i}gk\\ f\end{array}\right),\quad \boldsymbol{b}_{\pm }=\frac{1}{\sqrt{2h}|\unicode[STIX]{x1D714}_{\pm }||\boldsymbol{k}|}\left(\begin{array}{@{}c@{}}\unicode[STIX]{x1D714}_{\pm }k+\text{i}fl\\ \unicode[STIX]{x1D714}_{\pm }l-\text{i}fk\\ h|\boldsymbol{k}|^{2}\end{array}\right).\end{eqnarray}$$
These are orthonormal in the sense that
$$\begin{eqnarray}\langle \,\boldsymbol{b}_{i},\boldsymbol{b}_{j}\rangle _{\unicode[STIX]{x1D648}}=\unicode[STIX]{x1D6FF}_{ij},\quad i,j=0,\pm .\end{eqnarray}$$
The zero eigenvalue
$\unicode[STIX]{x1D714}_{0}$
corresponds to the vortical mode of the system. Since this mode is accounted for by the quasigeostrophic background flow, it does not appear in any of the following derivation. In addition to the right eigenvectors
$\boldsymbol{b}_{i}$
, we also use the left eigenvectors
$\boldsymbol{c}_{i}$
which satisfy
$$\begin{eqnarray}\boldsymbol{c}_{j}\unicode[STIX]{x1D647}=\text{i}\unicode[STIX]{x1D714}_{j}\boldsymbol{c}_{j},\quad \boldsymbol{c}_{j}=\boldsymbol{b}_{j}^{\ast }\unicode[STIX]{x1D648}\quad \text{and}\quad \boldsymbol{c}_{i}\boldsymbol{b}_{j}=\unicode[STIX]{x1D6FF}_{ij}.\end{eqnarray}$$
B.2 Scaled Wigner transform
It is convenient to rescale the space and time coordinates as
$(\boldsymbol{x},t)\mapsto (\boldsymbol{x}/\unicode[STIX]{x1D700},t/\unicode[STIX]{x1D700})$
, and define
$\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}(\boldsymbol{x},t):=\unicode[STIX]{x1D753}(\boldsymbol{x}/\unicode[STIX]{x1D700},t/\unicode[STIX]{x1D700})$
. Under this rescaling, the equations of motion (B 1) take the form
$$\begin{eqnarray}\unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}+[\unicode[STIX]{x1D647}(\unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{\boldsymbol{x}})+\sqrt{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D649}(\boldsymbol{x}/\unicode[STIX]{x1D700},\unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{\boldsymbol{x}},t/\sqrt{\unicode[STIX]{x1D700}})]\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}=0.\end{eqnarray}$$
In the new coordinates, the wavevector is
$O(\unicode[STIX]{x1D700}^{-1})$
. To account for this and to obtain a well-defined Wigner transform in the limit
$\unicode[STIX]{x1D700}\rightarrow 0$
, we use the rescaled definition
$$\begin{eqnarray}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}(\boldsymbol{x},\boldsymbol{k},t)=\unicode[STIX]{x1D700}^{-2}\unicode[STIX]{x1D652}(\boldsymbol{x},\boldsymbol{k}/\unicode[STIX]{x1D700},t)=\frac{1}{(2\unicode[STIX]{x03C0})^{2}}\int _{\mathbb{R}^{2}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{y}}\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}(\boldsymbol{x}-\unicode[STIX]{x1D700}\boldsymbol{y}/2,t)\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}^{\ast }(\boldsymbol{x}+\unicode[STIX]{x1D700}\boldsymbol{y}/2,t)\,\text{d}\boldsymbol{y}.\end{eqnarray}$$
We note that this function remains positive semi-definite as
$\unicode[STIX]{x1D700}$
goes to zero, which is important as it is closely related to the energy of the system. Indeed it can readily be shown that
$$\begin{eqnarray}\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}(\boldsymbol{x},\boldsymbol{k},t)\,\text{d}\boldsymbol{k}=\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}(\boldsymbol{x},t)\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}^{\ast }(\boldsymbol{x},t),\end{eqnarray}$$
and so the energy density is found from the Wigner transform of the wave fields as
$$\begin{eqnarray}{\mathcal{E}}(\boldsymbol{x},t)=\frac{1}{2}\int _{\mathbb{R}^{2}}\text{tr}(\unicode[STIX]{x1D648}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}(\boldsymbol{x},\boldsymbol{k},t))\,\text{d}\boldsymbol{k}.\end{eqnarray}$$
A useful alternative to (B 12) expresses the Wigner transform in terms of Fourier transforms as
$$\begin{eqnarray}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}(\boldsymbol{x},\boldsymbol{k},t)=\unicode[STIX]{x1D700}^{-2}\int _{\mathbb{R}^{2}}\text{e}^{\text{i}\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{x}}\hat{\unicode[STIX]{x1D753}}_{\unicode[STIX]{x1D700}}\left(-\frac{\boldsymbol{k}}{\unicode[STIX]{x1D700}}-\frac{\boldsymbol{p}}{2},t\right)\hat{\unicode[STIX]{x1D753}}_{\unicode[STIX]{x1D700}}^{\ast }\left(-\frac{\boldsymbol{k}}{\unicode[STIX]{x1D700}}+\frac{\boldsymbol{p}}{2},t\right)\,\text{d}\boldsymbol{p},\end{eqnarray}$$
where
$$\begin{eqnarray}\hat{\unicode[STIX]{x1D753}}_{\unicode[STIX]{x1D700}}(\boldsymbol{k},t)=\frac{1}{(2\unicode[STIX]{x03C0})^{2}}\int _{\mathbb{R}^{2}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}(\boldsymbol{x},t)\,\text{d}\boldsymbol{x}\quad \text{and}\quad \unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}(\boldsymbol{x},t)=\int _{\mathbb{R}^{2}}\text{e}^{-\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}\hat{\unicode[STIX]{x1D753}}_{\unicode[STIX]{x1D700}}(\boldsymbol{k},t)\,\text{d}\boldsymbol{k}.\end{eqnarray}$$
B.3 Wigner-transform evolution equation
We obtain an evolution equation for the Wigner transform by differentiating (B 12) with respect to
$t$
and substituting (B 11). This gives
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}(\boldsymbol{x},\boldsymbol{k},t)\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{\unicode[STIX]{x1D700}}{(2\unicode[STIX]{x03C0})^{2}}\int _{\mathbb{R}^{2}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{y}}\left(\left[\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}\left(t,\boldsymbol{x}-\frac{\unicode[STIX]{x1D700}\boldsymbol{y}}{2}\right)\right]\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}^{\ast }\left(t,\boldsymbol{x}+\frac{\unicode[STIX]{x1D700}\boldsymbol{y}}{2}\right)+\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}\left(t,\boldsymbol{x}-\frac{\unicode[STIX]{x1D700}\boldsymbol{y}}{2}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \qquad \times \left.\left[\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}^{\ast }\left(t,\boldsymbol{x}+\frac{\unicode[STIX]{x1D700}\boldsymbol{y}}{2}\right)\right]\right)\,\text{d}\boldsymbol{y}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{-1}{(2\unicode[STIX]{x03C0})^{2}}\int _{\mathbb{R}^{2}}\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{y}}\left[\unicode[STIX]{x1D647}(\unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{\boldsymbol{x}})+\sqrt{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D649}\left(\frac{\boldsymbol{x}}{\unicode[STIX]{x1D700}}-\frac{\boldsymbol{y}}{2},\unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{\boldsymbol{x}}\right)\right]\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}\left(t,\boldsymbol{x}-\frac{\unicode[STIX]{x1D700}\boldsymbol{y}}{2}\right)\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}^{\ast }\left(t,\boldsymbol{x}+\frac{\unicode[STIX]{x1D700}\boldsymbol{y}}{2}\right)\,\text{d}\boldsymbol{y}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\text{c.c.},\end{eqnarray}$$
where c.c. denotes the complex conjugate of the term preceding it. Rewriting
$\unicode[STIX]{x1D753}_{\unicode[STIX]{x1D700}}$
and
$\unicode[STIX]{x1D649}$
in terms of their Fourier transforms and making use of (B 15), we find that
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D700}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}(\boldsymbol{x},\boldsymbol{k},t)+\overbrace{\left[\unicode[STIX]{x1D647}\left(\text{i}\boldsymbol{k}+\frac{\unicode[STIX]{x1D700}}{2}\unicode[STIX]{x2202}_{\boldsymbol{x}}\right)\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}(\boldsymbol{x},\boldsymbol{k})+\text{c.c.}\right]}^{:={\mathcal{Q}}^{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\sqrt{\unicode[STIX]{x1D700}}\left[\,\underbrace{\int \text{e}^{-\text{i}\boldsymbol{p}\boldsymbol{\cdot }\unicode[STIX]{x1D743}}\widehat{\unicode[STIX]{x1D649}}\left(\boldsymbol{p},\text{i}\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)+\frac{\unicode[STIX]{x1D700}}{2}\unicode[STIX]{x2202}_{\boldsymbol{x}},\frac{t}{\sqrt{\unicode[STIX]{x1D700}}}\right)\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}\left(\boldsymbol{x},\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\,\text{d}\boldsymbol{p}+\text{c.c.}}_{:={\mathcal{P}}^{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}}\,\right]=0.\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
Scattering effects due to interactions of waves with the background flow are controlled by the third term,
${\mathcal{P}}^{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}$
. To derive it, we have introduced the fast-space variable
$\unicode[STIX]{x1D743}:=\boldsymbol{x}/\unicode[STIX]{x1D700}$
and the Fourier transform
$\widehat{\unicode[STIX]{x1D649}}$
of
$\unicode[STIX]{x1D649}$
defined by
$$\begin{eqnarray}\unicode[STIX]{x1D649}(\unicode[STIX]{x1D743},\cdot )=\int _{\mathbb{R}^{2}}\text{e}^{-\text{i}\boldsymbol{p}\boldsymbol{\cdot }\unicode[STIX]{x1D743}}\widehat{\unicode[STIX]{x1D649}}(\,\boldsymbol{p},\cdot )\,\text{d}\boldsymbol{p}.\end{eqnarray}$$
We now derive the asymptotic limit of (B 18) using a multiscale expansion. Defining the intermediate time variable
$\unicode[STIX]{x1D70F}:=t/\sqrt{\unicode[STIX]{x1D700}}$
to cater for the time dependence of the streamfunction, we expand
$$\begin{eqnarray}\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k},t,\unicode[STIX]{x1D70F})=\unicode[STIX]{x1D652}^{(0)}(\boldsymbol{x},\boldsymbol{k},t)+\sqrt{\unicode[STIX]{x1D700}}\unicode[STIX]{x1D652}^{(1)}(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k},t,\unicode[STIX]{x1D70F})+\unicode[STIX]{x1D700}\unicode[STIX]{x1D652}^{(2)}(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k},t,\unicode[STIX]{x1D70F})+O(\unicode[STIX]{x1D700}^{3/2}),\end{eqnarray}$$
where we have anticipated that the leading-order term depends on the slow variables only. The differential operators are then expanded as
$$\begin{eqnarray}\unicode[STIX]{x2202}_{\boldsymbol{x}}\mapsto \unicode[STIX]{x2202}_{\boldsymbol{x}}+\unicode[STIX]{x1D700}^{-1}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D743}},\quad \unicode[STIX]{x2202}_{t}\mapsto \unicode[STIX]{x2202}_{t}+\unicode[STIX]{x1D700}^{-1/2}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}},\end{eqnarray}$$
where
$\boldsymbol{x}$
and
$\unicode[STIX]{x1D743}$
,
$t$
and
$\unicode[STIX]{x1D70F}$
are treated as independent variables, leading to the expansion
$$\begin{eqnarray}{\mathcal{Q}}^{\unicode[STIX]{x1D700}}={\mathcal{Q}}_{0}+\unicode[STIX]{x1D700}{\mathcal{Q}}_{1}+O(\unicode[STIX]{x1D700}^{2}),\quad {\mathcal{P}}^{\unicode[STIX]{x1D700}}={\mathcal{P}}_{0}+\unicode[STIX]{x1D700}{\mathcal{P}}_{1}+O(\unicode[STIX]{x1D700}^{2})\end{eqnarray}$$
of the operators in (B 18). It turns out that only the leading-order term
${\mathcal{P}}_{0}$
is required for the derivation of the kinetic equation.
The operators in (B 22) can be written explicitly through their action on an arbitrary function
$Z(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k})$
:
$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{{\mathcal{Q}}}_{0}Z(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k})=\unicode[STIX]{x1D647}\left(\text{i}\boldsymbol{k}+{\textstyle \frac{1}{2}}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D743}}\right)Z(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k})+\text{c.c.} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{{\mathcal{Q}}}_{1}Z(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k})={\displaystyle \frac{1}{2\text{i}}}\left[\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D647}\left(\text{i}\boldsymbol{k}+{\displaystyle \frac{1}{2}}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D743}}\right)\right]\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}Z(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k})+\text{c.c.} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \widetilde{{\mathcal{P}}}_{0}Z(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k})=\displaystyle \int _{\mathbb{R}^{2}}\text{e}^{-\text{i}\boldsymbol{p}\boldsymbol{\cdot }\unicode[STIX]{x1D743}}\widehat{\unicode[STIX]{x1D649}}\left(\boldsymbol{p},\text{i}\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)+\frac{1}{2}\unicode[STIX]{x2202}_{\unicode[STIX]{x1D743}},\unicode[STIX]{x1D70F}\right)Z\left(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\,\text{d}\boldsymbol{p}+\text{c.c.}\qquad & \displaystyle\end{eqnarray}$$
We have decorated the operators with a tilde to highlight the presence of
$\unicode[STIX]{x2202}_{\unicode[STIX]{x1D743}}$
in their definition; the tildes will be removed whenever this dependence disappears.
Substituting the operators into (B 18) gives us the evolution equation for the Wigner function as
$$\begin{eqnarray}\left[\frac{1}{\unicode[STIX]{x1D700}}\widetilde{{\mathcal{Q}}}_{0}+\frac{1}{\sqrt{\unicode[STIX]{x1D700}}}\left(\widetilde{{\mathcal{P}}}_{0}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70F}}\right)+\left(\widetilde{{\mathcal{Q}}}_{1}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\right)\right]\unicode[STIX]{x1D652}^{\unicode[STIX]{x1D700}}(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k},t,\unicode[STIX]{x1D70F})=0.\end{eqnarray}$$
Introducing the expansion (B 20) then leads to a hierarchy of equations to be solved at each order in
$\unicode[STIX]{x1D700}$
.
The leading-order equation is
$$\begin{eqnarray}{\mathcal{Q}}_{0}\unicode[STIX]{x1D652}^{(0)}=\unicode[STIX]{x1D647}(\text{i}\boldsymbol{k})\unicode[STIX]{x1D652}^{(0)}(\boldsymbol{x},\boldsymbol{k},t)+\text{c.c.}=0.\end{eqnarray}$$
The eigenvalues of
$\unicode[STIX]{x1D647}$
are purely imaginary, so this equation is satisfied by taking
$\unicode[STIX]{x1D652}^{(0)}$
in the form of a linear combination of the eigenvectors of the dispersion matrix. Defining the matrices
$$\begin{eqnarray}\unicode[STIX]{x1D63D}_{j}(\boldsymbol{k})=\boldsymbol{b}_{j}(\boldsymbol{k})\boldsymbol{b}_{j}^{\ast }(\boldsymbol{k}),\end{eqnarray}$$
the leading order Wigner function is thus given by
$$\begin{eqnarray}\unicode[STIX]{x1D652}^{(0)}(\boldsymbol{x},\boldsymbol{k},t)=\mathop{\sum }_{j=\pm }a_{j}(\boldsymbol{x},\boldsymbol{k},t)\unicode[STIX]{x1D63D}_{j}(\boldsymbol{k}).\end{eqnarray}$$
The so-far undetermined amplitudes
$a_{j}(\boldsymbol{x},\boldsymbol{k},t)$
are real because the Wigner function is Hermitian.
At
$O(\unicode[STIX]{x1D700}^{-1/2})$
, we find
$$\begin{eqnarray}\widetilde{{\mathcal{Q}}}_{0}\unicode[STIX]{x1D652}^{(1)}(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k},t,\unicode[STIX]{x1D70F})=-{\mathcal{P}}_{0}\unicode[STIX]{x1D652}^{(0)}(\boldsymbol{x},\boldsymbol{k},t),\end{eqnarray}$$
where we have used that
$\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\unicode[STIX]{x1D652}^{(0)}=0$
. To solve (B 30), we rewrite
$\unicode[STIX]{x1D652}^{(1)}$
in terms of its Fourier transform with respect to
$\unicode[STIX]{x1D743}$
,
$$\begin{eqnarray}\unicode[STIX]{x1D652}^{(1)}(\boldsymbol{x},\unicode[STIX]{x1D743},\boldsymbol{k},t,\unicode[STIX]{x1D70F})=\int _{\mathbb{ R}^{2}}\text{e}^{-\text{i}\boldsymbol{p}\boldsymbol{\cdot }\unicode[STIX]{x1D743}}\widehat{\unicode[STIX]{x1D652}}^{(1)}(\boldsymbol{x},\boldsymbol{p},\boldsymbol{k},t,\unicode[STIX]{x1D70F})\,\text{d}\boldsymbol{p}.\end{eqnarray}$$
Substituting this into (B 30) yields
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D647}\left(\text{i}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\right)\widehat{\unicode[STIX]{x1D652}}^{(1)}(\,\boldsymbol{p},\boldsymbol{k})+\left[\unicode[STIX]{x1D647}\left(\text{i}\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\right)\widehat{\unicode[STIX]{x1D652}}^{(1)}(-\boldsymbol{p},\boldsymbol{k})\right]^{\ast }+\unicode[STIX]{x1D703}\widehat{\unicode[STIX]{x1D652}}^{(1)}(\,\boldsymbol{p},\boldsymbol{k})\nonumber\\ \displaystyle & & \displaystyle \quad =-\widehat{\unicode[STIX]{x1D649}}\left(\boldsymbol{p},\text{i}\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\right)\unicode[STIX]{x1D652}^{(0)}\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)-\left[\widehat{\unicode[STIX]{x1D649}}\left(-\boldsymbol{p},\text{i}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\right)\unicode[STIX]{x1D652}^{(0)}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\right]^{\ast },\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
where we have suppressed dependencies on
$\boldsymbol{x}$
,
$t$
and
$\unicode[STIX]{x1D70F}$
for conciseness. Following Ryzhik et al. (Reference Ryzhik, Papanicolaou and Keller1996), we have introduced a regularisation parameter
$\unicode[STIX]{x1D703}>0$
which will be taken to zero at a later stage.
Since the Wigner transform is Hermitian,
$\widehat{\unicode[STIX]{x1D652}}^{(1)}(\,\boldsymbol{p},\boldsymbol{k})=\widehat{\unicode[STIX]{x1D652}}^{(1)\ast }(-\boldsymbol{p},\boldsymbol{k})$
. Using this, expanding
$\unicode[STIX]{x1D652}^{(0)}$
according to (B 29), and pre- and post-multiplying the resulting expression by
$\boldsymbol{c}_{n}(\boldsymbol{k}-\boldsymbol{p}/2)$
and
$\boldsymbol{c}_{m}^{\ast }(\boldsymbol{k}+\boldsymbol{p}/2)$
(with
$\boldsymbol{c}_{n}$
the left eigenvector defined in (B 10)) gives
$$\begin{eqnarray}\displaystyle & & \displaystyle -\left(\text{i}\left(\unicode[STIX]{x1D714}_{n}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)-\unicode[STIX]{x1D714}_{m}\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\right)+\unicode[STIX]{x1D703}\right)\boldsymbol{c}_{n}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\widehat{\unicode[STIX]{x1D652}}^{(1)}(\,\boldsymbol{p},\boldsymbol{k})\boldsymbol{c}_{m}^{\ast }\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\nonumber\\ \displaystyle & & \displaystyle \quad =\mathop{\sum }_{i=\pm }a_{i}\left(\boldsymbol{k}\!+\frac{\boldsymbol{p}}{2}\right)\boldsymbol{c}_{n}\left(\boldsymbol{k}\!-\frac{\boldsymbol{p}}{2}\right)\widehat{\unicode[STIX]{x1D649}}\left(\boldsymbol{p},\text{i}\left(\boldsymbol{k}\!+\frac{\boldsymbol{p}}{2}\right)\right)\boldsymbol{b}_{i}\left(\boldsymbol{k}\!+\frac{\boldsymbol{p}}{2}\right)\boldsymbol{b}_{i}^{\ast }\left(\boldsymbol{k}\!+\frac{\boldsymbol{p}}{2}\right)\boldsymbol{c}_{m}^{\ast }\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\mathop{\sum }_{j=\pm }a_{j}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\boldsymbol{c}_{n}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\boldsymbol{b}_{j}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\boldsymbol{b}_{j}^{\ast }\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\widehat{\unicode[STIX]{x1D649}}^{\ast }\left(-\boldsymbol{p},\text{i}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\right)\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,\boldsymbol{c}_{m}^{\ast }\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right).\end{eqnarray}$$
It is convenient to extract the linear dependence of
$\widehat{\unicode[STIX]{x1D649}}$
on the streamfunction by defining a matrix
$\widehat{\unicode[STIX]{x1D650}}(\,\boldsymbol{p},\text{i}\boldsymbol{q})$
such that
$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D649}}(\,\boldsymbol{p},\text{i}\boldsymbol{q},\unicode[STIX]{x1D70F})=\widehat{\unicode[STIX]{x1D650}}(\,\boldsymbol{p},\text{i}\boldsymbol{q})\widehat{\unicode[STIX]{x1D713}}(\,\boldsymbol{p},\unicode[STIX]{x1D70F}).\end{eqnarray}$$
We now decompose
$\widehat{\unicode[STIX]{x1D652}}^{(1)}$
using the vectors
$\boldsymbol{b}_{i}(\boldsymbol{k})$
, which form a complete basis, as
$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D652}}^{(1)}(\boldsymbol{x},\boldsymbol{p},\boldsymbol{k},t,\unicode[STIX]{x1D70F})=\mathop{\sum }_{m,n=\pm }\unicode[STIX]{x1D6FC}_{mn}(\boldsymbol{x},\boldsymbol{p},\boldsymbol{k},t,\unicode[STIX]{x1D70F})\boldsymbol{b}_{n}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\boldsymbol{b}_{m}^{\ast }\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right).\end{eqnarray}$$
Using this along with the orthonormality of the eigenvectors we finally write the solution
$$\begin{eqnarray}\displaystyle & & \displaystyle \widehat{\unicode[STIX]{x1D652}}^{(1)}(\boldsymbol{x},\boldsymbol{p},\boldsymbol{k},t,\unicode[STIX]{x1D70F})=\mathop{\sum }_{m,n=\pm }\left[a_{m}\left(\boldsymbol{x},\boldsymbol{k}+\frac{\boldsymbol{p}}{2},t\right)\boldsymbol{c}_{n}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\widehat{\unicode[STIX]{x1D650}}\left(\boldsymbol{p},\text{i}\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\right)\boldsymbol{b}_{m}\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \quad +\left.a_{n}\left(\boldsymbol{x},\boldsymbol{k}-\frac{\boldsymbol{p}}{2},t\right)\boldsymbol{b}_{n}^{\ast }\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\widehat{\unicode[STIX]{x1D650}}^{\ast }\left(-\boldsymbol{p},\text{i}\left(\boldsymbol{k}-\frac{\boldsymbol{p}}{2}\right)\right)\boldsymbol{c}_{m}^{\ast }\left(\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\right]\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\frac{\boldsymbol{b}_{n}\left(\boldsymbol{k}-{\displaystyle \frac{\boldsymbol{p}}{2}}\right)\boldsymbol{b}_{m}^{\ast }\left(\boldsymbol{k}+{\displaystyle \frac{\boldsymbol{p}}{2}}\right)\widehat{\unicode[STIX]{x1D713}}(\,\boldsymbol{p},\unicode[STIX]{x1D70F})}{\text{i}\left(\unicode[STIX]{x1D714}_{m}\left(\boldsymbol{k}+{\displaystyle \frac{\boldsymbol{p}}{2}}\right)-\unicode[STIX]{x1D714}_{n}\left(\boldsymbol{k}-{\displaystyle \frac{\boldsymbol{p}}{2}}\right)\right)-\unicode[STIX]{x1D703}},\end{eqnarray}$$
where we have taken into account that
$\widehat{\unicode[STIX]{x1D713}}(\,\boldsymbol{p})=\widehat{\unicode[STIX]{x1D713}}^{\ast }(-\boldsymbol{p})$
. We note that this solution shows
$\unicode[STIX]{x1D652}^{(1)}$
is linear in the random field
$\unicode[STIX]{x1D713}$
.
The slow evolution of the leading-order Wigner function
$\unicode[STIX]{x1D652}^{(0)}$
is controlled by the
$O(1)$
term in the expansion of (B 26), given by
$$\begin{eqnarray}-\widetilde{{\mathcal{Q}}}_{0}\unicode[STIX]{x1D652}^{(2)}=(\widetilde{{\mathcal{P}}}_{0}+\unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}})\unicode[STIX]{x1D652}^{(1)}+({\mathcal{Q}}_{1}+\unicode[STIX]{x2202}_{t})\unicode[STIX]{x1D652}^{(0)}.\end{eqnarray}$$
We assume that the random streamfunction is a stationary process in
$\unicode[STIX]{x1D70F}$
and homogeneous in
$\unicode[STIX]{x1D743}$
, with zero mean,
$\langle \unicode[STIX]{x1D713}(\unicode[STIX]{x1D743},\unicode[STIX]{x1D70F})\rangle =0$
, and covariance
$$\begin{eqnarray}\langle \unicode[STIX]{x1D713}(\unicode[STIX]{x1D743},\unicode[STIX]{x1D70F})\unicode[STIX]{x1D713}(\unicode[STIX]{x1D743}^{\prime },\unicode[STIX]{x1D70F})\rangle =R(\unicode[STIX]{x1D743}-\unicode[STIX]{x1D743}^{\prime }),\end{eqnarray}$$
where
$\langle \cdot \rangle$
denotes an ensemble average, or equivalently an average over
$\unicode[STIX]{x1D743}$
. In terms of Fourier transforms, this implies that
$$\begin{eqnarray}\langle \widehat{\unicode[STIX]{x1D713}}(\,\boldsymbol{p})\widehat{\unicode[STIX]{x1D713}}(\,\boldsymbol{p}^{\prime })\rangle =\widehat{R}(\,\boldsymbol{p})\,\unicode[STIX]{x1D6FF}(\,\boldsymbol{p}+\boldsymbol{p}^{\prime }),\end{eqnarray}$$
where the streamfunction power spectrum
$\widehat{R}$
is the Fourier transform of
$R$
. The more familiar kinetic-energy spectrum is then
$$\begin{eqnarray}\widehat{E}(\boldsymbol{k})=|\boldsymbol{k}|^{2}\widehat{R}(\boldsymbol{k}).\end{eqnarray}$$
We now take the average of (B 37). The slow time derivative term on the right-hand side disappears since
$\langle \unicode[STIX]{x1D652}^{(1)}\rangle =0$
. Since
$\langle \unicode[STIX]{x2202}_{\unicode[STIX]{x1D709}}\unicode[STIX]{x1D652}^{(2)}\rangle =0$
,
$\langle \widetilde{{\mathcal{Q}}}_{0}\unicode[STIX]{x1D652}^{(2)}\rangle ={\mathcal{Q}}_{0}\langle \unicode[STIX]{x1D652}^{(2)}\rangle$
, where the removal of the tilde corresponds to setting
$\unicode[STIX]{x2202}_{\unicode[STIX]{x1D743}}$
to 0 in
${\mathcal{Q}}_{0}$
. This leads to
$$\begin{eqnarray}-{\mathcal{Q}}_{0}\unicode[STIX]{x1D652}^{(2)}=\langle \widetilde{{\mathcal{P}}}_{0}\unicode[STIX]{x1D652}^{(1)}+({\mathcal{Q}}_{1}+\unicode[STIX]{x2202}_{t})\unicode[STIX]{x1D652}^{(0)}\rangle ,\end{eqnarray}$$
an inhomogeneous version of (B 27).
The matrix
${\mathcal{Q}}_{0}$
has a non-trivial null space, spanned by the matrices
$\unicode[STIX]{x1D63D}_{j}(\boldsymbol{k})$
; the right-hand side of (B 41) must therefore satisfy a solvability condition. Since
$\text{i}{\mathcal{Q}}_{0}$
is self-adjoint with respect to the matrix inner product
$$\begin{eqnarray}\langle \!\langle \unicode[STIX]{x1D653},\unicode[STIX]{x1D654}\rangle \!\rangle :=\text{tr}(\unicode[STIX]{x1D648}\unicode[STIX]{x1D653}^{\ast }\unicode[STIX]{x1D648}\unicode[STIX]{x1D654}),\end{eqnarray}$$
this condition is obtained by applying
$\langle \!\langle \unicode[STIX]{x1D63D}_{j}(\boldsymbol{k}),\cdot \rangle \!\rangle$
to (B 41). We deal with the resulting terms one by one. First, by orthogonality and (B 29) we have
$$\begin{eqnarray}\langle \!\langle \unicode[STIX]{x1D63D}_{i},\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D652}^{(0)}\rangle \!\rangle =\mathop{\sum }_{j=\pm }(\unicode[STIX]{x2202}_{t}a_{j})\langle \!\langle \unicode[STIX]{x1D63D}_{i},\unicode[STIX]{x1D63D}_{j}\rangle \!\rangle =\unicode[STIX]{x2202}_{t}a_{i}(\boldsymbol{x},\boldsymbol{k},t).\end{eqnarray}$$
Next,
$$\begin{eqnarray}\displaystyle \langle \!\langle \unicode[STIX]{x1D63D}_{i},{\mathcal{Q}}_{1}\unicode[STIX]{x1D652}^{(0)}\rangle \!\rangle & = & \displaystyle \mathop{\sum }_{j=\pm }\frac{1}{2\text{i}}\langle \!\langle \unicode[STIX]{x1D63D}_{i},(\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D647}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a_{j})\unicode[STIX]{x1D63D}_{j}\rangle \!\rangle +\text{c.c.}\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{j=\pm }\frac{1}{2\text{i}}\langle \!\langle \unicode[STIX]{x1D63D}_{i},\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}(\unicode[STIX]{x1D647}\unicode[STIX]{x1D63D}_{j})-\unicode[STIX]{x1D647}\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D63D}_{j}\rangle \!\rangle \boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a_{j}+\text{c.c.}\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{j=\pm }\frac{1}{2\text{i}}\langle \!\langle \unicode[STIX]{x1D63D}_{i},\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}(\text{i}\unicode[STIX]{x1D714}_{j})\unicode[STIX]{x1D63D}_{j}+(\text{i}\unicode[STIX]{x1D714}_{j}-\unicode[STIX]{x1D647})\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D63D}_{j}\rangle \!\rangle \boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a_{j}+\text{c.c.}\nonumber\\ \displaystyle & = & \displaystyle \unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a_{i}(\boldsymbol{x},\boldsymbol{k},t).\end{eqnarray}$$
In order to evaluate the remaining term, we note that using (B 34) and (B 39), we have
$$\begin{eqnarray}\langle \widehat{\unicode[STIX]{x1D615}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\,\boldsymbol{p},\text{i}\boldsymbol{q})\widehat{\unicode[STIX]{x1D615}}_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}(\,\boldsymbol{p}^{\prime },\text{i}\boldsymbol{q}^{\prime })\rangle =\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\,\boldsymbol{p},\text{i}\boldsymbol{q})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D6FF}}(\,\boldsymbol{p}^{\prime },\text{i}\boldsymbol{q}^{\prime })\widehat{R}(\,\boldsymbol{p})\,\unicode[STIX]{x1D6FF}(\,\boldsymbol{p}+\boldsymbol{p}^{\prime }),\end{eqnarray}$$
where Greek indices are used for matrix elements to make the following derivation clearer, and summation over repeated Greek indices is implied.
Expanding all terms, and making use of the delta function in (B 45), we have
$$\begin{eqnarray}\displaystyle & & \displaystyle \langle \!\langle \unicode[STIX]{x1D63D}_{i},\langle \widetilde{{\mathcal{P}}}_{0}\unicode[STIX]{x1D652}^{(1)}\rangle \rangle \!\rangle \nonumber\\ \displaystyle & & \displaystyle =-\iint \text{e}^{\text{i}(\,\boldsymbol{p}+\boldsymbol{p}^{\prime })\boldsymbol{\cdot }\unicode[STIX]{x1D743}}\unicode[STIX]{x1D614}_{\unicode[STIX]{x1D708}\unicode[STIX]{x1D70C}}b_{\unicode[STIX]{x1D70C}}^{i}(\boldsymbol{k})b_{\unicode[STIX]{x1D70E}}^{i\ast }(\boldsymbol{k})\unicode[STIX]{x1D614}_{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D706}}\left\langle \widehat{\unicode[STIX]{x1D615}}_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D707}}\left(\boldsymbol{p},\text{i}\left(\boldsymbol{k}+\frac{\boldsymbol{p}-\boldsymbol{p}^{\prime }}{2}\right)\right)\right.\nonumber\\ \displaystyle & & \displaystyle \quad \times \left.\widehat{\unicode[STIX]{x1D61E}}_{\unicode[STIX]{x1D707}\unicode[STIX]{x1D708}}^{(1)}\left(\boldsymbol{p}^{\prime },\boldsymbol{k}+\frac{\boldsymbol{p}}{2}\right)\vphantom{\left(\boldsymbol{k}+\frac{\boldsymbol{p}-\boldsymbol{p}^{\prime }}{2}\right)}\right\rangle \,\text{d}\boldsymbol{p}\,\text{d}\boldsymbol{p}^{\prime }+\text{c.c.}\nonumber\\ \displaystyle & & \displaystyle =-\int \mathop{\sum }_{m,n=\pm }c_{\unicode[STIX]{x1D706}}^{i}(\boldsymbol{k})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D707}}(\,\boldsymbol{p},\text{i}(\boldsymbol{k}+\boldsymbol{p}))b_{\unicode[STIX]{x1D707}}^{n}(\boldsymbol{k}+\boldsymbol{p})\overbrace{c_{\unicode[STIX]{x1D70C}}^{m}(\boldsymbol{k})b_{\unicode[STIX]{x1D70C}}^{i}(\boldsymbol{k})}^{\unicode[STIX]{x1D6FF}^{im}}\widehat{R}(\,\boldsymbol{p})\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\frac{a_{m}(\boldsymbol{k})c_{\unicode[STIX]{x1D6FC}}^{n}(\boldsymbol{k}+\boldsymbol{p})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(-\boldsymbol{p},\text{i}\boldsymbol{k})b_{\unicode[STIX]{x1D6FD}}^{m}(\boldsymbol{k})+a_{n}(\boldsymbol{k}+\boldsymbol{p})(c_{\unicode[STIX]{x1D6FC}}^{m}(\boldsymbol{k})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\,\boldsymbol{p},\text{i}(\boldsymbol{k}+\boldsymbol{p}))b_{\unicode[STIX]{x1D6FD}}^{n}(\boldsymbol{k}+\boldsymbol{p}))^{\ast }}{\text{i}(\unicode[STIX]{x1D714}_{m}(\boldsymbol{k})-\unicode[STIX]{x1D714}_{n}(\boldsymbol{k}+\boldsymbol{p}))-\unicode[STIX]{x1D703}}\,\text{d}\boldsymbol{p}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\text{c.c.}\nonumber\\ \displaystyle & & \displaystyle \quad =2\unicode[STIX]{x1D703}\text{Re}\,\int \mathop{\sum }_{n=\pm }c_{\unicode[STIX]{x1D706}}^{i}(\boldsymbol{k})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D707}}(\boldsymbol{k}^{\prime }-\boldsymbol{k},\text{i}\boldsymbol{k}^{\prime })b_{\unicode[STIX]{x1D707}}^{n}(\boldsymbol{k}^{\prime })\widehat{R}(\boldsymbol{k}^{\prime }-\boldsymbol{k})\nonumber\\ \displaystyle & & \displaystyle \quad \times \,\frac{a_{i}(\boldsymbol{k})c_{\unicode[STIX]{x1D6FC}}^{n}(\boldsymbol{k}^{\prime })\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\boldsymbol{k}-\boldsymbol{k}^{\prime },\text{i}\boldsymbol{k})b_{\unicode[STIX]{x1D6FD}}^{i}(\boldsymbol{k})+a_{n}(\boldsymbol{k}^{\prime })(c_{\unicode[STIX]{x1D6FC}}^{i}(\boldsymbol{k})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\boldsymbol{k}^{\prime }-\boldsymbol{k},\text{i}\boldsymbol{k}^{\prime })b_{\unicode[STIX]{x1D6FD}}^{n}(\boldsymbol{k}^{\prime }))^{\ast }}{(\unicode[STIX]{x1D714}_{i}(\boldsymbol{k})-\unicode[STIX]{x1D714}_{n}(\boldsymbol{k}^{\prime }))^{2}+\unicode[STIX]{x1D703}^{2}}\,\text{d}\boldsymbol{k}^{\prime },\end{eqnarray}$$
where we have let
$\boldsymbol{k}^{\prime }:=\boldsymbol{k}+\boldsymbol{p}$
. Setting the regularisation parameter
$\unicode[STIX]{x1D703}\rightarrow 0$
, we have that
$\unicode[STIX]{x1D703}/(x^{2}+\unicode[STIX]{x1D703}^{2})\rightarrow \unicode[STIX]{x03C0}\unicode[STIX]{x1D6FF}(x)$
. This leads to a factor
$\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D714}_{i}(\boldsymbol{k})-\unicode[STIX]{x1D714}_{n}(\boldsymbol{k}^{\prime }))$
which indicates that scattering is restricted within a single branch of the dispersion relation, and so we may drop the sum over
$n$
.
In order to evaluate this expression, we define
$$\begin{eqnarray}c_{\unicode[STIX]{x1D706}}(\boldsymbol{k})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D707}}(\boldsymbol{k}^{\prime }-\boldsymbol{k},\text{i}\boldsymbol{k}^{\prime })b_{\unicode[STIX]{x1D707}}(\boldsymbol{k}^{\prime }):=\unicode[STIX]{x1D6FC}(\boldsymbol{k},\boldsymbol{k}^{\prime })+\text{i}\unicode[STIX]{x1D6FD}(\boldsymbol{k},\boldsymbol{k}^{\prime }).\end{eqnarray}$$
We find that
$$\begin{eqnarray}\unicode[STIX]{x1D6FC}(\boldsymbol{k},\boldsymbol{k}^{\prime })=\frac{\boldsymbol{k}^{\prime }\times \boldsymbol{k}}{\unicode[STIX]{x1D714}^{2}|\boldsymbol{k}|^{2}}[(f^{2}+\unicode[STIX]{x1D714}^{2})\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}^{\prime }-f^{2}|\boldsymbol{k}|^{2}],\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D6FD}(\boldsymbol{k},\boldsymbol{k}^{\prime })=\frac{f\unicode[STIX]{x1D714}}{\unicode[STIX]{x1D714}^{2}|\boldsymbol{k}|^{2}}[|\boldsymbol{k}\times \boldsymbol{k}^{\prime }|^{2}+\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}^{\prime }(|\boldsymbol{k}|^{2}-\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{k}^{\prime })],\end{eqnarray}$$
and we see that
$\unicode[STIX]{x1D6FC}(\boldsymbol{k},\boldsymbol{k}^{\prime })=-\unicode[STIX]{x1D6FC}(\boldsymbol{k}^{\prime },\boldsymbol{k})$
and
$\unicode[STIX]{x1D6FD}(\boldsymbol{k},\boldsymbol{k}^{\prime })=\unicode[STIX]{x1D6FD}(\boldsymbol{k}^{\prime },\boldsymbol{k})$
. Using these symmetries makes it straightforward to show that
$$\begin{eqnarray}\displaystyle & & \displaystyle \text{Re}\,(c_{\unicode[STIX]{x1D706}}(\boldsymbol{k})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D707}}(\boldsymbol{k}^{\prime }-\boldsymbol{k},\text{i}\boldsymbol{k}^{\prime })b_{\unicode[STIX]{x1D707}}(\boldsymbol{k}^{\prime }))(c_{\unicode[STIX]{x1D6FC}}(\boldsymbol{k}^{\prime })\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\boldsymbol{k}-\boldsymbol{k}^{\prime },\text{i}\boldsymbol{k})b_{\unicode[STIX]{x1D6FD}}(\boldsymbol{k}))\nonumber\\ \displaystyle & & \displaystyle \quad =-(\unicode[STIX]{x1D6FC}^{2}(\boldsymbol{k},\boldsymbol{k}^{\prime })+\unicode[STIX]{x1D6FD}^{2}(\boldsymbol{k},\boldsymbol{k}^{\prime })),\end{eqnarray}$$
and
$$\begin{eqnarray}\text{Re}\,(c_{\unicode[STIX]{x1D706}}(\boldsymbol{k})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D707}}(\boldsymbol{k}^{\prime }-\boldsymbol{k},\text{i}\boldsymbol{k}^{\prime })b_{\unicode[STIX]{x1D707}}(\boldsymbol{k}^{\prime }))(c_{\unicode[STIX]{x1D6FC}}(\boldsymbol{k})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}(\boldsymbol{k}^{\prime }-\boldsymbol{k},\text{i}\boldsymbol{k}^{\prime })b_{\unicode[STIX]{x1D6FD}}(\boldsymbol{k}^{\prime }))^{\ast }=\unicode[STIX]{x1D6FC}^{2}(\boldsymbol{k},\boldsymbol{k}^{\prime })+\unicode[STIX]{x1D6FD}^{2}(\boldsymbol{k},\boldsymbol{k}^{\prime }).\end{eqnarray}$$
We may thus define the differential scattering cross-section as
$$\begin{eqnarray}\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime }):=2\unicode[STIX]{x03C0}|c_{\unicode[STIX]{x1D706}}(\boldsymbol{k})\widehat{\unicode[STIX]{x1D61C}}_{\unicode[STIX]{x1D706}\unicode[STIX]{x1D707}}(\boldsymbol{k}^{\prime }-\boldsymbol{k},\text{i}\boldsymbol{k})b_{\unicode[STIX]{x1D707}}(\boldsymbol{k}^{\prime })|^{2}\widehat{R}(\boldsymbol{k}^{\prime }-\boldsymbol{k})\,\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D714}(\boldsymbol{k})-\unicode[STIX]{x1D714}(\boldsymbol{k}^{\prime }))\end{eqnarray}$$
such that the amplitudes satisfy the kinetic equation
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}a+\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}a=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })a(\boldsymbol{k}^{\prime })\,\text{d}\boldsymbol{k}^{\prime }-\unicode[STIX]{x1D6F4}(\boldsymbol{k})a(\boldsymbol{k}),\end{eqnarray}$$
with
$$\begin{eqnarray}\unicode[STIX]{x1D6F4}(\boldsymbol{k})=\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })\,\text{d}\boldsymbol{k}^{\prime }.\end{eqnarray}$$
Expanding the terms in (B 52), replacing the power spectrum with the energy spectrum using (B 40), and expressing the delta function instead in terms of wavevectors, we find that the differential scattering cross-section for the shallow-water system (2.2) is given by (2.9).
The conservation of the leading-order energy is established by noting that the equivalent of (B 14) for the scaled Wigner function is
$$\begin{eqnarray}\displaystyle {\mathcal{E}}(\boldsymbol{x},t) & = & \displaystyle \frac{1}{2}\text{tr}\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D648}[\unicode[STIX]{x1D652}^{(0)}+\unicode[STIX]{x1D700}^{1/2}\unicode[STIX]{x1D652}^{(1)}+O(\unicode[STIX]{x1D700})]\,\text{d}\boldsymbol{k}\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2}\int _{\mathbb{R}^{2}}\mathop{\sum }_{j=\pm }a_{j}(\boldsymbol{x},\boldsymbol{k},t)\,\text{d}\boldsymbol{k}+O(\unicode[STIX]{x1D700}^{1/2})\nonumber\\ \displaystyle & = & \displaystyle {\mathcal{E}}_{0}(\boldsymbol{x},t)+O(\unicode[STIX]{x1D700}^{1/2}).\end{eqnarray}$$
Integrating (B 53) with respect to
$\boldsymbol{k}$
and noting that the right-hand side vanishes because of the symmetry
$\unicode[STIX]{x1D70E}(\boldsymbol{k},\boldsymbol{k}^{\prime })=\unicode[STIX]{x1D70E}(\boldsymbol{k}^{\prime },\boldsymbol{k})$
, we find leading-order energy conservation
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}{\mathcal{E}}_{0}+\unicode[STIX]{x1D735}_{\!\!\boldsymbol{x}}\boldsymbol{\cdot }{\mathcal{F}}_{0}=0,\end{eqnarray}$$
with the leading-order energy flux
$$\begin{eqnarray}{\mathcal{F}}_{0}(\boldsymbol{x},t):=\int _{\mathbb{R}^{2}}[\unicode[STIX]{x1D735}_{\!\!\boldsymbol{k}}\unicode[STIX]{x1D714}(\boldsymbol{k})]a(\boldsymbol{x},\boldsymbol{k},t)\,\text{d}\boldsymbol{k}.\end{eqnarray}$$
Appendix C. Projection of simulation data onto modes
We describe how the energy density
$a(\boldsymbol{x},\boldsymbol{k},t)$
can be estimated from local Fourier transforms of the wave fields. Using the Fourier representation of the Wigner function given in (B 15) (with
$\unicode[STIX]{x1D700}=1$
), it is easily verified that the projection property
$$\begin{eqnarray}\int _{\mathbb{R}^{2}}\unicode[STIX]{x1D652}(\boldsymbol{x},-\boldsymbol{k},t)\,\text{d}\boldsymbol{x}=|\widehat{\unicode[STIX]{x1D753}}(\boldsymbol{k})|^{2}\end{eqnarray}$$
holds. In order to discriminate between the energy contributions from the different modes, we expand the fields in the eigenvectors basis according to
$$\begin{eqnarray}\widehat{\unicode[STIX]{x1D753}}(\boldsymbol{k},t)=\mathop{\sum }_{j=\pm }A^{(j)}(\boldsymbol{k},t)\boldsymbol{b}_{j}(\boldsymbol{k}),\end{eqnarray}$$
so that the energy for each wavenumber is given by
$$\begin{eqnarray}\tilde{{\mathcal{E}}}(\boldsymbol{k},t)={\textstyle \frac{1}{2}}\langle \,\widehat{\unicode[STIX]{x1D753}},\widehat{\unicode[STIX]{x1D753}}\rangle _{\unicode[STIX]{x1D648}}={\textstyle \frac{1}{2}}(|A^{(+)}(\boldsymbol{k},t)|^{2}+|A^{(-)}(\boldsymbol{k},t)|^{2}),\end{eqnarray}$$
with
$\langle \,\!\cdot \,,\,\cdot \rangle _{\unicode[STIX]{x1D648}}$
as defined in (B 6). Orthonomality of the eigenvectors means that we can extract the modal energy contributions by projection to find
$$\begin{eqnarray}|A^{(j)}(\boldsymbol{k},t)|^{2}=|\langle \,\boldsymbol{b}_{j}(\boldsymbol{k}),\widehat{\unicode[STIX]{x1D753}}(\boldsymbol{k},t)\rangle _{\unicode[STIX]{x1D648}}|^{2}.\end{eqnarray}$$
In terms of the energy density, the leading-order energy is given by
$$\begin{eqnarray}\tilde{{\mathcal{E}}}_{0}(\boldsymbol{k},t)=\frac{1}{2}\int _{\mathbb{R}^{2}}\mathop{\sum }_{j=\pm }a_{j}(\boldsymbol{x},-\boldsymbol{k},t)\,\text{d}\boldsymbol{x}.\end{eqnarray}$$
Thus, we may track the energy from the ‘
$+$
’ mode by projecting the Fourier transform of the wave fields:
$$\begin{eqnarray}\int _{\mathbb{R}^{2}}a_{+}(\boldsymbol{x},-\boldsymbol{k},t)\,\text{d}\boldsymbol{x}=\langle |\langle \,\boldsymbol{b}_{+}(\boldsymbol{k}),\widehat{\unicode[STIX]{x1D753}}(\boldsymbol{k},t)\rangle _{\unicode[STIX]{x1D648}}|^{2}\rangle =\langle |\unicode[STIX]{x1D714}\widehat{\unicode[STIX]{x1D702}}/|\boldsymbol{k}||^{2}\rangle ,\end{eqnarray}$$
with
$\langle \cdot \rangle$
denoting the ensemble average. This relates to the leading-order energy density of a single wave mode to the sea-surface height.
Open access
