1. Introduction
Atmospheric and oceanic inertia-gravity waves (IGWs) propagate in a complex turbulent flow which is in approximately geostrophic and hydrostatic balance. The inhomogeneities of this flow result in the scattering of IGWs which redistributes their energy across wavevector space. This process has long been thought to play a role in the energetics of the atmosphere and ocean and it has been modelled using a range of approximations (see Müller Reference Müller1976, Reference Müller1977; Watson Reference Watson1985; Müller et al. Reference Müller, Holloway, Henyey and Pomphrey1986; Savva, Kafiabad & Vanneste Reference Savva, Kafiabad and Vanneste2021; Young Reference Young2021).
Kafiabad, Savva & Vanneste (Reference Kafiabad, Savva and Vanneste2019, hereafter Reference Kafiabad, Savva and VannesteKSV) used multiscale asymptotics to show that the wave-action of linear IGWs propagating in a steady random geostrophic flow of much larger spatial scale evolves according to the diffusion equation
\begin{equation} \partial_t a + \boldsymbol{c} \boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{x}} a = \boldsymbol{\nabla}_{\boldsymbol{k}} \boldsymbol{\cdot} \left( \boldsymbol{\mathsf{D}} \boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{k}} a \right) + F. \end{equation}
Here, 
$a(\boldsymbol {x}, \boldsymbol {k}, t)$ is the wave-action density in the 
$(\boldsymbol {x}, \boldsymbol {k})$ phase space, 
$\boldsymbol {k}$ is the wavevector, 
$\boldsymbol {c} = \boldsymbol {\nabla }_{\boldsymbol {k}} \omega$ is the intrinsic group velocity of IGWs, and 
$F(\boldsymbol {x}, \boldsymbol {k}, t)$ is a forcing term. The IGW intrinsic frequency
\begin{equation} \omega = (f^2 \cos^2 \theta + N^2 \sin^2 \theta)^{1/2}, \end{equation}
with 
$f < N$ the Coriolis and buoyancy frequencies, depends on the angle 
$\theta$ between 
$\boldsymbol {k}$ and the vertical. The 
$\boldsymbol {k}$-dependent diffusivity tensor 
$\boldsymbol{\mathsf{D}}$ is given in components by
\begin{equation} {\mathsf{D}}_{ij} = k_m k_n \int_0^\infty \langle \partial_{x_i} U_n(\boldsymbol{x}) \partial_{x_j} U_m({\boldsymbol{x} - \boldsymbol{c}s})\rangle \, \text{d} s, \end{equation}
where 
$\langle {\cdot } \rangle$ denotes ensemble average and 
$\boldsymbol {U}$ is the flow velocity field, with prescribed homogeneous statistics. A striking prediction of the diffusion equation (1.1) is that forced IGWs have a stationary spectrum scaling with wavenumber as 
$k^{-2}$, consistent with observed atmospheric mesoscale spectra (Gage & Nastrom Reference Gage and Nastrom1986; Lindborg Reference Lindborg1999) and oceanic submesoscale spectra (Callies & Ferrari Reference Callies and Ferrari2013). This provides support to the interpretation of the dynamics in these ranges as dominated by almost linear IGWs (Dewan Reference Dewan1979; VanZandt Reference VanZandt1982; Bühler, Callies & Ferrari Reference Bühler, Callies and Ferrari2014; Callies, Ferrari & Bühler Reference Callies, Ferrari and Bühler2014; Callies, Bühler & Ferrari Reference Callies, Bühler and Ferrari2016). (The nature of the dynamics and level of nonlinearity in the atmospheric mesoscales is still a subject of debate; see Li & Lindborg (Reference Li and Lindborg2018) and references therein for a contrasting view.)
Crucially, the assumption of time-independent flow implies that the diffusivity tensor satisfies 
$\boldsymbol{\mathsf{D}} \boldsymbol {\cdot } \boldsymbol {c} = 0$, as shown in Reference Kafiabad, Savva and VannesteKSV. Thus, noting 
$\boldsymbol{\mathsf{D}}$ is symmetric, the diffusive flux 
$\boldsymbol{\mathsf{D}} \boldsymbol {\cdot } \boldsymbol {\nabla }_{\boldsymbol {k}} a$ is perpendicular to 
$\boldsymbol {c} = \boldsymbol {\nabla }_{\boldsymbol {k}} \omega$ and hence the diffusion of wave action is restricted to a constant-frequency surface, namely a cone 
$\theta = \mathrm {const}$. This prediction is the direct consequence of the assumed linearity and time independence. Simulations of the nonlinear Boussinesq equations reported by Reference Kafiabad, Savva and VannesteKSV nonetheless indicate that it applies to a good approximation to small-Rossby-number flows, because their time scale is asymptotically larger than the IGW propagation time scale. This is illustrated in figure 1 which shows the result of a forced nonlinear Boussinesq simulation similar to Reference Kafiabad, Savva and VannesteKSV's (see § 3.2 for details): the energy density in wavevector space is confined close to the constant-
$\theta$ cone corresponding to the forcing frequency.
Figure 1. IGW energy spectrum 
$e$ as a function of the horizontal and vertical wavenumbers 
$(k_h,k_z)$ in the forced Boussinesq simulation described in § 3.2. The wavenumbers are scaled by the forcing wavenumbers 
$(k_{h*},k_{z*})$ indicated by the white crosses. The cone corresponding to the forcing frequency is indicated by the solid lines. The energy 
$e$ is rescaled by a characteristic value. The dashed lines indicate the boundary of the region of validity of the diffusion approximation (see Appendix A).
However, Dong, Bühler & Smith (Reference Dong, Bühler and Smith2020) suggest that the slow diffusion of wave action across constant-frequency surfaces that results from slow flow time dependence causes significant transfer of wave action from low to high frequency and demonstrate this for IGWs in rotating shallow water. The relevance of this result to three-dimensional IGWs is unclear. It is therefore an open question whether flow time dependence can radically alter the phenomenology of IGW diffusion by geostrophic turbulence, possibly on time scales much longer than the length of the simulations reported in Reference Kafiabad, Savva and VannesteKSV and in figure 1.
We address this question in this paper by revisiting Reference Kafiabad, Savva and VannesteKSV to account for the slow time dependence of the geostrophic flow. Our starting point is the McComas & Bretherton (Reference McComas and Bretherton1977) diffusivity
\begin{equation} {\mathsf{D}}_{ij} = k_m k_n \int_0^\infty \langle \partial_{x_i} U_n(\boldsymbol{x}, t) \partial_{x_j} U_m({\boldsymbol{x} - \boldsymbol{c}s}, t-s)\rangle \, \mathrm{d} s, \end{equation}
which applies to flows with arbitrary time dependence and was originally derived for wave–wave interactions in the induced diffusion regime. This diffusivity reduces to (1.3) in the time-independent case. (See Dong et al. (Reference Dong, Bühler and Smith2020) for a derivation using multiscale asymptotics.) Under the assumption of slow time dependence, encapsulated by a small parameter 
$\varepsilon$ – the ratio of the geostrophic flow velocity to the IGW group speed – we approximate (1.4) and solve the associated diffusion equation asymptotically to obtain the equilibrium action distribution resulting from a steady single-frequency forcing. The results show that the action remains localised within an 
$O(\varepsilon )$-thick boundary layer around the cone corresponding to the forcing frequency. This indicates that the diffusion of three-dimensional IGWs is largely unaffected by the slow time dependence of geostrophic turbulence. In particular, the 
$k^{-2}$ equilibrium spectrum found by Reference Kafiabad, Savva and VannesteKSV can be recovered by integration of the solution across the boundary layer. We confirm the main theoretical predictions by comparison with a high-resolution simulation of the nonlinear Boussinesq equations as shown in figure 1.
2. Approximation of the diffusivity tensor
In this section we approximate the diffusivity in (1.4) taking advantage of the slow time dependence of the geostrophic flow. Introducing the velocity correlation tensor 
${\varPi}_{mn}(\boldsymbol {y},s) = \langle U_m(\boldsymbol {x} + \boldsymbol {y}, t + s)U_n(\boldsymbol {x}, t)\rangle$ we rewrite (1.4) as
\begin{equation} {\mathsf{D}}_{ij} ={-}\frac{1}{2}k_m k_n \int_{-\infty}^\infty \frac{\partial^2 {\varPi}_{mn}}{\partial y_i \partial y_j} (\boldsymbol{c}s,s) \, \mathrm{d} s, \end{equation}
where we extend the integration range to 
$(-\infty, \infty )$ using that 
$k_m k_n {\varPi}_{mn}(-\boldsymbol {c}s,-s) = k_m k_n {\varPi}_{mn}(\boldsymbol {c}s,s)$. In terms of the wavevector–frequency spectrum 
$\hat {{\varPi }}_{mn}$ defined via the Fourier transform
\begin{equation} {\varPi}_{mn}(\boldsymbol{x}, t) = \int_{\mathbb{R}^4} \hat{{\varPi}}_{mn}(\boldsymbol{K}, {\varOmega}) \exp({\mathrm{i}(\boldsymbol{K}\boldsymbol{\cdot} \boldsymbol{x} -{\varOmega}t)})\, \mathrm{d} \boldsymbol{K} \,\text{d} {\varOmega} \end{equation}this becomes
\begin{equation} {\mathsf{D}}_{ij} = {\rm \pi}k_m k_n \int_{\mathbb{R}^4}K_i K_j \hat{{\varPi}}_{mn}(\boldsymbol{K}, {\varOmega}) \delta(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c}- {\varOmega}) \, \mathrm{d} \boldsymbol{K} \,\text{d} {\varOmega} \end{equation}
on using 
$\int _\mathbb {R} \exp ({\mathrm {i}(\boldsymbol {K}\boldsymbol {\cdot } \boldsymbol {c} -{\varOmega })s}) \, \mathrm {d} s = 2 {\rm \pi}\delta ( \boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c}- {\varOmega })$. Using the spherical polar coordinates 
$(k,\theta,\phi )$ for 
$\boldsymbol {k}$ and 
$(K,{\varTheta },{\varPhi })$ for 
$\boldsymbol {K}$ (lowercase symbols for IGW 
$\boldsymbol {k}$-space and uppercase symbols for geostrophic flow 
$\boldsymbol {K}$-space), we compute
\begin{equation} k_m k_n \hat{{\varPi}}_{mn}(\boldsymbol{K}, {\varOmega}) = (k_1 K_2 - k_2 K_1)^2 E_\psi(\boldsymbol{K},{\varOmega}) = 2 k^2 \sin^2 \theta \sin^2 \gamma E(\boldsymbol{K}, {\varOmega}), \end{equation}
where 
$E_\psi$ is the spectrum of the stream function 
$\psi$ of the geostrophic flow (that is, the Fourier transform of 
$\langle \psi (\boldsymbol {x} + \boldsymbol {y}, t+s) \psi (\boldsymbol {x}, t) \rangle$),
\begin{equation} \gamma = {\varPhi} - \phi, \end{equation}
and 
$E(\boldsymbol {K},{\varOmega })= K^2 \sin ^2 {\varTheta } E_\psi /2$ is the geostrophic flow kinetic-energy spectrum. Substituting (2.4)–(2.5) into (2.3) yields
\begin{equation} {\mathsf{D}}_{ij} = 2 {\rm \pi}k^2 \sin^2 \theta \int_{\mathbb{R}^4}K_i K_j \sin^2 \gamma E(\boldsymbol{K}, {\varOmega}) \delta(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c}- {\varOmega}) \, \mathrm{d} \boldsymbol{K} \,\text{d} {\varOmega}. \end{equation} Following Reference Kafiabad, Savva and VannesteKSV, we assume that the flow is isotropic in the horizontal so that 
$E(\boldsymbol {K}, {\varOmega })$ is independent of 
${\varPhi }$. In spherical polar coordinates, several components of 
$\boldsymbol{\mathsf{D}}$ vanish. To see this, we replace 
${\varPhi }$ by 
$\gamma$ as an integration variable in (2.6) and express 
$\boldsymbol {K}$ in the local spherical basis 
$(\boldsymbol {e}_k, \boldsymbol {e}_\theta, \boldsymbol {e}_\phi )$ associated with 
$\boldsymbol {k}$. Thus, we write
\begin{align} \boldsymbol{K} = K \sin {\varTheta}((\sin \theta \cos \gamma + \cot {\varTheta} \cos \theta) \boldsymbol{e}_k + (\cos \theta \cos \gamma - \cot {\varTheta} \sin \theta)\boldsymbol{e}_\theta +\sin \gamma \boldsymbol{e}_\phi). \end{align}
We can now use the parity of the integrand in (2.6) with respect to 
$\gamma$, noting that 
$\delta (\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c}- {\varOmega })$ is even since 
$\boldsymbol {c} = \boldsymbol {\nabla }_{\boldsymbol {k}} \omega (\theta )$ implies that 
$\boldsymbol {c} \parallel \boldsymbol {e}_\theta$ hence 
$\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c} = \boldsymbol {e}_\theta \boldsymbol {\cdot } \boldsymbol {K} \, c$. The parity of the integrands giving the components 
${\mathsf{D}}_{kk} = \boldsymbol {e}_k \boldsymbol {\cdot } \boldsymbol{\mathsf{D}} \boldsymbol {\cdot } \boldsymbol {e}_k$, etc. of 
$\boldsymbol{\mathsf{D}}$ is then determined by the parity of pairwise products of 
$\boldsymbol {e}_k \boldsymbol {\cdot } \boldsymbol {K}$, 
$\boldsymbol {e}_\theta \boldsymbol {\cdot } \boldsymbol {K}$ and 
$\boldsymbol {e}_\phi \boldsymbol {\cdot } \boldsymbol {K}$. We conclude from this that the only non-zero components of 
$\boldsymbol{\mathsf{D}}$ are 
${\mathsf{D}}_{kk}, {\mathsf{D}}_{k \theta } = {\mathsf{D}}_{\theta k}, {\mathsf{D}}_{\theta \theta }$ and 
${\mathsf{D}}_{\phi \phi }$. Thus, diffusion in the azimuthal direction depends only on azimuthal gradients of action and is decoupled from the 
$k$ and 
$\theta$ directions.
We now restrict our attention to flows that are slowly time dependent in the sense that their typical frequencies 
${\varOmega }$ and wavevectors 
$\boldsymbol {K}$ satisfy 
${\varOmega } \ll \boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c}$. For realistic turbulent flows, 
${\varOmega } \sim UK$, hence this condition is equivalent to the condition 
$U \ll c$ that underpins the diffusion approximation (1.1), the limitation of which is discussed in Appendix A. To make the smallness of 
${\varOmega }$ relative to 
$\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c}$ explicit, we introduce a bookkeeping parameter 
$\varepsilon \ll 1$ to mark out asymptotically small terms. The delta function 
$\delta (\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c} - {\varOmega })$ in (2.6) becomes 
$\delta (\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c} - \varepsilon {\varOmega })$ and can be expanded as
\begin{equation} \delta(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c} - \varepsilon {\varOmega}) = \delta(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c}) - \varepsilon \delta'(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c}) {\varOmega} + \varepsilon^2 \delta''(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c}) {\varOmega}^2/2 + O(\varepsilon^3). \end{equation}
Using this alongside the evenness of 
$E(\boldsymbol {K}, {\varOmega })$ in 
${\varOmega }$ leads to the approximation
\begin{equation} \boldsymbol{\mathsf{D}} = \boldsymbol{\mathsf{D}}^{(0)} + \varepsilon^2 \boldsymbol{\mathsf{D}}^{(1)} +O(\varepsilon^4). \end{equation}Here
\begin{equation} {\mathsf{D}}^{(0)}_{ij} = 2 {\rm \pi}k^2 \sin^2 \theta \int_{\mathbb{R}^3}K_i K_j \sin^2 \gamma E(\boldsymbol{K}) \delta(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c}) \, \mathrm{d} \boldsymbol{K}, \end{equation}
where 
$E(\boldsymbol {K})$ is the geostrophic flow kinetic energy spectrum marginalised over frequencies, recovers the diffusivity of time-independent flows obtained by Reference Kafiabad, Savva and VannesteKSV (up to a factor 
$(2{\rm \pi} )^3$ corresponding to a different Fourier transform convention, see (2.2)). For a horizontally isotropic geostrophic flow, 
${\mathsf{D}}^{(0)}_{ij}$ has two non-zero components in spherical polar coordinates, namely
$$\begin{gather} {\mathsf{D}}^{(0)}_{kk} = \frac{4 {\rm \pi}k^3 \omega \sin^2\theta}{(N^2 - f^2) |\cos^5 \theta|}\int_{-\infty}^{\infty} \int_{\theta}^{{\rm \pi} - \theta} K^3 \cos^2{\varTheta} (\cot^2\theta - \cot^2 {\varTheta})^{1/2}E(K,{\varTheta}) \, \mathrm{d} K \,\text{d} {\varTheta}, \end{gather}$$
$$\begin{gather}{\mathsf{D}}^{(0)}_{\phi\phi} = \frac{4 {\rm \pi}k^3 \omega \sin^4\theta}{(N^2 - f^2) |\cos^5 \theta|}\int_{-\infty}^{\infty} \int_{\theta}^{{\rm \pi} - \theta} K^3 \sin^2{\varTheta} (\cot^2\theta - \cot^2 {\varTheta})^{3/2}E(K,{\varTheta}) \, \mathrm{d} K \,\text{d} {\varTheta}. \end{gather}$$
(These equations are (A 13) in Reference Kafiabad, Savva and VannesteKSV, up to the 
$(2{\rm \pi} )^3$ factor and a typographical correction in the lower limit of 
$\varTheta$.)
The leading-order correction to (2.10) induced by the slow flow time dependence is
\begin{equation} {\mathsf{D}}^{(1)}_{ij} = {\rm \pi}k^2 \sin^2 \theta \int_{\mathbb{R}^3}K_i K_j \sin^2 \gamma A(\boldsymbol{K}) \delta''(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c}) \, \mathrm{d} \boldsymbol{K}, \end{equation}and depends on the geostrophic-flow acceleration spectrum
\begin{equation} A(\boldsymbol{K}) = \int_{\mathbb{R}} E(\boldsymbol{K}, {\varOmega}){\varOmega}^2 \, \mathrm{d} {\varOmega}, \end{equation}a natural measure of the flow's unsteadiness.
It turns out that the only dynamically significant component of 
$\boldsymbol{\mathsf{D}}^{(1)}$ is 
${\mathsf{D}}^{(1)}_{\theta \theta }$, corresponding to across-cone diffusion, on which we now concentrate. Contracting (2.12) twice with 
$\boldsymbol {e}_\theta = \boldsymbol {c}/c$, we obtain
\begin{equation} {\mathsf{D}}^{(1)}_{\theta \theta}= \frac{\boldsymbol{c} \boldsymbol{\cdot} \boldsymbol{\mathsf{D}}^{(1)} \boldsymbol{\cdot} \boldsymbol{c}}{c^2} = \frac{{\rm \pi} k^2 \sin^2 \theta}{c^2} \int_{\mathbb{R}^3} (\boldsymbol{K}\boldsymbol{\cdot} \boldsymbol{c})^2 \sin^2 \gamma A(\boldsymbol{K}) \delta''(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c})\, \mathrm{d} \boldsymbol{K}. \end{equation}Noting that
\begin{equation} \int_\mathbb{R} x^2f(x)\delta''(x)\, \mathrm{d}\kern0.7pt x = 2 \int_\mathbb{R} f(x)\delta(x)\, \mathrm{d}\kern0.7pt x \end{equation}
for any smooth 
$f(x)$ reduces (2.14) to
\begin{equation} {\mathsf{D}}^{(1)}_{\theta \theta} = \frac{2 {\rm \pi}k^2 \sin^2 \theta}{c^2} \int_{\mathbb{R}^3} \sin^2 \gamma A(\boldsymbol{K}) \delta(\boldsymbol{K} \boldsymbol{\cdot} \boldsymbol{c})\, \mathrm{d} \boldsymbol{K}. \end{equation}
Representing 
$\boldsymbol {K}$ in the polar spherical coordinates 
$(K,{\varTheta },{\gamma })$ and expanding 
$\boldsymbol {K} \boldsymbol {\cdot } \boldsymbol {c}$ using (2.7) gives
$$\begin{gather} {\mathsf{D}}^{(1)}_{\theta \theta} = \frac{2 {\rm \pi}k^2 \sin^2 \theta}{c^2} \int_{0}^{\infty} \,\mathrm{d} K\int_{0}^{\rm \pi} \,\text{d} {\varTheta} \int_{-{\rm \pi}}^{\rm \pi} \,\text{d} \gamma K^2 \sin {\varTheta} \sin^2 \gamma A(K,{\varTheta}) \nonumber\\ \times \delta\left(K c \sin {\varTheta} \cos\theta \left(\cos \gamma - \frac{\cot {\varTheta}}{\cot \theta}\right)\right), \end{gather}$$
where we use horizontal isotropy to write 
$A(\boldsymbol {K})= A(K,{\varTheta })$. Under the change of variable 
$\zeta = \cos \gamma$ this simplifies into
\begin{align} {\mathsf{D}}^{(1)}_{\theta \theta} = \frac{4 {\rm \pi}k^2 \sin^2 \theta}{c^3 |\cos \theta|} \int_{0}^{\infty}\,\text{d} K\int_{0}^{\rm \pi} \,\text{d} {\varTheta} \int_{- 1}^{1} \,\text{d} \zeta \, K (1-\zeta^2)^{1/2} A(K,{\varTheta}) \delta\left(\zeta - \frac{\cot {\varTheta}}{\cot \theta}\right), \end{align}
where the factor of 2 arises from the evenness of 
$\cos \gamma$. Only values of 
${\varTheta }$ for which 
$|\cot {\varTheta }/\cot \theta | < 1$ contribute to the integral, which reduces the integration range to 
$(\theta, {\rm \pi}- \theta )$. Integrating over 
$\zeta$ then yields
\begin{equation} {\mathsf{D}}^{(1)}_{\theta \theta} = \frac{4 {\rm \pi}k^2 \sin^2 \theta}{c^3 |\cos \theta|} \int_{0}^{\infty}\,\text{d} K\int_{\theta}^{{\rm \pi} - \theta} \,\text{d} {\varTheta} \, K \left(1-\left(\frac{\cot {\varTheta}}{\cot \theta}\right)^2\right)^{1/2} A(K, {\varTheta}). \end{equation}Substituting in
\begin{equation} c = |\boldsymbol{\nabla}_{\boldsymbol{k}} \omega| = \frac{1}{k}\partial_\theta (f^2 \cos^2 \theta + N^2 \sin^2 \theta)^{1/2} = \frac{(N^2 - f^2)|\sin\theta \cos\theta|}{\omega k}, \end{equation}and rearranging gives the final form
\begin{equation} {\mathsf{D}}^{(1)}_{\theta \theta} = \frac{4 {\rm \pi}\omega^3 k^5}{(N^2 - f^2)^3 |\cos^5 \theta| } \int_{0}^{\infty}\int_{\theta}^{{\rm \pi} - \theta} K (\cot^2\theta- \cot^2 {\varTheta})^{1/2} A(K,{\varTheta})\, \text{d} K\,\text{d} {\varTheta}. \end{equation} In summary, the diffusivity with time-dependent geostrophic flow has three significant components: 
${\mathsf{D}}_{kk} = {\mathsf{D}}^{(0)}_{kk}$ and 
${\mathsf{D}}_{\phi \phi } = {\mathsf{D}}^{(0)}_{\phi \phi }$ given by (2.11) and dependent on the energy spectrum of the geostrophic flow, and 
${\mathsf{D}}_{\theta \theta } = \varepsilon ^2 {\mathsf{D}}^{(1)}_{\theta \theta }$ given by (2.21) and dependent on the flow acceleration spectrum. The small, non-zero 
${\mathsf{D}}_{\theta \theta }$ for non-vanishing flow acceleration captures the weak across-cone diffusion pointed out by Dong et al. (Reference Dong, Bühler and Smith2020).
3. Equilibrium spectrum
3.1. Solution of the steady diffusion equation
We now focus on the response to the spatially homogeneous, azimuthally isotropic steady forcing
\begin{equation} F(\boldsymbol{k}) = \delta(k-k_*) \delta(\theta - \theta_*) \end{equation}
corresponding to a single IGW frequency. (The response to a forcing with arbitrary dependence on 
$k$ and 
$\theta$ can be obtained by integration.) We aim to show that the action density reaches an equilibrium 
$a(k,\theta )$ that is localised near 
$\theta = \theta _*$ – in other words, that the frequencies remain close to the forcing frequency for all time. This is in contrast with the two-dimensional case of Dong et al. (Reference Dong, Bühler and Smith2020) for which no such localised equilibrium exists.
For ease of interpretation, we replace the action density by the energy density 
$e(k,\theta ) = 2{\rm \pi} k^2 \sin \theta \omega a(k, \theta )$, such that 
$e \, \text {d}k \,\text {d}\theta$ is the energy contained in the box 
$[k, k + \text {d}k]$ and 
$[\theta, \theta + \text {d} \theta ]$. Equation (1.1) then reduces to
$$\begin{gather} \partial_k\left( k^2 \left({\mathsf{D}}^{(0)}_{kk} + \varepsilon^2 {\mathsf{D}}^{(1)}_{kk} \right) \partial_k \frac{e}{k^2} + \varepsilon^2 \frac{\sin\theta \omega}{k} {\mathsf{D}}^{(1)}_{k \theta} \partial_\theta \frac{e}{\sin\theta \omega} \right)\nonumber\\ + \varepsilon^2 \omega k \partial_\theta \left(\frac{1}{\omega} {\mathsf{D}}^{(1)}_{k \theta}\partial_k \frac{e}{k^2} + \frac{\sin \theta}{k^3} {\mathsf{D}}^{(1)}_{\theta \theta} \partial_\theta \frac{e}{\sin \theta \omega} \right) ={-}\delta(k-k_*) \delta(\theta - \theta_*), \end{gather}$$
where we ignore unimportant prefactors on the right-hand side. We seek solutions localised in 
$\theta$ in a boundary layer of thickness 
$\varepsilon$ around 
$\theta _*$, assuming
\begin{equation} \sigma = (\theta-\theta_*)/\varepsilon = O(1). \end{equation}
To leading-order in 
$\varepsilon$, (3.2) reduces to
\begin{equation} \partial_k\left( k^2 {\mathsf{D}}^{(0)}_{kk}(\theta_*) \partial_k \frac{e}{k^2}\right) + \frac{1}{k^2} {\mathsf{D}}^{(1)}_{\theta \theta}(\theta_*) \partial_{\sigma \sigma} e ={-}\delta(k - k_*)\delta(\sigma), \end{equation}
ignoring again a prefactor on the right-hand side (in this case 
$1/\varepsilon$). Note that 
${\mathsf{D}}^{(1)}_{\theta \theta }$ is the only correction to the diffusivity tensor induced by flow time dependence that appears in (3.4). (This also applies to anisotropic IGWs in the sense that 
$\partial _\phi a \not = 0$.) This correction appears at leading order, even though the corresponding diffusivity 
$\varepsilon ^2 {\mathsf{D}}^{(1)}_{\theta \theta }$ is small, because of the large gradients in 
$\theta$ of the solution.
We make the dependence on 
$k$ of the diffusivity components 
${\mathsf{D}}^{(0)}_{kk}$ and 
${\mathsf{D}}^{(1)}_{\theta \theta }$ in (2.11a) and (2.21) explicit by writing
\begin{equation} {\mathsf{D}}^{(0)}_{kk} = Q(\theta) k^3 \quad \text{and} \quad {\mathsf{D}}^{(1)}_{\theta \theta} = R(\theta) k^5. \end{equation}Under the change of variables
\begin{equation} e = \bar{e}/(Q_*R_*)^{1/2} \quad \text{and} \quad \sigma = \bar{\sigma} (R_{*}/Q_*)^{1/2}, \end{equation}
where 
$Q_* = Q(\theta _*)$ and 
$R_* = R(\theta _*)$, (3.4) becomes
\begin{equation} k^3 \partial_{kk}\bar{e} + k^2 \partial_k \bar{e} - 4 k \bar{e} + k^3 \partial_{\bar{\sigma}\bar{\sigma}} \bar{e} ={-} \delta(k - k_*)\delta(\bar{\sigma}). \end{equation}In the following, we drop the overbars for simplicity.
We now solve the rescaled problem (3.7). Taking a Fourier transform in 
$\sigma$, with 
$l$ the corresponding Fourier variable, we find
\begin{equation} k^3 \partial_{kk}\hat{e} + k^2\partial_{k}\hat{e} -4k\hat{e} - k^3 l^2 \hat{e} ={-}\frac{\delta(k-k_*)}{2{\rm \pi}}, \end{equation}where the hat denotes the Fourier transform. The solution to the homogeneous problem can be written in terms of modified Bessel functions (DLMF 2022, Chapter 10), leading to the piecewise expression
\begin{equation} \hat{e}(k,l) =\left\{ \begin{array}{@{}ll} A(l)\mathrm{I}_2(|l|k) + B(l)\mathrm{K}_2(|l|k), & \text{for} \ 0 < k < k_*, \\ C(l)\mathrm{I}_2(|l|k) + D(l)\mathrm{K}_2(|l|k), & \text{for} \ k > k_*, \end{array} \right. \end{equation}
where 
$\mathrm {I}$ and 
$\mathrm {K}$ are modified Bessel functions of the first and second kind, and 
$A, B, C$ and 
$D$ are so far arbitrary functions of 
$l$. These functions are determined by the boundary and jump conditions. Finiteness as 
$k \to 0$ and 
$k \to \infty$ requires that 
$B(l)=C(l)=0$. Imposing continuity at 
$k_*$ and the jump 
$[\partial _k \hat e]_{k_*^-}^{k_*^+} = -1/(2 {\rm \pi}k_*^3)$ then gives
\begin{equation} \begin{pmatrix} A\\ D \end{pmatrix} = \frac{1}{2 {\rm \pi}\mathcal{W}\{\mathrm{K}_2(|l|k_*), \mathrm{I}_2(|l|k_*)\}|l|k_*^3} \begin{pmatrix} \mathrm{K}_2(|l|k_*)\\ \mathrm{I}_2(|l|k_*) \end{pmatrix} = \frac{1}{2{\rm \pi} k_*^2} \begin{pmatrix} \mathrm{K}_2(|l|k_*)\\ \mathrm{I}_2(|l|k_*) \end{pmatrix}, \end{equation}
where 
$\mathcal {W}$ is the Wronskian and we use that 
$\mathcal {W}\{\mathrm {K}_2(z), \mathrm {I}_2(z)\} = {1}/{z}$ (DLMF 2022, Eq. (10.28.2)). Hence the solution in Fourier space is
\begin{equation} \hat{e}(k,l) =\frac{1}{2 {\rm \pi}k_*^2} \times \left\{ \begin{array}{@{}ll} \mathrm{K}_2(|l|k_*)\mathrm{I}_2(|l|k), & \text{for} \ 0 < k < k_*, \\ \mathrm{I}_2(|l|k_*)\mathrm{K}_2(|l|k), & \text{for} \ k > k_* . \end{array} \right. \end{equation} We invert the Fourier transform. As 
$\hat {e}$ is symmetric in 
$l$, the inverse of (3.11) is
\begin{equation} e(k, \sigma) = \frac{1}{{\rm \pi} k_*^2} \times \left\{ \begin{array}{@{}ll} \int_{0}^{\infty}\mathrm{K}_2(lk_*)\mathrm{I}_2(lk)\cos(\sigma l) \, \mathrm{d} l, & \text{for} \ 0 < k < k_*,\\ \int_{0}^{\infty}\mathrm{I}_2(lk_*)\mathrm{K}_2(lk)\cos(\sigma l) \, \mathrm{d} l, & \text{for} \ k > k_*. \end{array} \right. \end{equation}This can be evaluated exactly using Equation (4), § 6.672 of Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2014),
\begin{equation} \int_0^\infty \mathrm{K}_{v}(ax) \mathrm{I}_{v}(bx) \cos(cx) \, \mathrm{d}\kern0.7pt x = \frac{1}{2(ab)^{1/2}}\mathrm{Q}_{v - 1/2} \left(\frac{a^2 + b^2 + c^2}{2ab} \right), \end{equation}
which holds provided that 
$\mathrm {Re}(a) > |\mathrm {Re}(b)|$ and 
$\mathrm {Re}(v)> -1/2$. Here, 
$\mathrm {Q}_{v - 1/2}$ is the Legendre function of the second kind (DLMF 2022, Chapter 14). Clearly, the condition on 
$v$ is satisfied for (3.12). For 
$0< k< k_*$, 
$a = k_* > k = b$; for 
$k > k_*$, 
$a = k > k_* = b$. Therefore, the condition on 
$a$ and 
$b$ also holds. Due to the symmetry of the solution under exchanges of 
$a$ and 
$b$, both integrals in (3.12) are equivalent, leading to
\begin{equation} e(k, \sigma) = \frac{1}{2{\rm \pi} k_{*}^{5/2}k^{1/2}}\mathrm{Q}_{3/2} \left(\frac{k_*^{2} + k^{2} + \sigma^2}{2k_* k} \right). \end{equation} Equation (3.14) is the main result of the paper. It gives the form of the equilibrium distribution of IGW energy forced at a single wavenumber and frequency, accounting for the time dependence of the turbulence. Since 
$\mathrm {Q}_{3/2}$ decays rapidly as its argument increases, (3.14) shows that the IGW energy is localised within an 
$O(\varepsilon )$ layer around the constant frequency cone 
$\theta = \theta _*$ (recall (3.3)). Note that 
$e(k,\sigma )$ has a mild, logarithmic singularity as 
$\sigma \to 0$ for 
$k=k_*$.
We illustrate the form of the energy spectrum predicted by (3.14) in figure 2. Here, 
$e$ scaled by 
$k_*^3$ is plotted against the scaled angle 
$\sigma /k_*$ for a few values of non-dimensionalised total wavenumber 
$k/k_*$. In figure 3, 
$e$ is shown as a function of horizontal and vertical wavenumber and is scaled to approximately match the energy level of the simulation in § 3.2. The value of 
$R_*/Q_*$ is also required for figure 3 and is chosen to match simulation results.
Figure 2. IGW energy spectrum 
$e$ in (3.14) scaled by 
$k_*^3$ as a function of the scaled angle 
$\sigma /k_*$ for a few values of non-dimensionalised total wavenumber 
$k/k_*$.
Figure 3. IGW energy spectrum 
$e$ in (3.14) as a function of horizontal and vertical wavenumbers 
$(k_h, k_z)$. The wavenumbers are scaled by the forcing wavenumbers 
$(k_{h*}, k_{z*})$ indicated by the white crosses. The cone corresponding to the forcing frequency is indicated by the solid lines. The parameters 
$Q_*$ and 
$R_*$ are chosen to match the simulation results in § 3.2 (cf. figure 1).
A useful approximation to (3.14) is obtained from the asymptotics of the Legendre function for large argument,
\begin{equation} e(k, \sigma) \sim \frac{3}{16} k^{{-}3}\left(1 + \frac{k_{*}^{2} + \sigma^2}{k^2}\right)^{{-}5/2}, \end{equation}
which applies for 
$k \to 0$, 
$k \to \infty$ or 
$\sigma \to \infty$. In particular, it makes it possible to characterise the angular localisation of the energy by the power law
\begin{equation} e(k, \sigma) \sim \frac{3}{16} \frac{k^2}{\sigma^5} \quad \text{as} \ \sigma \rightarrow \infty. \end{equation}
Equation (3.15) further shows that at fixed 
$\sigma$, that is, at fixed angle 
$\theta$ or frequency, 
$e(k,\sigma ) \propto k^{-3}$ as 
$k \to \infty$, and 
$e(k,\sigma ) \propto k^{2}$ for 
$k \ll k_*$.
Another limit of interest deduced from (3.15) is
\begin{equation} e(k, \sigma) \sim \frac{3}{16} k^{{-}3}\left(1 + \left(\frac{\sigma}{k}\right)^2\right)^{{-}5/2} \quad \text{as} \ k \to \infty, \ \sigma/k=O(1), \end{equation}
which shows that the spectrum broadens in 
$\sigma$ like 
$k$. Consequently, integration of (3.15) across angles results in a spectrum decaying like 
$k^{-2}$. In fact, the integrated spectrum is exactly proportional to 
$k^{-2}$ for 
$k > k_*$: indeed, integration of (3.7) with respect to 
$\bar {\sigma }$ recovers the equation found by Reference Kafiabad, Savva and VannesteKSV for time-independent flows, with solution proportional to 
$k^{-2}$ for 
$k > k_*$ and 
$k^2$ for 
$k< k_*$.
In dimensional terms, the thickness of the boundary layer around the cone is proportional to the square root of the ratio 
$R_*/Q_*$ (see (3.6b)), which roughly amounts to the ratio of the flow acceleration variance to its energy, and can be interpreted as the relevant flow frequency. This increases when the flow becomes more transient resulting in a thicker boundary layer.
3.2. Comparison with Boussinesq simulations
We compare the analytical prediction (3.14) with the results of a high-resolution three-dimensional Boussinesq simulation. We solve the non-hydrostatic Boussinesq equations using a dealiased pseudospectral code adopted from that in Waite & Bartello (Reference Waite and Bartello2006). A third-order Adams–Bashforth scheme with timestep 
$0.0044/f$ is employed for time integration. The triply periodic domain 
$[0, 2 {\rm \pi}]^2 \times [0, 2 {\rm \pi}f/N]$ is discretised with 
$2304^3$ grid points. A hyperdissipation of the form 
$\nu _h (\partial _x^2+\partial _y^2)^4 + \nu _z \partial _{z}^8$, with 
$\nu _h = 7.8 \times 10^{-23}$ and 
$\nu _z = 7.1 \times 10^{-35}$ (in dimensionless units, with the domain size as reference length and 
$f^{-1}$ as reference time), is implemented in the momentum and buoyancy equations. We take 
$N/f=32$, a representative value of mid-depth ocean stratification. We initialise the simulation with a fully developed geostrophic turbulent flow, which is the output of a decaying quasigeostrophic model with the initial energy spectrum proportional to 
$\exp {( -(((K_h^2 + f^2 K_z^2/N^2)^{1/2}-24)/10)^2)}$. This model is run until the energy spectrum fills the spectral space, peaking at 
$K_h = 4$ and scaling approximately as 
$K_h^{-3}$ and 
$K_z^{-3}$. The flow parameters are selected such that the Rossby number based on the vertical vorticity 
$\zeta$ is 
${Ro} = \langle \zeta ^2 \rangle ^{1/2}/f = 0.11$. Throughout the simulation, an Ornstein–Uhlenbeck forcing with short correlation time (three timesteps) is applied to the linear wave modes with 
$(k_{h*},k_{z*})=(12,221)$ corresponding to the fixed IGW frequency of 
$2f$. This relatively low frequency is chosen so that the aspect ratio of the IGWs is similar to the aspect ratio 
$N/f$ of the geostrophic flow and thus the IGWs are well resolved with the anisotropic grid we use. The simulation is performed until 
$t = 160/f$ by which time the statistics are approximately stationary. We separate IGWs from the mean flow (both for forcing and extracting energy spectra) using the normal-mode decomposition of Bartello (Reference Bartello1995).
We compare the functional form implied by (3.14) with the spectrum 
$e(k_h,k_z)$ obtained in the simulation. This involves fitting two parameters, one that fixes the scale of 
$e$ and corresponds to strength of the forcing, and the other that fixes the scale of 
$\sigma$ and corresponds to 
$(R_*/Q_*)^{1/2}$ (see (3.6b)). We estimate these two parameters by matching the simulation spectrum as a function of 
$\theta - \theta _*$ for 
$k \gtrsim 5 k_*$ as shown in figure 4. These values of 
$k$ are large enough for the perturbation induced by the non-ideal nature of the forcing in the simulation to be negligible, and for discretisation effects to play only a minor role. A difficulty, evident in figure 4, is that the simulation spectrum is not symmetric. We attribute this to an edge effect caused by the proximity of the IGW frequency 
$\omega = 2 f$ to the minimum allowable frequency 
$\omega = f$, and to the breakdown of the diffusion approximation when 
$\omega$ is close to 
$f$ (see Appendix A for details). (The forcing frequency cone has a small opening angle, 
$\theta _* = \tan ^{-1} (k_{h*}/k_{z*}) \approx 3^\circ$, a feature obscured by the anisotropic scaling of the axes in figures 1 and 3.) We therefore carry out the parameter fitting based on the parts of the curves in figure 4 right of their maxima. We further allow for an offset of 
$\theta - \theta _*$, likely the result of the coarse discretisation of the wavevector in the forcing region.
Figure 4. IGW energy spectrum 
$e$ versus 
$\theta -\theta _*$ for several values of 
$k/k_*$: comparison between simulation results (solid lines) and analytical prediction (3.14) (dashed lines). The scalings of 
$e$ and 
$\sigma \propto \theta - \theta _*$ are chosen for the analytical prediction to best match the simulation data (scaling 
$\sigma$ corresponds to estimating 
$(R_*/Q_*)^{1/2}$, see (3.3) and (3.6b)).
The prediction of (3.14) with the two fitted parameters is shown by the dashed curves in figure 4. The agreement with the numerical results is good: (3.14) captures the localisation of 
$e$ and the general form of its decrease with 
$\theta - \theta _*$ at different values of 
$k$. (We emphasise that the same two parameters are used for all the curves.) A complementary view is provided by figure 5 which shows 
$e$ obtained in the simulation as a function of 
$\theta - \theta _*$ (figure 5a) and of 
$k/k_*$ (figure 5b) in log–log coordinates. The power laws 
$\sigma ^{-5}$ (equivalent to 
$(\theta - \theta _*)^{-5}$), 
$k^{-3}$ and 
$k^2$ derived in (3.15)–(3.16) from (3.14) are shown in their range of expected validity. The 
$\sigma ^{-5}$ and 
$k^{-3}$ power laws are consistent with the data albeit over a limited wavenumber range. We regard this as a reasonable match given the difficulties in capturing such rapid decay in a numerical model, and the pollution by the forcing. The 
$k^{2}$ power law is a poorer match. This is to be expected since the spatial scale-separation assumption between IGWs and geostrophic flow that underpins the diffusion equation (1.1) is not satisfied for wavenumbers smaller than the forcing wavenumber. The numerical spectrum for small 
$k$ is also strongly affected by discretisation effects. Note that the abrupt drop in the tail of spectra in figure 5(b) comes from the truncation of data due to storage limitation; the total energy spectrum shows a smooth transition to dissipation range (not shown).
Figure 5. IGW energy spectrum 
$e$ from simulation data in log–log coordinates: (a) as a function of 
$\theta -\theta _*$ for several values of 
$k/k_*$; and (b) as a function of 
$k/k_*$ for several values of 
$\theta -\theta _*$. Predicted power laws are indicated by dashed lines.
Overall, the simulation results compare as well with (3.14) as can be expected given the numerical challenges posed by the finite resolution, non-ideal forcing and an IGW signal that has both low amplitude and decreases rapidly with 
$k$ and 
$\theta - \theta _*$. We note that it is in principle possible to compute the scaling parameter 
$(R_*/Q_*)^{1/2}$ from simulation data using the explicit expressions for 
$R_*$ and 
$Q_*$ deduced from (2.11a), (2.21) and (3.5a,b). Reference Kafiabad, Savva and VannesteKSV evaluate 
$Q_*$ based on the energy spectrum of the geostrophic flow they estimate from simulation data. An analogous evaluation of 
$R_*$ requires the acceleration spectrum of the geostrophic flow. We leave this computation for future work.
4. Discussion
This paper is part of a sequence of works that apply techniques of waves in random media to address the role of the geostrophic flow in shaping the energy distribution of atmospheric and oceanic IGWs (Danioux & Vanneste Reference Danioux and Vanneste2016; Savva & Vanneste Reference Savva and Vanneste2018; Kafiabad et al. Reference Kafiabad, Savva and Vanneste2019; Savva et al. Reference Savva, Kafiabad and Vanneste2021). Their main assumption is that the flow is weak enough to be regarded as a small perturbation to what would otherwise be IGWs propagating in a medium at rest. The perturbation, physically refraction, can be interpreted as arising from resonant triadic interactions involving two IGW modes and a geostrophic (or vortical) mode – these are known as ‘catalytic interactions’ in recognition of the fact that the geostrophic mode is left unaffected (Lelong & Riley Reference Lelong and Riley1991; Bartello Reference Bartello1995). The present paper further assumes that the IGWs have spatial scales much smaller than the flow scales. In this case, the impact of the flow, modelled as a random field, on the IGWs is a diffusion of wave action in wavevector space. (This is the induced diffusion regime considered by McComas & Bretherton (Reference McComas and Bretherton1977) in the context of wave–wave interactions.) Reference Kafiabad, Savva and VannesteKSV examined this process in some detail and showed, in particular, that it leads to IGW characteristics such as a 
$k^{-2}$ stationary spectrum that are consistent with atmospheric and oceanic observations.
To obtain these results, Reference Kafiabad, Savva and VannesteKSV treated the geostrophic flow as time independent, on the grounds that it evolves on a time scale much longer than the IGW periods. With this assumption, the geostrophic mode has a zero frequency. The (resonant) catalytic interactions therefore involve two IGW modes with exactly the same frequency, and wave action exchanges are restricted to a constant-frequency surface in wavevector space. Here, we revisit this assumption by taking the geostrophic flow to be slowly evolving. In this case, the catalytic interactions are between a low frequency geostrophic mode and two IGWs with slightly different frequencies, and action diffuses across the constant-frequency surface. The question is therefore whether this leads to qualitative changes in the statistics of IGWs, for instance by enabling IGW frequencies to diffuse freely and wave action to spread unimpeded across wavevector space (as was recently shown to be the case for two-dimensional waves by Dong et al. (Reference Dong, Bühler and Smith2020)). The answer is no: we show that the stationary spectrum established by forcing single-frequency IGWs is localised within a boundary layer close to the cone of constant frequency associated with the forcing. Thus, even in the infinite-time limit corresponding to this stationary response, the time dependence of the geostrophic flow has only a minor impact on the IGW scattering. Hence, the conclusions of Reference Kafiabad, Savva and VannesteKSV drawn by neglecting the time dependence hold for realistic slowly evolving flows. In particular, scattering by geostrophic flow does not control the frequency distribution of IGWs which, in the absence of other mechanisms, is determined by the forcing or initial conditions. This is only strictly true over a finite range of wavenumbers 
$k$, since the thickness of the boundary layer increases with 
$k$ (see (3.17)). However, at large 
$k$ the hypotheses of weak flow and linear waves also break down (see Reference Kafiabad, Savva and VannesteKSV) and may have a larger impact than the flow time dependence (see Appendix A for a discussion of the restriction on 
$k$ imposed by the weak-flow hypothesis).
It is worth commenting on the sharp difference between the conclusion drawn here for three-dimensional IGWs in a three-dimensional geostrophic flow and that drawn by Dong et al. (Reference Dong, Bühler and Smith2020) in a two-dimensional set up. This difference stems from the very different geometry of the constant-frequency surfaces which are compact in dimension two (circles) and non-compact in dimension three (cones). In the compact case, an initial distribution of action quickly relaxes to become uniform on constant-frequency circles, then slowly spreads across these circles because of the flow time dependence. The flux of action perpendicular to the constant-frequency circles is small, but it allows for the wave frequencies to change without restriction over long time scales. In contrast, for the (non-compact) cones of the three-dimensional case, there is a non-zero action flux along cones, even in the absence of flow time dependence, corresponding to a forward cascade towards small scales. The flux across cones introduced by the slow time dependence of the geostrophic flow acts therefore only as a small perturbation which barely affects the (non-equilibrium) stationary spectrum at finite distances along the cones.
Acknowledgement
We thank O. Bühler and the anonymous referees for their valuable comments.
Funding
M.R.C. was supported by the MAC-MIGS Centre for Doctoral Training under grant EP/S023291/1 of the UK Engineering & Physical Sciences Research Council (EPSRC). H.A.K. and J.V. were supported by EPSRC, grant EP/W007436/1. J.V. was also supported by the UK Natural Environment Research Council, grant NE/W002876/1.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are openly available in the Geophysical Fluid Dynamics collection of Edinburgh DataShare at https://doi.org/10.7488/ds/3490.
Appendix A. Limitation of the diffusion approximation
The diffusion approximation (1.1) on which (3.14) and Reference Kafiabad, Savva and VannesteKSV rely is valid for 
$U \ll c$. Defining the velocity-based Rossby number, 
${Ro} = U K_h/f$ (rather than the vorticity-based definition of § 3.2), and using 
$c$ in (2.20), we can rewrite this condition as
\begin{equation} \frac{k_h}{K_h} \ll \frac{\left((N/f)^2 - 1\right)\sin^2 \theta \cos \theta}{{Ro} \left((N/f)^2 \sin^2 \theta + \cos^2 \theta\right)^{1/2}}, \end{equation}
where 
$k_h = k \sin \theta$ is the horizontal wavenumber and we have taken 
$0 \leq \theta \leq {\rm \pi}/2$ without loss of generality. Figure 6 displays the right-hand side of (A1) against 
$\theta$ for a range of values of 
$N/f$ typical of the ocean and atmosphere. The figure shows that, for realistic, small Rossby numbers (
${Ro} \in [10^{-2},10^{-1}]$), the range of 
$k_h$ over which the diffusion approximation is valid extends to 20–200 times the typical flow wavenumber 
$K_h$ for all IGWs except those with frequencies very close to 
$f$ (
$\theta =0$) and 
$N$ (
$\theta ={\rm \pi} /2$). (A scattering theory tailored to IGWs with frequencies close to 
$f$, that is, near-inertial waves, is developed in Danioux & Vanneste (Reference Danioux and Vanneste2016).)
Figure 6. Upper bound of 
${Ro}\, k_h/K_h$ given by the right-hand side of (A1) as a function of 
$\theta$ for a range of 
$N/f$ values including our simulation value, 
$N/f = 32$.
To determine the range of 
$k_h$ and 
$k_z$ for which condition (A1) is met in our simulation, we recast (A1) in terms of 
$k/k_*$ as used in figures 4 and 5 to obtain
\begin{equation} \frac{k}{k_*} = \frac{k_h \sin \theta_*}{k_{h*} \sin \theta} = \frac{k_z \cos \theta_*}{k_{z*} \cos \theta} \ll \frac{\left((N/f)^2 - 1\right)\sin \theta \cos \theta}{{Ro} \left((N/f)^2 \sin^2 \theta + \cos^2 \theta\right)^{1/2}} \frac{K_h}{k_*}. \end{equation}
The simulation parameters are: 
$N/f = 32$, 
$k_* = 221.3$ and 
$K_h = 4$. The velocity-based Rossby number is estimated to be 
${Ro} = 0.05$. Using these parameters we compute the curve in the 
$(k_h,k_z)$-plane where (A2) is satisfied as an equality and show the result in figure 7. The two lobes labelled C and D indicate the region of validity of the diffusion approximation. The rectangles labelled A (also shown in the inset) and B show the ranges of 
$k_h$ and 
$k_z$ used in figure 1 and resolved in the simulation, respectively. This confirms that the diffusion approximation applies to the typical wavenumbers considered in our analysis. However, because of the rapid change of 
$k_z$ as 
$\theta$ decreases from 
$\theta _*$, the diffusion approximation can be expected to break down around 
$\theta - \theta _* \approx -0.03$ in figure 4. This likely explains the mismatch between theoretical prediction and simulation results to the left of the curves’ maxima in the figure.
Figure 7. Region of validity of the diffusion approximation in the simulation: inequality (A2) holds in regions C and D, where the diffusion approximation applies; rectangle A (also in the inset) reproduces figure 1 in showing the energy density 
$e$; rectangle B shows the ranges of 
$k_h$ and 
$k_z$ resolved in the simulation.
Open access
