2.1. Derivation

We define a barotropic background geostrophic velocity $V(x)$ and stratification $N^2(z) := -(g/\rho _0)(\partial \bar {\rho }/\partial z)$ , where $\bar {\rho }(z)$ is the background density, $\rho _0$ is the Boussinesq reference density and $g$ the gravitational acceleration. Focusing on NIWs with small aspect ratios, we linearise the hydrostatic, inviscid, adiabatic Boussinesq equations about this background state and make the simplifying assumption that the dynamics are independent of $y$ . The linearised equations are

(2.1a) \begin{align} \frac {\partial u}{\partial t} - fv + \frac {1}{\rho _0}\frac {\partial p}{\partial x} = 0, \end{align}

(2.1b) \begin{align} \frac {\partial v}{\partial t} + u\frac {\partial V}{\partial x} + fu = 0, \end{align}

(2.1c) \begin{align} \frac {1}{\rho _0}\frac {\partial p}{\partial z} - b = 0, \end{align}

(2.1d) \begin{align} \frac {\partial b}{\partial t} + w N^2 = 0, \end{align}

(2.1e) \begin{align} \frac {\partial u}{\partial x} + \frac {\partial w}{\partial z} = 0, \end{align}

where $u,v,w$ are the perturbation velocities, $b := g(\bar {\rho } - \rho )/\rho _0$ is the perturbation buoyancy defined from the density $\rho$ , $p$ is the perturbation pressure and $f$ is the Coriolis frequency, taken to be constant under a traditional $f$ -plane approximation.

By defining a streamfunction $\psi$ such that $u = -\partial \psi / \partial z$ , $w = \partial \psi / \partial x$ and systematically eliminating variables, (2.1) can be reduced to a single equation, i.e.

(2.2a) \begin{equation} \left (\frac {\partial ^2 \,}{\partial t^2} + f_{\textit{eff}}^2\right )\frac {\partial ^2 \psi }{\partial z^2} + N^2\frac {\partial ^2 \psi }{\partial x^2} = 0, \end{equation}

where

(2.2b) \begin{equation} f_{\textit{eff}}^2(x) := f\left (f + \frac {\partial V}{\partial x}\right ) \end{equation}

is the square of the effective Coriolis frequency. Note that with variations in $y$ neglected, the along-filament momentum equation is a conservation equation for the absolute momentum $M := fx + V + v$ . Furthermore, $f_{\textit{eff}}^2(x)$ is proportional to the gradient of the background absolute momentum $V + fx$ . Equation (2.2a ) is separable and so we write

(2.3a) \begin{equation} \psi = \sum _m \mathcal{X}_m(x,t)\mathcal{Z}_m(z), \end{equation}

with the vertical modes defined by

(2.3b) \begin{equation} c_m^2\frac {\partial ^2 \,}{\partial z^2}\mathcal{Z}_m = -N^2\mathcal{Z}_m, \end{equation}

where the mode speed $c_m$ is a constant of separation. Taking $N^2$ to be constant and projecting $\psi$ onto sine modes over a domain of depth $L_z$ , $\mathcal{Z}_m = \sin {m\pi z/L_z}$ , thereby satisfying no-penetration boundary conditions, we have

(2.4) \begin{equation} c_m = \frac {NL_z}{m\pi }. \end{equation}

Alternatively, if the vertical structure is described by a plane wave with wavenumber $k_z$ then $c_m = Nk_z^{-1}$ . Before § 5, we are not particularly interested in the structure of the vertical modes, only in the value of $c_m$ , which increases with increasing stratification or vertical scale. Note that in these vertical mode decompositions we have filtered out the barotropic mode.

With the vertical structure set, the equation for the horizontal structure is

(2.5) \begin{equation} \left [\frac {\partial ^2 \,}{\partial t^2} - c_m^2\frac {\partial ^2 \,}{\partial x^2} + f_{\textit{eff}}^2(x)\right ]\mathcal{X} = 0 ,\end{equation}

where we have dropped the subscript $m$ on $\mathcal{X}$ to lighten the notation. This is the Klein–Gordon equation. First derived as a relativistic wave equation, the Klein–Gordon equation can also be interpreted as describing classical waves in an elastic medium where, in this case, the term involving $f_{\textit{eff}}$ represents elasticity. Both interpretations provide insight into the roles of the background flow and vertical structure of the waves in determining the behaviour of the solutions. For waves with frequencies much greater than $f_{\textit{eff}}$ , the elastic term is negligible and the Klein–Gordon equation reduces to the wave equation with wave speed $c_m$ . Furthermore, the characteristics of (2.5) have slope $c_m$ (in $x,t$ space) and so $c_m$ is the rate at which the Klein–Gordon equation propagates information.

We solve the Klein–Gordon equation subject to the initial conditions

(2.6a,b) \begin{align} \mathcal{X} = 1, \quad \frac {\partial \mathcal{X}}{\partial t} = 0. \end{align}

Such uniform initial conditions are a popular theoretical device (e.g. Balmforth et al. Reference Balmforth, Llewellyn and Young1998; Asselin & Young Reference Asselin and Young2020; Asselin et al. Reference Asselin, Thomas, Young and Rainville2020; Kafiabad, Vanneste & Young Reference Kafiabad, Vanneste and Young2021) motivated by the scenario in which a large-scale wind event impulsively excites ageostrophic motions in the upper ocean. Crucially, the spatial structure of the waves develops through interactions with the background flow rather than being prescribed by the initial conditions.

2.2. The YBJ approximation

Young & Ben Jelloul (Reference Young and Jelloul1997; YBJ) model the evolution of NIWs via a multiple time-scale expansion. They formally justify the expansion by assuming that the Rossby number of the background flow and the wave Burger number $Bu_w := c_m^2/f^2L_w^2$ , where $L_w$ is the horizontal scale of the waves, are small. For the barotropic case with a projection onto vertical modes, the leading-order horizontal wave velocities are given by

(2.7) \begin{equation} u_m + \mathrm{i}v_m = -f^2c_m^2\mathcal{A}\mathrm{e}^{-\mathrm{i}ft}, \end{equation}

where the complex amplitude $\mathcal{A}$ varies on a slow time scale (Balmforth et al. Reference Balmforth, Llewellyn and Young1998). Furthermore, the leading-order vertical velocity, buoyancy and pressure may also be expressed in terms of $\mathcal{A}$ . The evolution equation for $\mathcal{A}$ is

(2.8) \begin{equation} \frac {\partial \mathcal{A}}{\partial t} + \frac {\partial (\varPsi , \mathcal{A})}{\partial (x, y)} + \mathrm{i}\frac {1}{2}\zeta \mathcal{A} = \mathrm{i}\frac {c_m^2}{2f}{\nabla} ^2\mathcal{A}, \end{equation}

where $\varPsi$ is the geostrophic steamfunction, $\partial (\boldsymbol{\cdot },\boldsymbol{\cdot })/\partial (x,y)$ is the Jacobian, $\zeta := {\nabla} ^2 \varPsi \equiv \partial V/\partial x - \partial U/\partial y$ is the vertical vorticity of the background flow and ${\nabla} ^2 := \partial ^2/\partial x^2 + \partial ^2/\partial y^2$ is the horizontal Laplacian. Imposing $\partial /\partial y \equiv 0$ , (2.8) reduces to

(2.9) \begin{equation} 2f\mathrm{i}\frac {\partial \mathcal{A}}{\partial t} = -c_m^2\frac {\partial ^2\mathcal{A}}{\partial x^2} + f\frac {\partial V}{\partial x}\mathcal{A}, \end{equation}

which is the time-dependent Schrödinger equation in one spatial dimension.

To reduce the Klein–Gordon equation (2.5) to the time-dependent Schrödinger equation (2.9), we simply let

(2.10a) \begin{equation} \mathcal{X} = \mathcal{A}\mathrm{e}^{-\mathrm{i}ft} + \mathrm{c.c.} \end{equation}

such that

(2.10b) \begin{equation} \frac {\partial ^2 \,}{\partial t^2}\mathcal{X} = \left [\left (-f^2 -2\mathrm{i}f\frac {\partial \,}{\partial t} + \frac {\partial ^2 \,}{\partial t^2}\right )\mathcal{A}\right ]\mathrm{e}^{-\mathrm{i}ft} + \mathrm{c.c.} \end{equation}

and neglect the $\partial ^2 \mathcal{A}/\partial t^2$ term, consistent with the multiple scales approximation underpinning the YBJ expansion. Thus, the validity of the YBJ approximation in this problem hinges on whether or not we are justified in dropping the $\partial ^2 \mathcal{A}/\partial t^2$ term. This approximation has a direct analogue in quantum field theory where it is used to recover the Schrödinger equation from the Klein–Gordon equation in the non-relativistic limit (e.g. Sterman Reference Sterman1993). Throughout the paper we comment on the validity of the YBJ approximation as a function of the parameters describing the different set-ups we consider.

3.1. Problem set-up

To proceed further we must now specify the background flow $V(x)$ . In our first set-up, we consider a background state defined by a velocity scale and a single length scale. That is, a cyclonic filament in an unbounded domain (figure 1 a) with geostrophic velocity gradient,

(3.1) \begin{equation} \frac {\partial V}{\partial x} = \Delta \kern-1pt V \frac {1}{L_f}\mathcal{F}\left (\frac {x}{L_f}\right )\!, \end{equation}

where $L_f$ is the width of the filament, $\Delta\! V \gt 0$ is the change in geostrophic velocity over the filament and $\mathcal{F}(\eta )$ is a positive function of total integral 1. Since the filament is cyclonic, we expect to observe the radiation of NIWs out of the filament analogous to the results of Kafiabad et al. (Reference Kafiabad, Vanneste and Young2021) who considered axisymmetric vortices in the low-Rossby-number limit.

The boundary conditions are $\mathcal{X} \to \cos {ft}$ as $|x| \to \infty$ . Clearly, the solution is not normalisable, i.e. it has infinite energy, but in this aspect the set-up is similar to a scattering problem in quantum mechanics. Furthermore, we always have the option of reframing our problem in terms of $\mathcal{Y} := \mathcal{X} - \cos {ft}$ . This converts the set-up into a forced problem with friendlier boundary conditions, $\mathcal{Y} \to 0$ as $|x| \to \infty$ , and homogeneous initial conditions, $\mathcal{Y} = \partial \mathcal{Y}/\partial t = 0$ .

Since we are considering a radiation problem, a quantity of considerable interest is the time scale, $T$ , over which this radiation occurs. To quantify this, we define

(3.2a) \begin{equation} T_{10} := \min \left \{\ t \ \left | \ \mathcal{X}(0,t)^2 + f^{-2}\frac {\partial \mathcal{X}}{\partial t}(0,t)^2 \lt \frac {1}{10}^2\right.\right \} \end{equation}

and the YBJ equivalent

(3.2b) \begin{equation} T_{10}^{(\textit{YBJ}\kern1pt)} := \min \left \{\ t \ \left | \ 4|\mathcal{A}(0,t)|^2 \lt \frac {1}{10}^2\right.\right \}\!. \end{equation}

Here $T_{10}$ is the time at which this particular norm of the solution at the centre of the filament first drops to one tenth of its initial value. The factor of 4 in $T_{10}^{(\textit{YBJ}\kern1pt)}$ comes from the fact that the initial value of $\mathcal{A}$ is ${1}/{2}$ .

3.2. The parameter space

With background flow (3.1), the problem is defined by four dimensional parameters: the Coriolis frequency $f$ , the mode speed $c_m$ , the filament width $L_f$ and the filament strength $\Delta \kern-1pt V$ . Therefore, after choosing time and length scales to non-dimensionalise the problem, we are left with a 2-D parameter space. Non-dimensionalising using the available time and length scales,

(3.3a,b) \begin{align} \tilde {t} := ft, \quad \tilde {x} := \frac {x}{L_f}, \end{align}

the Klein–Gordon equation (2.5) becomes

(3.4) \begin{equation} \left [\frac {\partial ^2\,}{\partial \tilde {t}^2} - Bu_m\frac {\partial ^2\,}{\partial \tilde {x}^2} + 1 + \textit{Ro}_{\!f}\mathcal{F}(\tilde {x})\right ]\mathcal{X} = 0. \end{equation}

The two parameters appearing here are the filament Burger and Rossby numbers:

(3.5a,b) \begin{align} Bu_m := \frac {c_m^2}{f^2L_f^2}, \quad \textit{Ro}_{\!f} := \frac {\Delta \kern-1pt V}{fL_f}. \end{align}

However, this pair of non-dimensional numbers is not necessarily the best choice to span the parameter space. Indeed, we find that the ratio of the filament Burger and Rossby numbers,

(3.6) \begin{equation} \gamma _m := \frac {Bu_m}{\textit{Ro}_{\!f}} \equiv \frac {c_m^2}{fL_f\Delta \kern-1pt V}, \end{equation}

can be used to distinguish qualitatively different dynamical regimes. A heuristic explanation of why $\gamma _m$ is an important parameter is that the spatial structure of the waves is set by the dispersive, $Bu_m\partial ^2/\partial \tilde {x}^2$ , and refractive, $\textit{Ro}_{\!f}\mathcal{F}(\tilde {x})$ , terms. As a result, $\gamma _m := Bu_m / \textit{Ro}_{\!f}$ determines whether the spatial scale of the waves is short or long compared with the width of the filament. In particular, we find that, for $\gamma _m \ll 1$ , the spatial scale of the waves is $\sqrt {\gamma _m}L_f \ll L_f$ whereas, for $\gamma _m \gg 1$ , the spatial scale of the waves is $\gamma _m L_f \gg L_f$ . When $\gamma _m$ is large, the waves may penetrate across the cyclonic filament in a manner analogous to quantum tunnelling. Therefore, we call $\gamma _m$ the ‘tunnelling parameter’. Critically, these length scale considerations determine which asymptotic techniques we may employ to investigate the problem. Namely, when $\gamma _m \ll 1$ , we use ray-tracing results derived from a WKBJ approximation whereas when $\gamma _m \gg 1$ , we approximate the filament as a delta-function.

To make the dependence on $\gamma _m$ explicit, and to aid the asymptotic analysis by expressing the filament as an order 1 function of an order 1 parameter, we divide (3.4) by $\textit{Ro}_{\!f}$ to give

(3.7a) \begin{equation} \left [\textit{Ro}_{\!f}^{-1}\left (1 + \frac {\partial ^2\,}{\partial \tilde {t}^2}\right ) - \gamma _m\frac {\partial ^2\,}{\partial \tilde {x}^2} + \mathcal{F}(\tilde {x})\right ]\mathcal{X} = 0. \end{equation}

Applying the same manipulations to the time-dependent Schrödinger equation (2.9) gives

(3.7b) \begin{equation} \left [-2\mathrm{i}\textit{Ro}_{\!f}^{-1}\frac {\partial \,}{\partial \tilde {t}} - \gamma _m\frac {\partial ^2\,}{\partial \tilde {x}^2} + \mathcal{F}(\tilde {x})\right ]\mathcal{A} = 0 \end{equation}

where the filament Rossby number can be absorbed into a rescaled time. Thus, in YBJ theory, $\gamma _m$ is the only dynamically interesting parameter (e.g. Young & Ben Jelloul Reference Young and Jelloul1997; Danioux, Vanneste & Bühler Reference Danioux, Vanneste and Bühler2015; Asselin & Young Reference Asselin and Young2019; Conn, Callies & Lawrence Reference Conn, Callies and Lawrence2025). This parameter dependence is neatly illustrated by the radiation time scale $T_{10}$ , that is, the time scale over which the solution decays at the centre of the filament. We plot the radiation time scale, computed from numerical simulations of a Gaussian filament (see § 3.5), in units of inertial periods (figure 2 b). There is a simple monotonic structure along lines of $\gamma _m = \mathrm{constant}$ . Indeed, (3.7b ) implies a linear dependence on $\textit{Ro}_{\!f}$ in log space. Whereas, on other lines with $\gamma _m$ varying, we observe a more complex non-monotonic structure. The same plot for the full Klein–Gordon equation (figure 2 a) displays similar parameter dependence but with more complicated behaviour when $\textit{Ro}_{\!f} \gt 1$ .

Figure 2. Radiation time $\tilde {T}_{10}/2\pi := fT_{10}/2\pi$ (3.2) in inertial periods computed from numerical simulations of a Gaussian filament (§ 3.5) as a function of $\gamma _m$ and $\textit{Ro}_{\!f}$ for both the Klein–Gordon (a) and YBJ (b) problems. Lines of constant $Bu_m$ and $\alpha _m^2$ are overlaid in grey and brown, respectively. The diagram inset in (b) indicates the distinguished limits summarised in table 3. White regions are excluded as, for these parameter values, waves radiated from the filament loop around the finite numerical domain and return to the filament before the $T_{10}$ criterion is met.

The critical importance of the parameter $\gamma _m$ motivates using it as one of the two non-dimensional numbers spanning the parameter space. Writing the Klein–Gordon equation in the form (3.7a ) suggests using the filament Rossby number as the second non-dimensional number and by default this is what we do (see, e.g. figures 2 and 4). However, this is not the only choice. Indeed, there are at least two other interesting choices that warrant comment. Of the four dimensional parameters defining this problem two are velocity scales, $c_m$ and $\Delta \kern-1pt V$ , and the length and time scales may be combined to form a third velocity scale, $fL_f$ . We can define three non-dimensional parameters that are independent of one of these velocity scales. The filament Rossby number (which is independent of $c_m$ ) is the first of these, the second is the filament Burger number (which is independent of $\Delta \kern-1pt V$ ) and the third is

(3.8) \begin{equation} \alpha _m := \frac {\Delta \kern-1pt V}{c_m} \end{equation}

that is independent of $f$ and $L_f$ . These three parameters are related through $\gamma _m$ :

(3.9) \begin{equation} Bu_m \gamma _m^{-1} \equiv \textit{Ro}_{\!f} \equiv \alpha _m^2 \gamma _m. \end{equation}

By fixing one of these three parameters and then sending $\gamma _m \to 0$ or $\gamma _m \to \infty$ , six physically interesting distinguished limits can be reached. For example, suppose that we wish to study the strongly stratified limit, then we fix the values of $f$ , $L_f$ and $\Delta \kern-1pt V$ before sending $c_m \to \infty$ . However, this is simply the distinguished limit $\gamma _m \to \infty$ with $\textit{Ro}_{\!f}$ fixed. Alternatively, the sharp filament limit $L_f \to 0$ with $f$ , $c_m$ and $\Delta \kern-1pt V$ fixed is the distinguished limit $\gamma _m \to \infty$ with $\alpha _m$ fixed. The six distinct distinguished limits are summarised in table 3 and figure 2(b).

Table 3. Distinguished limits achieved by fixing three of the four dimensional parameters and sending the fourth to $0$ or $\infty$ .

3.3. The WKBJ approximation – $\gamma _m \ll 1$

We begin by considering the regime in which the WKBJ approximation is valid, namely $\gamma _m \ll 1$ . Here, we recall some standard ray-tracing results. However, in Appendix A we give a formal WKBJ derivation in the limit $\gamma _m \to 0$ considering the three distinguished limits summarised in table 3. The solution is expressed in terms of a slowly varying amplitude and rapidly varying phase

(3.10) \begin{equation} \mathcal{X} = A\mathrm{e}^{\mathrm{i}\theta } + \mathrm{c.c.} \end{equation}

Then, with non-dimensional frequency $\tilde {\omega } := -\partial \theta / \partial \tilde {t}$ , which given the steady background flow is conserved along rays, and local wavenumber $\tilde {k} := \partial \theta / \partial \tilde {x}$ , (3.7a ) admits the dispersion relation

(3.11a) \begin{equation} \tilde {\omega }^2 - 1 = \textit{Ro}_{\!f}\big [\mathcal{F}(\tilde {x}) + \gamma _m \tilde {k}^2\big ] \end{equation}

and the group velocity is

(3.11b) \begin{equation} \tilde {c}_g := \frac {\partial \tilde {\omega }}{\partial \tilde {k}} = \textit{Ro}_{\!f}\gamma _m\frac {\tilde {k}}{\tilde {\omega }}. \end{equation}

Notably, as the waves radiate away from the centre of the filament, to conserve $\tilde {\omega }$ , $\tilde {k}$ becomes $O(\gamma _m^{-1/2})$ , i.e. the waves have spatial scale $\sqrt {\gamma _m}L_f$ . The YBJ versions of (3.11) are

(3.12a,b) \begin{align} \tilde {\omega }^{(\textit{YBJ}\kern1pt)} = 1 + \frac {1}{2} \textit{Ro}_{\!f}\big [\mathcal{F}(\tilde {x}) + \gamma _m \tilde {k}^2\big ], \quad \tilde {c}_g^{(\textit{YBJ}\kern1pt)} = \textit{Ro}_{\!f}\gamma _m \tilde {k}. \end{align}

In both cases the transport equation is

(3.13) \begin{equation} \frac {\partial \,}{\partial \tilde {t}}(A^2) + \frac {\partial \,}{\partial \tilde {x}}\big (\tilde {c}_gA^2\big) = 0. \end{equation}

We note that $\tilde {\omega }^{(\textit{YBJ}\kern1pt)}$ is the expansion of $\tilde {\omega }$ to $O(\textit{Ro}_{\!f})$ and that $\tilde {c}_g^{(\textit{YBJ}\kern1pt)}$ is (3.11b ) but using the leading-order, i.e. $\tilde {\omega } = 1$ , expression for the frequency. Consequently, if we apply the YBJ approximation when high frequency waves are being radiated then the group velocity will be overpredicted and the wave energy radiated too rapidly.

We now compute expressions for $T_{10}$ and $T_{10}^{(\textit{YBJ}\kern1pt)}$ under the WKBJ approximation. We start by computing the travel time of a ray originating at $\tilde {x} = \tilde {x}_0$ . Using the dispersion relation (3.11a ) to eliminate $\tilde {k}$ from the group velocity (3.11b ), we have

(3.14a) \begin{equation} {\tilde {c}_g}^2 = \left (\frac {\,\mathrm{d}\tilde {x}}{\,\mathrm{d}\tilde {t}}\right )^2 = \textit{Ro}_{\!f}\gamma _m \frac {\tilde {\omega }^2 - (1 + \textit{Ro}_{\!f}\mathcal{F}(\tilde {x}))}{\tilde {\omega }^2}. \end{equation}

As the initial conditions are uniform, we have $\tilde {k} = 0$ at $\tilde {t} = 0$ and, thus, the frequency satisfies

(3.14b) \begin{equation} \tilde {\omega }^2 = 1 + \textit{Ro}_{\!f}\mathcal{F}(\tilde {x}_0). \end{equation}

It then follows that the travel time $\tilde {\tau }$ of a ray from $\tilde {x}_0 \gt 0$ to $\tilde {x} \gt \tilde {x}_0$ is

(3.14c) \begin{equation} \tilde {\tau }(\tilde {x};\tilde {x}_0) = \textit{Ro}_{\!f}^{-\frac {1}{2}}\gamma _m^{-\frac {1}{2}}\int _{\tilde {x}_0}^{\tilde {x}}\sqrt {\frac {\textit{Ro}_{\!f}^{-1} + \mathcal{F}(\tilde {x}_0)}{\mathcal{F}(\tilde {x}_0) - \mathcal{F}(\tilde {x}')}}\,\mathrm{d}\tilde {x}'. \end{equation}

The equivalent YBJ calculation gives

(3.15) \begin{equation} \tilde {\tau }^{(\textit{YBJ}\kern1pt)}(\tilde {x};\tilde {x}_0) = \textit{Ro}_{\!f}^{-1}\gamma _m^{-\frac {1}{2}}\int _{\tilde {x}_0}^{\tilde {x}}\frac {1}{\sqrt {\mathcal{F}(\tilde {x}_0) - \mathcal{F}(\tilde {x}')}}\,\mathrm{d}\tilde {x}', \end{equation}

which is the leading-order term of (3.14c ) in the small $\textit{Ro}_{\!f}$ limit.

To compute $T_{10}$ , we must consider rays close to the centre of the filament where $\mathcal{F}(\tilde {x}) = \mathcal{F}(0) + ({1}/{2})\mathcal{F}''(0)\tilde {x}^2 + \mathrm{h.o.t.}$ . Here, $\mathcal{F}'' := \partial ^2 \mathcal{F}/\partial \tilde {x}^2$ . Inserting this into (3.14c ) we get

(3.16) \begin{equation} \textit{Ro}_{\!f}^{\frac {1}{2}}\gamma _m^{\frac {1}{2}}\tilde {\tau }(\tilde {x};\tilde {x}_0) \approx \int _{\tilde {x}_0}^{\tilde {x}}\sqrt {\frac {\textit{Ro}_{\!f}^{-1} + \mathcal{F}(0)}{\frac {1}{2}\mathcal{F}''(0)\left (\tilde {x}_0^2 - \tilde {x}'^2\right )}}\,\mathrm{d}\tilde {x}' = \sqrt {\frac {\textit{Ro}_{\!f}^{-1} + \mathcal{F}(0)}{\frac {1}{2}\mathcal{F}''(0)}}\mathrm{arccosh}\left (\frac {\tilde {x}}{\tilde {x}_0}\right )\!. \end{equation}

From this expression we see that rays originating near the centre of the filament fan out uniformly in the sense that $\tau (\tilde {x};\tilde {x_0}) = \tau (\tilde {x}/\tilde {x}_0)$ . The transport equation (3.13) implies that the total wave amplitude squared between any two rays is conserved. That is, if $a(\tilde {t})$ and $b(\tilde {t})$ are two rays with initial positions $a_0$ and $b_0$ then

(3.17) \begin{equation} \frac {\mathrm{d}\,}{\,\mathrm{d}\tilde {t}}\int _{a(\tilde {t})}^{b(\tilde {t})}A^2\,\mathrm{d}\tilde {x} = 0\ \; \implies\; \frac {\mathrm{d}\,}{\mathrm{d}\tilde {t}}\int _{a_0}^{b_0}A^2\frac {\mathrm{d}\tilde {x}}{\mathrm{d}\tilde {x}_0}\mathrm{d}\tilde {x}_0 = 0\ \; \implies\, \int _{a_0}^{b_0}A^2\frac {\mathrm{d}\tilde {x}}{\mathrm{d}\tilde {x}_0}\mathrm{d}\tilde {x}_0 = \mathrm{const.} \end{equation}

Not only does this hold for all rays but these rays are fanning out uniformly and, hence, $\,\mathrm{d}\tilde {x}/\mathrm{d}\tilde {x}_0$ is independent of $\tilde {x}_0$ . It therefore follows that, for the amplitude to drop by a factor of 10, we require $\tilde {x} / \tilde {x}_0 = 10^2 = 100$ . Thus, we have

(3.18a) \begin{equation} \tilde {T}_{10} := fT_{10} = \textit{Ro}_{\!f}^{-\frac {1}{2}}\gamma _m^{-\frac {1}{2}}\sqrt {\frac {\textit{Ro}_{\!f}^{-1} + \mathcal{F}(0)}{\frac {1}{2}\mathcal{F}''(0)}}\mathrm{arccosh}(100) \end{equation}

and

(3.18b) \begin{equation} \tilde {T}^{(\textit{YBJ}\kern1pt)}_{10} := fT^{(\textit{YBJ}\kern1pt)}_{10} = \textit{Ro}_{\!f}^{-1}\gamma _m^{-\frac {1}{2}}\frac {1}{\sqrt {\frac {1}{2}\mathcal{F}''(0)}}\mathrm{arccosh}(100).\end{equation}

Comparing these ray-tracing results for the full Klein–Gordon and YBJ versions of the problem, we conclude that we may divide the $\gamma _m \ll 1$ regime into two subregimes (IA and IB) dependent on the Rossby number. Regime IA is the regime $\gamma _m \ll 1$ , $\textit{Ro}_{\!f} \ll 1$ where the YBJ approximation is valid. On the other hand, regime IB is the regime $\gamma _m \ll 1$ , $\textit{Ro}_{\!f} \gg 1$ and the YBJ approximation is not valid. These conclusions regarding the validity of the YBJ approximation are corroborated by the formal asymptotic analysis in Appendix A. In the two subregimes the radiation time scale, $\tilde {T} := fT$ , scales as

(3.19) \begin{equation} \tilde {T} \sim \begin{cases} \gamma _m^{-\frac {1}{2}}\textit{Ro}_{\!f}^{-1}, & \textit{Ro}_{\!f} \ll 1, \\ \gamma _m^{-\frac {1}{2}}\textit{Ro}_{\!f}^{-\frac {1}{2}}, & \textit{Ro}_{\!f} \gg 1. \end{cases} \end{equation}

In both subregimes the radiation time scale is a decreasing function of both the Rossby number and tunnelling parameter.

3.4. The delta-function limit – $\gamma _m \gg 1$

We now consider $\gamma _m \gg 1$ . From (3.7a ), we see that, to leading order, NIWs will not vary on the filament scale or the dispersive term, $\gamma _m\partial ^2/\partial \tilde {x}^2$ , would be unbalanced. Thus, our previous choice of spatial non-dimensionalisation, by $L_f$ , is no longer appropriate. We choose a new non-dimensionalisation by defining

(3.20a,b) \begin{align} \check {x} := \frac {x}{L_m}, \quad L_m := \gamma _m L_f \equiv \frac {c_m^2}{f\Delta \kern-1pt V}. \end{align}

Multiplying (3.7a ) through by $\gamma _m$ we have

(3.21a) \begin{equation} \left [\alpha _m^{-2}\left (1 + \frac {\partial ^2\,}{\partial \tilde {t}^2}\right ) - \frac {\partial ^2\,}{\partial \check {x}^2} + \gamma _m\mathcal{F}(\gamma _m\check {x})\right ]\mathcal{X} = 0. \end{equation}

This choice of non-dimensionalisation has achieved two things. Firstly, we have cleared the coefficient of the spatial derivatives and, secondly, the filament is now in the form of a nascent delta-function. As $\gamma _m \to \infty$ ,

(3.21b) \begin{equation} \gamma _m\mathcal{F}(\gamma _m\check {x}) \to \delta (\check {x}). \end{equation}

We first study the distinguished limit in which $\alpha _m$ is held constant. This is the most natural interpretation of the delta-function limit as it corresponds to sending $L_f \to 0$ with $c_m$ , $f$ and $\Delta \kern-1pt V$ held fixed. We consider the other distinguished limits summarised in table 3 at the end of the section.

Integrating over the delta-function and requiring that the solution is even, we have the jump condition $\partial \mathcal{X}/\partial \check {x} = ({1}/{2}) \mathcal{X}$ at $\check {x} = 0^+$ . Redimensionalising, the jump condition is $L_m\partial \mathcal{X} / \partial x = ({1}/{2})\mathcal{X}$ . Thus, we may interpret $L_m$ as the length scale imposed on the problem by the delta-function filament through the jump condition.

In Appendix B we derive a solution in the delta-function limit. For $\tilde {t} \lt \alpha _m^{-1}|\check {x}|$ , the solution is simply $\mathcal{X}(\check {x},\tilde {t}) = \cos {\tilde {t}}$ . Note that in ( $\check {x}$ , $\tilde {t}$ ) coordinates, the characteristics of (3.21a ) have slope $\pm \alpha _m^{-1}$ and so the points $|\check {x}| \gt \alpha _m\tilde {t}$ are beyond the influence of the filament. For $\tilde {t} \geqslant \alpha _m^{-1}|\check {x}|$ , the solution is given by

(3.22a) \begin{equation} \mathcal{X}(\check {x},\tilde {t}) = \cos {\tilde {t}} - \frac {1}{2}\alpha _m \int _0^{\tilde {t} - \alpha _m^{-1}|\check {x}|}J_0\left (\sqrt {(\tilde {t} - \tilde {t}')^2 - \alpha _m^{-2}\check {x}^2}\right )\mathcal{X}(0,\tilde {t}')\,\mathrm{d}\tilde {t}', \end{equation}

where $J_0$ is the zeroth-order Bessel function of the first kind. The solution at any point may be computed from the history of the solution at the centre of the filament. The solution at the centre of the filament is

(3.22b) \begin{equation} \mathcal{X}(0,\tilde {t}) = \frac {2}{\pi }\int _0^\infty \frac {1}{1 + u^2}\cos \left (\tilde {t}\sqrt {1 + \frac {1}{4}\alpha _m^2u^2}\right )\mathrm{d}u. \end{equation}

For fixed $\tilde {t}$ , as $\alpha _m \to 0$ , this reduces to the elementary integral $(2/\pi )\cos {\tilde {t}}\int _0^\infty (1 + u^2)^{-1}\mathrm{d}u = \cos {\tilde {t}}$ . Evaluating (3.22b ) numerically, we find that, for larger $\alpha _m$ , i.e. stronger filaments, the decay at the centre of the filament is more rapid (figure 3). Furthermore, we find very good agreement between the analytic delta-function solution and numerical simulations (see § 3.5) with $\gamma _m = 5$ .

Figure 3. Solutions $\mathcal{X}(0,\tilde {t})$ to the Klein–Gordon equation at $x = 0$ for (a) $\alpha _m^2 = 0.25$ and (b) $\alpha _m^2 = 4$ . Solid lines are from the numerical solutions. The dashed black line shows the analytic solution (3.22b ) and the dashed grey line shows the stationary phase approximation (3.23). The time axis is in inertial periods.

For large $\tilde {t}$ , (3.22b ) is amenable to the stationary phase approximation and we find that

(3.23) \begin{equation} \mathcal{X}(0,\tilde {t}) \sim \sqrt {\frac {8}{\pi \alpha _m^2 \tilde {t}}}\cos \left (\tilde {t} + \frac {\pi }{4}\right )\!. \end{equation}

Here, large $\tilde {t}$ means $\tilde {t} \gg \max (1,\alpha _m^{-2})$ . For $O(1)$ values of $\alpha _m$ , the stationary phase approximation is very good on an inertial time scale (e.g. figure 3 b where $\alpha _m = 2$ ). Furthermore, when $\tilde {t}$ is large, the solution is an inertial oscillation with a slowly decaying amplitude. This is a situation in which the YBJ approximation should be expected to perform well. However, this provides no guarantee that the YBJ approximation will correctly capture the small $\tilde {t}$ behaviour.

To compare to § 3.3, we again compute the radiation time scale. Here we have $\tilde {T} \sim \alpha _m^{-2} \equiv \textit{Ro}_{\!f}^{-1}\gamma _m$ , which is once again a decreasing function of the Rossby number but now an increasing function of the tunnelling parameter. More precisely,

(3.24) \begin{equation} \tilde {T}_{10} = 10^2 \times \frac {8}{\pi }\alpha _m^{-2} = \frac {800}{\pi }\textit{Ro}_{\!f}^{-1}\gamma _m. \end{equation}

However, if $\alpha _m^2$ is large then this result may be inaccurate as $\tilde {T}_{10}$ may be too small for the stationary phase approximation to be valid.

Finally, we consider the other distinguished limits. In particular, we consider $\gamma _m \to \infty$ with $\textit{Ro}_{\!f}$ fixed. The case with $Bu_m$ fixed can be handled in exactly the same way. The complicating factor is that $\alpha _m^{-2} \equiv \textit{Ro}_{\!f}^{-1} \gamma _m \to \infty$ as $\gamma _m \to \infty$ . The solution is to utilise substitution (2.10), $\mathcal{X} = \mathcal{A}\mathrm{e}^{-\mathrm{i}\tilde {t}} + \mathrm{c.c.}$ ,

(3.25a) \begin{equation} \left [\textit{Ro}_{\!f}^{-1}\gamma _m\left (\frac {\partial ^2\,}{\partial \tilde {t}^2} - 2\mathrm{i}\frac {\partial \,}{\partial \tilde {t}}\right ) - \frac {\partial ^2\,}{\partial \check {x}^2} + \gamma _m\mathcal{F}(\gamma _m\check {x})\right ]\mathcal{A} = 0, \end{equation}

and then rescale time. We define $\check {t} := \textit{Ro}_{\!f}\gamma _m^{-1}\tilde {t} \equiv \alpha _m^2\tilde {t}$ such that

(3.25b) \begin{equation} \left [\textit{Ro}_{\!f}\gamma _m^{-1}\frac {\partial ^2\,}{\partial \check {t}^2} - 2\mathrm{i}\frac {\partial \,}{\partial \check {t}} - \frac {\partial ^2\,}{\partial \check {x}^2} + \gamma _m\mathcal{F}(\gamma _m\check {x})\right ]\mathcal{A} = 0. \end{equation}

The leading-order solution is thus given by the delta-function solution to the YBJ equation, i.e. the $\partial ^2/\partial \check {t}^2$ term does not appear until $O(\gamma _m^{-1})$ . That being said, the full Klein–Gordon delta-function solution is still informative in this case as it just corresponds to the inclusion of a higher-order term. Finally, we note that the temporal scaling $\check {t} \equiv \alpha _m^2\tilde {t}$ introduced here in the YBJ solution is the same scaling that arose in the stationary phase approximation (3.23).

Importantly, this scaling analysis means that the delta-function approximation not only applies in the distinguished limit $L_f \to 0$ with $c_m$ , $f$ and $\Delta \kern-1pt V$ fixed but also in the limit $c_m \to \infty$ with $f$ , $L_f$ and $\Delta \kern-1pt V$ fixed. In both cases, the radiation time scale is $\tilde {T} \sim \alpha _m^2 = \textit{Ro}_{\!f}\gamma _m^{-1}$ , although in the latter case this time scale is guaranteed to be large. In other words, when the stratification is sufficiently strong or the vertical wavelength sufficiently large, the waves, to leading order, have no structure on the filament scale and it takes a very long time for the presence of the filament to be felt regardless of the strength of the filament.

3.5. Numerical solutions

To explore the intermediate regime $\gamma _m = O(1)$ , validate our scalings for the decay at the centre of the filament and further assess the validity of the YBJ approximation, we numerically solve the Klein–Gordon and time-dependent Schrödinger equations for a Gaussian filament $\mathcal{F}(\eta ) = \exp (-\eta ^2/2)/\sqrt {2\pi }$ (figure 1 a). Details of the numerical schemes can be found in Appendix C. We present the results using the non-dimensionalisation of § 3.3, $\tilde {t} = ft$ , $\tilde {x} = x/L_f$ , and use $\gamma _m$ and $\textit{Ro}_{\!f}$ as coordinates for the parameter space. For the numerical solutions, the ‘unbounded’ domain has periodic boundary conditions at $\tilde {x} = \pm 5000$ and we use $2^{18}$ grid points.

Table 4. Spatial non-dimensionalisations, approximations and scalings for the decay time scale at the centre of the filament in the different regimes of the unbounded radiation problem.

From a suite of 12 Klein–Gordon simulations (figure 4), we observe, for fixed $\textit{Ro}_{\!f}$ , a transition in behaviour from wave radiation for small $\gamma _m$ to a more spatially coherent decaying response for large $\gamma _m$ . We also observe that, for fixed $\gamma _m$ , increasing $\textit{Ro}_{\!f}$ leads to more rapid decay in the filament. However, the more interesting dependence is on $\gamma _m$ . The scalings derived earlier and summarised in table 4 predict that the radiation time scale is a decreasing function of $\gamma _m$ , $\tilde {T}\sim \gamma _m^{-1/2}$ , for small $\gamma _m$ and an increasing function $\tilde {T}\sim \gamma _m$ for large $\gamma _m$ . We should therefore expect the most rapid radiation for some intermediate value of $\gamma _m$ . This is exactly what we observe in the simulations with the most rapid radiation occurring for the $\gamma _m = 0.1$ simulations. For a large Rossby number, this radiation can be extremely rapid. For example, for $\gamma _m = 0.1$ , $\textit{Ro}_{\!f} = 10$ (figure 4 b), the solution at the centre of the filament becomes negligibly small in less than two inertial periods. It should be noted that these values are quite reasonable in some regions of the world’s oceans. With $\textit{Ro}_{\!f} = 10$ the geostrophic velocity gradient at the centre of the filament is $\partial V/\partial x = (10/\sqrt {2\pi }) f \approx 4.0 f$ , well within the range of values simulated by Schlichting et al. (Reference Schlichting, Qu, Kobashi and Hetland2023) for example.

Figure 4. Hovmöller plots of $\mathcal{X}$ showing the radiation by an unbound cyclonic filament for varying $\textit{Ro}_{\!f}$ and $\gamma _m$ . Dashed white lines indicate rays travelling at the mode speed $\tilde {t} = \pm \tilde {x}/\sqrt {\textit{Ro}_{\!f}\gamma _m}$ . The time axis is in inertial periods.

Plotting $\tilde {T}_{10}$ , computed by running numerical simulations until the stopping criterion (3.2) is reached, as a function of $\gamma _m$ for various $\textit{Ro}_{\!f}$ (figures 2 and 5 a) once again highlights the monotonic dependence on Rossby number but non-monotonic dependence on the tunnelling parameter. However, there is a significant and sharp drop in $\tilde {T}_{10}$ as we transition to intermediate values of $\gamma _m$ . In particular, the transition into regime II, which, for $\textit{Ro}_{\!f} \ll 1$ , occurs at $\gamma _m \approx 0.4$ , is abrupt and $\tilde {T}_{10}$ is very sensitive to the value of $\gamma _m$ around this point. To validate our predictions for the radiation time scale, we should be able to make the curves collapse by appropriately rescaling the axes. Therefore, we let $Y = \textit{Ro}_{\!f}^a \tilde {T}_{10}$ and $X = \textit{Ro}_{\!f}^b \gamma _m$ , where $a$ and $b$ are exponents to be determined. In each of the three subregimes we have scalings for $\tilde {T}$ in terms of $\gamma _m$ and $\textit{Ro}_{\!f}$ (table 4). Eliminating $\gamma _m$ gives the following relationships between $Y$ and $X$ for regimes IA, IB and II respectively:

(3.26a–c) \begin{align} Y^2 \sim \textit{Ro}_{\!f}^{(2a + b - 2)}X^{-1}, \quad Y^2 \sim \textit{Ro}_{\!f}^{(2a + b - 1)}X^{-1}, \quad Y \sim \textit{Ro}_{\!f}^{(a - b - 1)}X. \end{align}

To make the curves collapse, we must eliminate $\textit{Ro}_{\!f}$ . For small $\textit{Ro}_{\!f}$ , subregimes IA and II apply and, thus, we require $a = 1$ and $b = 0$ . Whereas for large $\textit{Ro}_{\!f}$ , regimes IB and II apply and we require $a = 2/3$ and $b = -1/3$ . Applying these rescalings to the axes (figures 5 b and 5 c) we find that the curves do indeed collapse. Furthermore, by plotting the asymptotic predictions for $\tilde {T}_{10}$ from (3.18) and (3.24) on top of the numerical results we find excellent agreement.

Figure 5. Time scale for the decay of velocity at the centre of the filament as a function of $\gamma _m$ for various $\textit{Ro}_{\!f}$ . In (a) the time axis is in inertial periods. By appropriately rescaling the axes for small (b) and large (c) $\textit{Ro}_{\!f}$ the curves may be made to collapse. Dashed black lines indicate the radiation time scales (3.18) computed for regimes IA and IB. Dotted black lines indicate the radiation time scale (3.24) predicted by the stationary phase approximation for regime II.

Figure 6. Hovmöller plots of the difference between the YBJ and Klein-Gordon solutions for varying $\textit{Ro}_{\!f}$ and $\gamma _m$ . Dashed white lines indicate rays travelling at the mode speed $\tilde {t} = \pm \tilde {x}/\sqrt {\textit{Ro}_{\!f}\gamma _m}$ . The time axis is in inertial periods.

To further assess the validity of the YBJ approximation as a function of $\gamma _m$ and $\textit{Ro}_{\!f}$ , we plot the difference between the YBJ approximation and Klein–Gordon solutions in figure 6. For small $\textit{Ro}_{\!f}$ , the YBJ approximation performs well for all values of $\gamma _m$ . For larger $\textit{Ro}_{\!f}$ , when $\gamma _m$ is small, we observe differences. Notably, the YBJ solution radiates waves outside of the region bounded by the mode speed, $\tilde {t} = \pm \tilde {x}/\sqrt {\textit{Ro}_{\!f}\gamma _m}$ (white dashed lines in figure 6). This is a manifestation not just of the overprediction of the group velocity highlighted in § 3.3 but also of the fact that the Klein–Gordon equation is hyperbolic and propagates information at a finite speed (the mode speed) whereas the time-dependent Schrödinger equation is only first-order in time and can propagate information at any speed. However, when $\gamma _m = 10$ , the YBJ approximation performs well even for $\textit{Ro}_{\!f} = 10$ . This is consistent with our analysis in § 3.4 as for $\gamma _m = \textit{Ro}_{\!f} = 10$ , $\alpha _m^2 = 1$ is not large. Furthermore, the YBJ prediction for the radiation time scale (figure 2) shows good agreement with the Klein–Gordon prediction when $\gamma _m \gg 1$ and $\alpha _m \lt 1$ . These results and the earlier analysis suggest that the YBJ approximation is valid when $\gamma _m \ll 1$ and $\textit{Ro}_{\!f} \ll 1$ or $\gamma _m \gg 1$ and $\alpha _m \ll 1$ . If we assume that the validity criterion smoothly transitions across $\gamma _m = O(1)$ , as figure 2 suggests it should, then we can say that the YBJ approximation should be valid whenever

(3.27) \begin{equation} \textit{Ro}_{\!f} \ll 1 + \gamma _m. \end{equation}

3.6. Summary

In this section we considered the lateral radiation of a single vertical mode with no initial horizontal structure out of a cyclonic filament into an unbounded domain. The problem is described by a 2-D parameter space. There are two regimes distinguished by the tunnelling parameter $\gamma _m$ . For small $\gamma _m$ , a WKBJ approximation is permissible and we identify the filament Rossby number $\textit{Ro}_{\!f}$ as the most useful choice for the second non-dimensional parameter. In particular, for small $\textit{Ro}_{\!f}$ , the YBJ approximation is valid. For large $\gamma _m$ , the filament may be treated as a delta-function and we derive an analytic solution. In this regime the most useful choice for the second parameter is $\alpha _m^2 = \textit{Ro}_{\!f}/\gamma _m$ , which determines the validity of the YBJ approximation. A stationary phase approximation reveals that at large times the solution in the filament decays as $\alpha _m^{-1}\tilde {t}^{-({1}/{2})}$ . A particularly important result comes from the distinguished limit $\gamma _m \to \infty$ , $\alpha _m^2 \to 0$ with $\textit{Ro}_{\!f}$ held constant. In this limit, which may be achieved by sending the mode speed to infinity, we find that the filament, however large the Rossby number, is only felt by the inertial oscillations on very long time scales. Finally, we find that the decay at the centre of the filament is fastest for intermediate values ( $0.1 - 1$ ) of $\gamma _m$ and, for large but realistic values of the filament Rossby number, this decay can occur on inertial time scales. This ‘Goldilocks’ effect where the most efficient radiation for a given background flow (i.e. fixed values of $\Delta \kern-1pt V$ , $L_f$ , $f$ and $N$ ), occurs for waves with just-the-right vertical scale is a recurring theme both of this paper and of previous studies of NIW-mean flow interactions (Balmforth et al. Reference Balmforth, Llewellyn and Young1998; Klein & Llewellyn Smith Reference Klein and Smith2001; Danioux et al. Reference Danioux, Vanneste and Bühler2015).

4.1. Reduction to the radiation problem

If we may ignore the boundary conditions then this problem is equivalent to the unbounded radiation problem considered in § 3. However, we must account for the change in the background flow away from the filament. Separating out the constant anticyclonic part of the background flow, the Klein–Gordon equation (2.5) can be written as

(4.3a) \begin{equation} \left [\left (1 + \textit{Ro}_{ac} + \frac {\partial ^2\,}{\partial \tilde {t}^2}\right ) - \frac {c_m^2}{f^2}\frac {\partial ^2\,}{\partial x^2} + \textit{Ro}_{\!f} \mathcal{F}\left (\frac {x}{L_f}\right )\right ]\mathcal{X} = 0. \end{equation}

Now we rescale time by defining $\hat {t} := \sqrt {1 + \textit{Ro}_{ac}}\tilde {t} \equiv f_{\textit{eff}}^{(ac)}t$ , where

(4.3b) \begin{equation} f_{\textit{eff}}^{(ac)} := f\sqrt {1 + \textit{Ro}_{ac}} \lt f \end{equation}

is the effective Coriolis frequency (2.2b ) of the anticyclonic region. Finally, we proceed as in § 3.3 by non-dimensionalising space by $L_f$ , $\tilde {x} := x/L_f$ , and dividing (4.3a ) by $\textit{Ro}_{\!f}$ to get

(4.3c) \begin{equation} \left [\frac {1 + \textit{Ro}_{ac}}{\textit{Ro}_{\!f}}\left (1 + \frac {\partial ^2\,}{\partial \hat {t}^2}\right ) - \gamma _m\frac {\partial ^2\,}{\partial \tilde {x}^2} + \mathcal{F}(\tilde {x})\right ]\mathcal{X} = 0. \end{equation}

This is identical to (3.7a ) up to a redefinition of the Rossby number. Crucially, the key dimensionless number, the tunnelling parameter $\gamma _m$ , appears unchanged.

Figure 7. Hovmöller plots showing the evolution of NIWs interacting with a cyclonic filament in an otherwise anticyclonic flow. The parameters have been chosen to be dynamically equivalent to figure 4( a – d ). Dashed white lines indicate rays travelling at the free wave speed. The left time axis is normalised by the period of a wave with frequency equal to the effective Coriolis frequency of the anticyclonic region. The right time axis is in inertial periods.

Consider the top row of figure 4 in which we looked at the unbounded radiation problem with $\textit{Ro}_{\!f} = 10$ for various $\gamma _m$ . We can reproduce those results in the current set-up, assuming the boundary conditions are not important, by using the same $\gamma _m$ values and choosing $\textit{Ro}_{ac}$ and $\textit{Ro}_{\!f}$ such that $\textit{Ro}_{\!f} / (1 + \textit{Ro}_{ac}) = 10$ . We do this for $\textit{Ro}_{ac} = -1/6$ , $\textit{Ro}_{\!f} = 50/6\; \implies \xi = 25$ (figures 7 a– 7 d) and $\textit{Ro}_{ac} = -({1}/{2})$ , $\textit{Ro}_{\!f} = 5 \implies \xi = 5$ (figures 7 e–7 h). We observe that figures 7(a) and 7(b) are essentially identical to figures 4(a) and 4(b). Note that this comparison is aided by the fact that in figure 4 we plot the solutions for $|\tilde {x}| \leqslant 25$ and that in figures 7(a)–7(d) we impose periodic boundary conditions at $\tilde {x} = \pm 25$ . Furthermore, figure 7(e) is the same as figures 4(a) and 7(a) but zoomed in along the $\tilde {x}$ axis as the periodic boundary conditions are enforced at $\tilde {x} = \pm 5$ . This similitude is a manifestation of the dynamical equivalence of these two problems when the boundaries are not influencing the flow.

In order for the boundary conditions to not affect the flow we require that $L_x$ be much larger than the important length scales of the radiation problem. First, we have $L_x \gg L_f$ by construction. Another important length scale, derived from the mode speed, is the Rossby radius, $c_m/f$ , which determines the distance that the Klein–Gordon equation can propagate on an inertial time scale. If $L_x \gg c_m/f$ then it takes a long time for filament to feel the influence of the boundary conditions. Moving from figures 7(a) to 7(d) or from figures 7(e) to 7(h), corresponds to an increase in $c_m$ with the other parameters held constant. As $c_m$ is increased, the solutions diverge from figures 4(a) to 4(d) and tend towards uniform inertial oscillations. Notably, these oscillations are inertial, i.e. they have dimensional period $2\pi /f$ not $2\pi / f_{\textit{eff}}^{(ac)}$ , and do not decay in time. The other key length scale in the radiation problem was $L_m$ and we find that the ratio of $L_m$ to $L_x$ is the key non-dimensional number in the following sections.

4.2. Horizontal modes

We now consider the cases for which the boundary conditions are important. Here and throughout the rest of the paper, it is most convenient to non-dimensionalise space by the half-width of the domain $L_x$ . Therefore, we define

(4.4a) \begin{equation} \hat {x} := \frac {x}{L_x} \end{equation}

and we non-dimensionalise the lateral geostrophic shear with a factor of $|\textit{Ro}_{ac}|$ , i.e.

(4.4b) \begin{equation} \frac {\partial \hat {V}}{\partial \hat {x}} := |\textit{Ro}_{ac}|^{-1}\frac {1}{f}\frac {\partial V}{\partial x} = -1 + 2\xi \mathcal{F}(\xi \hat {x}), \end{equation}

where again $\xi := L_x/L_f$ . This implies that we have non-dimensionalised the geostrophic velocity using the velocity scale $({1}/{2})\Delta \kern-1pt V$ . Once more using the temporal non-dimensionalisation $\tilde {t} := ft$ , we write the Klein–Gordon equation as

(4.5) \begin{equation} \left [|\textit{Ro}_{ac}|^{-1}\left (1 + \frac {\partial ^2\,}{\partial \tilde {t}^2}\right ) - \varGamma _m \frac {\partial ^2\,}{\partial \hat {x}^2} + \frac {\partial \hat {V}}{\partial \hat {x}}\right ]\mathcal{X}. \end{equation}

The three dimensionless parameters describing the problem are

(4.6a–c) \begin{align} |\textit{Ro}_{ac}| := \frac {1}{2}\frac {\Delta \kern-1pt V}{f L_x}, \quad \varGamma _m := 2\frac {c_m}{\Delta \kern-1pt V}\frac {c_m}{fL_x} \equiv 2\frac {L_m}{L_x}, \quad \xi := \frac {L_x}{L_f} \gg 1. \end{align}

The tunnelling parameter $\gamma _m$ , filament Rossby number $\textit{Ro}_{\!f}$ and $\alpha _m$ may be expressed in terms of these parameters as

(4.7a–c) \begin{align} \gamma _m := \frac {L_m}{L_f} \equiv \frac {1}{2}\xi \varGamma _m, \quad \textit{Ro}_{f} \equiv 2\xi |\textit{Ro}_{ac}|, \quad \alpha _m^2 := \frac {\textit{Ro}_{\!f}}{\gamma _m} \equiv \frac {\Delta \kern-1pt V^2}{c_m^2} \equiv 4\frac {|\textit{Ro}_{ac}|}{\varGamma _m}. \end{align}

The imposition of periodic boundary conditions allows us to expand $\mathcal{X}$ into discrete normal modes

(4.8a) \begin{equation} \mathcal{X} = \sum _{n=0}^\infty \mathcal{X}_n \mathrm{e}^{\mathrm{i}\tilde {\omega }_n\tilde {t}} + \mathrm{c.c.},\end{equation}

where each mode satisfies

(4.8b) \begin{equation} \left [|\textit{Ro}_{ac}|^{-1}\big (1 - \tilde {\omega }^2_n\big ) - \varGamma _m \frac {\partial ^2}{\partial \hat {x}^2} + \frac {\partial \hat {V}}{\partial \hat {x}}\right ]\mathcal{X}_n. \end{equation}

We have a periodic Sturm–Liouville eigenvalue problem, more precisely the one-dimensional time-independent Schrödinger equation, which we may rearrange into the standard form

(4.9) \begin{equation} \left [- \varGamma _m \frac {\partial ^2\,}{\partial \hat {x}^2} + \frac {\partial \hat {V}}{\partial \hat {x}}\right ]\mathcal{X}_n = \lambda _n\mathcal{X}_n, \end{equation}

where $\lambda _n$ is the ${n}$ th eigenvalue. The problem is mathematically equivalent to the well-studied quantum mechanics problem of a particle in a one-dimensional crystal lattice. The eigenvalue problem is independent of $|\textit{Ro}_{ac}|$ . However, the relationship between the frequency and the eigenvalues

(4.10a) \begin{equation} \tilde {\omega }_n^2 = 1 + |\textit{Ro}_{ac}|\lambda _n \end{equation}

does depend on $|\textit{Ro}_{ac}|$ . If we were to make the YBJ approximation (2.10) then we would arrive at the exact same eigenvalue problem, but the frequency, including the carrier inertial oscillation, would be given by

(4.10b) \begin{equation} \tilde {\omega }_n^{(\textit{YBJ}\kern1pt)} = 1 + \frac {1}{2}|\textit{Ro}_{ac}|\lambda _n.\end{equation}

Furthermore, consider the role of the parameter $\xi$ . It only appears in (4.9) through $\partial \hat {V}/\partial \hat {x}$ . It is a parameter that determines the shape, but not strength (which is determined by $|\textit{Ro}_{ac}|$ ), of the background lateral shear. There are some results that hold for any $\partial \hat {V}/\partial \hat {x}$ . In these cases the only remaining parameter is $\varGamma _m$ .

One such result is a useful expression for the derivative of the eigenvalues with respect to $\varGamma _m$ . Let $\langle \boldsymbol{\cdot }\rangle := ({1}/{2})\int _{-1}^1\mathrm{d}\hat {x}$ be the domain average. Consider :

(4.11) \begin{equation} -\left \langle \mathcal{X}_n\frac {\partial ^2 \mathcal{X}_n}{\partial \hat {x}^2}\right \rangle + \left \langle \mathcal{X}_n\left [- \varGamma _m \frac {\partial ^2}{\partial \hat {x}^2} + \frac {\partial \hat {V}}{\partial \hat {x}}\right ]\frac {\partial \mathcal{X}_n}{\partial \varGamma _m}\right \rangle = \left \langle \lambda _n\mathcal{X}_n\frac {\partial \mathcal{X}_n}{\partial \varGamma _m}\right \rangle + \frac {\partial \lambda _n}{\partial \varGamma _m}\big \langle {\mathcal{X}_n}^2\big \rangle.\end{equation}

The linear operator is self-adjoint, with respect to the inner product defined by the domain average, and so the second term on the left cancels the first term on the right. Manipulating what is left gives

(4.12) \begin{equation} \frac {\partial \lambda _n}{\partial \varGamma _m} = \frac {\big \langle {\mathcal{X}_n'}^2\big \rangle }{\big \langle {\mathcal{X}_n}^2\big \rangle }, \end{equation}

where $\mathcal{X}_n' := \partial \mathcal{X}_n / \partial \hat {x}$ . As a corollary, we note that all the eigenvalues, and hence frequencies, are increasing functions of $\varGamma _m$ .

4.3. The minimum frequency mode

4.3.1. Bounds

The minimum frequency mode, i.e. the eigenmode with the smallest eigenvalue, is the zeroth mode and is special in that the solution has no zero crossings. This allows us to derive bounds on the minimum frequency. Since $\mathcal{X}_0'$ is periodic, $\mathcal{X}_0'' := \partial ^2 \mathcal{X}_0 / \partial \hat {x}^2$ must take both signs somewhere in the domain. It then follows from (4.9) that $\partial \hat {V} / \partial \hat {x} - \lambda _0$ must also take both signs. This implies the lower bound $\lambda _0 \gt \min _{\hat {x}} \partial \hat {V} / \partial \hat {x} = -1$ . To derive an upper bound, divide (4.9) by $\mathcal{X}_0$ and average to get

(4.13) \begin{equation} -\varGamma _m\left \langle \frac {1}{\mathcal{X}_0}\frac {\partial ^2 \mathcal{X}_0}{\partial \hat {x}^2}\right \rangle + \left \langle \frac {\partial \hat {V}}{\partial \hat {x}}\right \rangle = \lambda _0\ \implies \lambda _0 = -\varGamma _m \left \langle \left (\frac {\mathcal{X}_0'}{\mathcal{X}_0}\right )^2\right \rangle + \left \langle \frac {\partial \hat {V}}{\partial \hat {x}}\right \rangle \lt \left \langle \frac {\partial \hat {V}}{\partial \hat {x}}\right \rangle = 0 ,\end{equation}

where the second equality follows from integration by parts and the periodic boundary conditions. Together these bounds are

(4.14) \begin{equation} -1 \lt \lambda _0 \lt 0\ \implies 1 + \textit{Ro}_{ac} \lt \tilde {\omega }_0^2 \lt 1. \end{equation}

These are again examples of results that hold for arbitrary periodic background states. Redimensionalising, we can in general say that the minimum frequency $\omega _0$ must satisfy

(4.15) \begin{equation} \min _x f_{\textit{eff}}^2(x) \lt \omega _0^2 \lt \underset {x}{\mathrm{mean}} \ f_{\textit{eff}}^2(x). \end{equation}

We now consider a couple of cases where we can make analytic progress with a focus on the minimum frequency mode.

4.3.2. The WKBJ regime – $\gamma _m \ll 1$

In the regime $\gamma _m \equiv ({1}/{2})\xi \varGamma _m \ll 1$ we may again utilise a WKBJ approximation. In the interests of brevity, we assume that the filament is symmetric and so we need only solve over $0 \lt \hat {x} \lt 1$ . If we further assume that the filament has a single maximum, like the Gaussian filament, then we have a classic two turning point eigenvalue problem as found in Bender & Orszag (Reference Bender and Orszag1999). We define the slowly varying wavenumber for mode $n$ ,

(4.16) \begin{equation} \hat {k}_n(\hat {x};\lambda _n) := \sqrt {\frac {|-1 + 2\xi \mathcal{F}(\xi \hat {x}) - \lambda _n|}{\varGamma _m}}, \end{equation}

that is determined by the eigenvalue and the local (non-dimensional) vorticity of the background flow. Enforcing even boundary conditions the solution is

(4.17) \begin{equation} \mathcal{X}_{2n} = \begin{cases} C_1 \hat {k}_{2n}^{-1/2}\cosh \left (\int _0^{\hat {x}} \hat {k}_{2n}(\hat {x}')\mathrm{d}\hat {x}'\right )\!, & 0 \leqslant \hat {x} \lt \hat {x}_* \\ C_2 \hat {k}_{2n}^{-1/2}\cos \left (\int _{\hat {x}}^1 \hat {k}_{2n}(\hat {x}')\mathrm{d}\hat {x}'\right )\!, & \hat {x}_* \lt \hat {x} \leqslant 1, \end{cases} \end{equation}

where $\hat {x}_*$ is the turning point defined by $\hat {k}_{2n}(\hat {x}_*;\lambda _{2n}) = 0$ assuming one exists. The eigenvalues are determined by the connection formula across the turning point,

(4.18a) \begin{align} \frac {1}{2}\mathrm{e}^{-2A_{2n}} = \tan \left (\frac {\pi }{4} - B_{2n}\right )\!, \end{align}

where

(4.18b,c) \begin{align} A_{2n} := \int _0^{\hat {x}_*}\hat {k}_{2n}(x')\mathrm{d}x', \quad B_{2n} := \int _{\hat {x}_*}^1\hat {k}_n(x')\,\mathrm{d}x'. \end{align}

A brief derivation of these formulae is given in Appendix D. We also show that $B_{2n}$ is bounded. In particular, $B_{2n} \lt n\pi + \pi /4$ . Rewriting (4.18b , c ), we have

(4.19) \begin{equation} \sqrt {\varGamma _m}B_{2n} = \int _{\hat {x}_*}^1\sqrt {\lambda _{2n} + 1 - 2\xi \mathcal{F}(\xi x')}\,\mathrm{d}x'. \end{equation}

As $\varGamma _m \to 0$ , this integral must vanish and the only way this can occur is if $\lambda _{2n} \to -1$ . Therefore, all the modes, including the zeroth mode, will, for vanishingly small $\varGamma _m$ , achieve the lower bound for the minimum frequency derived earlier. For analytic filament structure functions $\mathcal{F}$ , including the Gaussian filament, this also requires $\hat {x}_* \to 1$ and the region in which the solution is wavelike becomes increasingly localised. Furthermore, in this limit the $A_{2n}$ integral is very large and the zeroth mode has $B_0 \to \pi /4$ and in the anticyclonic region $\mathcal{X}_0(\hat {x}) = \cos (\pi /4 - \int _{\hat {x}_*}^{\hat {x}} \hat {k}\mathrm{d}x')$ . One can check that the low horizontal modes become exponentially small in the centre of the filament (D7, D10 and figure 8 d). We also note that, as $\varGamma _m \to 0$ , the leading-order correction to the minimum frequency is $O(|\textit{Ro}_{ac}|)$ (figure 8 a) and that there will be $O(|\textit{Ro}_{ac}|)$ differences between the YBJ approximation to the frequency and the exact value (figure 8 b).

Figure 8. Minimum frequency eigenvalues and eigenmodes for a Gaussian vorticity filament. (a) Eigenvalue as a function of $\varGamma _m$ for different $\xi$ including theoretical results as $\xi \to \infty$ . The right axis is the frequency squared (4.10a ) for $|\textit{Ro}_{ac}| = 0.5$ . Dashed grey lines indicate the upper and lower bounds (4.14). (b) Minimum frequency computed using the Klein–Gordon equation (solid line, 4.10a ) and YBJ approximation (dashed line, 4.10b ) for $\xi = 10$ and $|\textit{Ro}_{ac}| = 0.5$ . (c) Large $\varGamma _m$ behaviour of the eigenvalues. Using the delta-function solution $\lambda _0\varGamma _m$ tends to a constant $-1/3$ (grey dashed line). (d) Structure of the lowest frequency modes for $\xi = 10$ and four values of $\varGamma _m$ .

Figure 9. Projection of the uniform initial condition $\mathcal{X} = 1$ onto even horizontal modes for $\xi = 10$ and various $\varGamma _m$ . The first five horizontal modes and their sum are plotted for points in (a) the WKBJ regime, $\varGamma _m = 0.005 \implies \gamma _m = 0.1$ , and (b) the tunnelling regime, $\varGamma _m = 0.5\ \implies \gamma _m = 10$ . (c) Energy content of the even horizontal modes and (d) frequency squared for $|\textit{Ro}_{ac}| = 0.5$ for four values of $\varGamma _m$ .

In this regime we have made no assumption about the value of $\xi$ , indeed the WKBJ approach works for arbitrary $\partial \hat {V}/\partial \hat {x}$ . We only require $\gamma _m \ll 1$ . However, since $\xi \gg 1$ , in this regime we must have $\varGamma _m \lll 1$ . Unless we also have $|\textit{Ro}_{ac}| \lll 1$ , $c_m^2/f^2L_x^2 \equiv |\textit{Ro}_{ac}|\varGamma _m \ll 1$ and, thus, we are in the regime discussed in § 4.1 in which $L_x$ is large compared with both $R_m$ and $L_m$ and the boundary conditions are not important. This is reflected in the fact that the eigenvalues are densely packed near to the minimum and that the projection onto horizontal modes is a very inefficient method for representing the solution (figures 9 a and 9 d).

4.3.3. The strong dispersion regime – $\varGamma _m \gg 1$

For large $\varGamma _m$ , the spatial derivatives in (4.9) – the dispersive term – are multiplied by a large parameter that acts to suppress spatial variations in the zeroth mode. For the higher modes, the large contribution from the spatial derivatives can be balanced by a large eigenvalue. However, the eigenvalue of the zeroth mode is bounded (4.14) and as $\varGamma _m$ gets larger the zeroth mode becomes more uniform. Young & Ben Jelloul (Reference Young and Jelloul1997), and later Conn et al. (Reference Conn, Callies and Lawrence2025), derive a result, known as the ‘strong dispersion approximation’, that carries over to the zeroth mode of the Klein-Gordon equation for arbitrary $\partial \hat {V}/\partial \hat {x}$ in this regime. The idea is to look for a small correction to the uniform state, i.e. the first-order term in an expansion in $\varGamma _m^{-1}$ .

Let $\mathcal{X}_0 = 1 +\varGamma _m^{-1}\mathcal{X}_0^{(1)} + O(\varGamma _m^{-2})$ and $\lambda _0 = 0 +\varGamma _m^{-1}\lambda _0^{(1)} + O(\varGamma _m^{-2})$ . At $O(1)$ , (4.9) gives

(4.20) \begin{equation} -\frac {\partial ^2 \mathcal{X}_0^{(1)}}{\partial \hat {x}^2} + \frac {\partial \hat {V}}{\partial \hat {x}}\boldsymbol{\cdot }1 = 0. \end{equation}

Therefore,

(4.21a,b) \begin{align} \frac {\partial \mathcal{X}_0}{\partial \hat {x}} = \varGamma _m^{-1}\hat {V} + O\big(\varGamma _m^{-1}\big), \quad \mathcal{X}_0 = 1 - \varGamma _m^{-1}\int ^{\hat {x}} \hat {V} \mathrm{d}\hat {x}' + O\big(\varGamma _m^{-2}\big). \end{align}

The leading-order correction to the eigenvalue may be found by considering the secularity condition at $O(\varGamma _m^{-1})$ , but also follows from (4.11),

(4.22a) \begin{equation} \frac {\partial \lambda _0}{\partial \varGamma _m^{-1}} = -\varGamma _m^2 \frac {\partial \lambda _0}{\partial \varGamma _m} = -\varGamma _m^2\frac {\big \langle {\mathcal{X}_0'}^2\big \rangle }{\big \langle {\mathcal{X}_0}^2\big \rangle } = - \langle \hat {V}^2 \rangle + O\big(\varGamma _m^{-1}\big). \end{equation}

Thus,

(4.22b) \begin{equation} \lambda _0 = -\langle \hat {V}^2 \rangle \varGamma _m^{-1} + O\big(\varGamma _m^{-2}\big). \end{equation}

Unwrapping the non-dimensionalisation, the leading-order correction to the frequency is given by

(4.23) \begin{equation} \omega _0^2 = f^2\big (1 + c_m^{-2} \langle \hat {V}^2 \rangle \big ). \end{equation}

We have recovered Young & Ben Jelloul (Reference Young and Jelloul1997)’s result that the leading-order correction is proportional to the average kinetic energy of the background flow.

The leading-order correction to the minimum frequency is $O(\Delta \kern-1pt V^2/c_m^2 \equiv \alpha _m^2)$ . Similar to the $\gamma _m \to \infty$ limit in § 3, we find that in the $\varGamma _m \to \infty$ limit the parameter determining the leading-order behaviour is $\alpha _m$ , which does not depend on the length scales of the background flow. However, unlike in § 3 we cannot look at the distinguished limit $\varGamma _m \to \infty$ with $\alpha _m$ held fixed as this would require $|\textit{Ro}_{ac}| \to \infty$ . This is unphysical as it results in an inertially unstable background flow (for $|\textit{Ro}_{ac}| \gt 1$ , the minimum frequency squared can be negative (4.10a )). This limitation means that the leading-order correction to the minimum frequency will always be small in the strong dispersion regime and, hence, the YBJ approximation will always be good.

4.4. Delta-function filament

So far we have considered the WKBJ regime, $\gamma _m \equiv ({1}/{2})\xi \varGamma _{m} \ll 1$ , and the minimum frequency mode in the strong dispersion regime, $\varGamma _m \gg 1$ . By construction $\xi \gg 1$ , but otherwise these regimes place no restriction on the value of $\xi$ . Indeed these arguments hold with very few restrictions on $\partial \hat {V}/\partial \hat {x}$ . We now consider a delta-function filament, which requires $L_f \ll L_m$ ( $\gamma _m \gg 1$ ) and $L_f \ll L_x$ ( $\xi \gg 1$ ) but places no restriction on the ratio of $L_m$ to $L_x$ . The delta-function filament is therefore reached through the distinguished limit $\xi \to \infty$ with $\varGamma _m$ held fixed. The WKBJ, strong dispersion and delta-function regimes are summarised in table 5.

Table 5. Summary of the different regimes considered in the periodic problem.

Away from the filament the even modes have the form

(4.24a) \begin{equation} \mathcal{X}_{2n} = \cos {K_{2n}(1 - \hat {x})} \end{equation}

with

(4.24b) \begin{equation} \lambda _{2n} = -1 + \varGamma _m K_{2n}^2. \end{equation}

The jump condition is $\varGamma _m\mathcal{X}_n' = \mathcal{X}_n$ at $\hat {x} = 0^+$ . Therefore, $K_{2n}$ satisfies

(4.24c) \begin{equation} K_{2n}\tan {K_{2n}} = \varGamma _m^{-1}, \end{equation}

where $K_{2n}$ is the unique solution in $(n\pi ,(n+1/2)\pi )$ . For any given filament, the delta-function approximation will fail for sufficiently large $n$ when the length scale of the waves becomes too short.

We now consider the limiting behaviours as $\varGamma _m \to 0$ and $\varGamma _m \to \infty$ . For $\varGamma _m \ll 1$ , let $K_{2n} = (n + ( {1}/{2}))\pi - \varGamma _m k_{2n} + O(\varGamma _m^2)$ . With this expansion, $K_{2n}\tan {K_{2n}} = (n + ({1}/{2}))\pi / (\varGamma _m k_{2n}) + O(1)$ . Therefore, as $\varGamma _m \to 0$ ,

(4.25a,b) \begin{align} K_{2n} = \left (n \!+\! \frac {1}{2}\right )\pi - \varGamma _m \frac {1}{\left (n \!+\! \frac {1}{2}\right )\pi } \!+\! O\big (\varGamma _m^2\big ), \,\, \lambda _{2n} = -1 \!+\! \varGamma _m \left (n + \frac {1}{2}\right )^2\pi ^2 + O\big (\varGamma _m^2\big ). \end{align}

In particular, $\lambda _0 \to -1$ as $\varGamma _m \to 0$ and, thus, even in the delta-function limit we can achieve the lower bound on the minimum frequency (figure 8 a). However, to be in the delta-function limit $\gamma _m \gg 1$ when $\varGamma _m \ll 1$ we must have $\xi \ggg 1$ . Furthermore, there are differences in the structure of the eigenmodes in the WKBJ and delta-function regimes as $\varGamma _m \to 0$ . For example, consider the phase of the zeroth mode as it enters the filament. In the delta-function solution, as $\hat {x} \to 0$ , $\mathcal{X}_0 \to \cos {\pi / 2} = 0$ . Whereas in the WKBJ solution, at $\hat {x} = \hat {x}_*$ , $\mathcal{X}_0 \sim \cos B_0 = \cos \pi /4$ , which represents a phase shift of $\pi /4$ .

For large $\varGamma _m$ , we need different expansions for $n = 0$ and $n \geqslant 1$ , highlighting that the zeroth mode is special. First, $n \geqslant 1$ , where we let $K_{2n} = n\pi + \varGamma _m^{-1}k_{2n} + O(\varGamma _m^{-2})$ . Working through the expansion, we find that, for $\varGamma _m \gg 1$ , $n \geqslant 1$ ,

(4.26a,b) \begin{align} K_{2n} = n\pi + \varGamma _m^{-1} \frac {1}{n\pi } + O\big(\varGamma _m^{-2}\big), \quad \lambda _{2n} = \varGamma _m n^2\pi ^2 + 1 + O\big(\varGamma _m^{-1}\big).\end{align}

Note that the eigenvalues are $O(\varGamma _m)$ and thus large. To leading order, the frequency is given by

(4.26c) \begin{align} \omega _{2n}^2 = n^2\pi ^2|\textit{Ro}_{ac}|\varGamma _m = n^2\pi ^2 \frac {c_m ^2}{f^2L_x^2}, \end{align}

which is independent of the properties of the background flow with the exception of the length scale $L_x$ . Here $L_x$ is setting the length scale of the waves $L_w = L_x/n\pi$ and the frequency is determined by the wave Burger number $(c_m/fL_w)^2$ .

For, $n = 0$ , we expand $K_0^2$ in powers of $\varGamma _m^{-1}$ , $K_0^2 = \varGamma _m^{-1} a_1 + \varGamma _m^{-2} a_2 + O(\varGamma _m^{-3})$ . We have $K_0 \tan {K_0} = K_0^2 + K_0^4/3 + O(K_0^6) = \varGamma _m^{-1}a_1 + \varGamma _m^{-2}(a_2 + a_1^2/3) + O(\varGamma _m^{-3})$ . Therefore, for $\varGamma _m \gg 1$ ,

(4.27a,b) \begin{align} K_{0}^2 = \varGamma _m^{-1} -\frac {1}{3}\varGamma _m^{-2} + O\big(\varGamma _m^{-3}\big), \quad \lambda _{0} = -\frac {1}{3}\varGamma _m^{-1} + O\big(\varGamma _m^{-2}\big). \end{align}

This not only achieves the upper bound $\lambda _0 \to 0$ as $\varGamma _m \to \infty$ (figure 8 a), but we observe that the wavenumber $K_0 \to 0$ as well. The solution is tending towards a uniform state. Moreover, the coefficient $- {1}/{3}$ in (4.27a , b ) is consistent with the strong dispersion approximation, i.e. for the sawtooth geostrophic velocity $\hat {V} = \mathrm{sign}(\hat {x})(1 - \hat {x})$ associated with the delta-function filament, we have $\langle \hat {V}^2\rangle = {1}/{3}$ .

The utility of the delta-function filament is not only that it provides an analytical result for intermediate values of $\varGamma _m$ but also that it is a good model for sharp submesoscale filaments. Solving the eigenvalue problem numerically for the Gaussian filament, we observe very good agreement between the eigenvalues computed for $\xi = 10, 20$ and the theoretical values assuming a delta-function (figure 8 a). This is important as these values are realisable in submesoscale flows, for example, taking $|\textit{Ro}_{ac}| = 0.5$ and $\xi = 10$ , the peak vorticity is only $\partial V/\partial x = 3.5 f$ , and thus, the delta-function filament is a relevant model for oceanic applications.

4.5. Lateral refraction

In this problem we are considering a discrete spectrum of horizontal modes. Lateral refraction and the associated horizontal flux of NIW energy are determined by the time scale over which the different modes dephase, the structure of each mode and their relative energy content. To understand this process, we look at how the laterally uniform initial condition $\mathcal{X} = 1$ projects onto horizontal modes. Conveniently, the Sturm–Liouville problem has a trivial weight function and, thus, the projection onto horizontal modes satisfies a particularly simple form of Parseval’s theorem,

(4.28) \begin{equation} \langle \mathcal{X}^2 \rangle = 1 = \sum _n \big \langle \mathcal{X}_n^2 \big \rangle , \end{equation}

where $\langle \boldsymbol{\cdot }\rangle$ again denotes a lateral average. For a given vertical mode, $\mathcal{X}$ is proportional to the across-filament velocity and we may interpret $\langle \mathcal{X}_n^2 \rangle$ as the fraction of the across-filament kinetic energy contained in the $n$ th horizontal mode. In the strong dispersion regime, $\varGamma _m \gg 1$ , the initial condition projects almost entirely onto the zeroth mode. The kinetic energy in the higher modes drops off exponentially (figures 9 b and 9 c) and, hence, the solution (e.g. figure 7 h) remains nearly uniform even as the different modes dephase. Even for smaller $\varGamma _m$ , the zeroth mode dominates the energy content. For example, in the case $\varGamma _m = 0.005$ , $\xi = 10$ , more than 65 % of the energy is in the zeroth mode. However, many horizontal modes are required to reproduce the initial condition (figures 9 a and 9 c).

The time scale over which the different horizontal modes dephase is determined by the differences in their frequencies (figure 9 d). While these differences are strongly dependent on $\varGamma _m$ and $|\textit{Ro}_{ac}|$ (differences in $\tilde {\omega }_n^2$ scale linearly with $|\textit{Ro}_{ac}|$ ), the refraction time scales can be of the order of an inertial period (e.g. figure 7 b). This is much faster than the refraction time scales at small Rossby numbers (e.g. Balmforth et al. Reference Balmforth, Llewellyn and Young1998; Danioux et al. Reference Danioux, Vanneste and Bühler2015; Rocha, Wagner & Young Reference Rocha, Wagner and Young2018; Asselin et al. Reference Asselin, Thomas, Young and Rainville2020), as one might expect given the large Rossby numbers in the filament we are considering. This is an important result because in the submesoscale regime the background flows also evolve on inertial time scales. However, we find that $\zeta$ -refraction can occur even faster. The rapid lateral refraction and their domination of the energy content means that the zeroth modes dictate the long time dynamics of the NIWs. With this in mind we now consider how the modes propagate in the vertical.

5.1. Group velocity

5.1.1. General expressions

With these choices we can consider the vertical group velocity of a horizontal mode:

(5.1) \begin{equation} c_{g,n} := \frac {\partial \omega _n}{\partial k_z}. \end{equation}

However, we can use $\varGamma _m \equiv c_m^2 / (f^2L_x^2|\textit{Ro}_{ac}|)$ , (4.10a ) and (4.11) to express the group velocity as

(5.2) \begin{equation} c_{g,n} \equiv \frac {\partial \varGamma _m}{\partial k_z}\frac {\partial \omega _n}{\partial \varGamma _m} \equiv \frac {-2N^2|\textit{Ro}_{ac}|^{-1}}{k_z^3L_x^2f^2}\frac {1}{2\omega _n}\frac {\partial \omega _n^2}{\partial \varGamma _m} \equiv \frac {-N^2}{\omega _nk_z^3}\frac {1}{L_x^2}\frac {\partial \lambda _n}{\partial \varGamma _m} \equiv -\frac {N^2}{\omega _nk_z^3}\frac {1}{L_x^2}\frac {\big \langle {\mathcal{X}_n^\prime }^2\big \rangle }{\big \langle {\mathcal{X}_n}^2\big \rangle }. \end{equation}

This final expression is reminiscent of the ray-tracing expression for the group velocity and would be the same if we replaced $\langle {\mathcal{X}_n^\prime }^2\rangle /\langle {\mathcal{X}_n}^2\rangle$ by $(k_xL_x)^2$ , where $k_x$ is the horizontal wavenumber. Furthermore, this expression for the group velocity holds for arbitrary $\partial \hat {V}/\partial \hat {x}$ . One of the key motivations for studying $\zeta$ -refraction is the intuition that it is necessary to shrink the horizontal scales of NIWs in order to allow them to propagate rapidly (Gill Reference Gill1984; Balmforth et al. Reference Balmforth, Llewellyn and Young1998). This intuition is based upon the plane wave dispersion relation and group velocity but we can see that it holds more generally: the generation of horizontal gradients is critical to the vertical propagation of NIWs.

This expression also holds under the YBJ approximation if we replace $\omega _n$ by $f$ ,

(5.3) \begin{equation} c_{g,n}^{(\textit{YBJ}\kern1pt)} \equiv -\frac {N^2}{f k_z^3}\frac {1}{L_x^2}\frac {\big \langle {\mathcal{X}_n^\prime }^2\big \rangle }{\big \langle {\mathcal{X}_n}^2\big \rangle }. \end{equation}

As a result, YBJ theory always underpredicts the group velocity of the minimum frequency mode, since $\omega _0 \lt f$ , but overpredicts the group velocity of the super-inertial high modes.

Here we ask: For a given background flow, what wavenumber maximises the magnitude of the vertical group velocity? It is natural to non-dimensionalise and express the group velocity in terms of $\varGamma _m$ , $|\textit{Ro}_{ac}|$ and $\xi$ . For given $|\textit{Ro}_{ac}|$ and $\xi$ , we then find the optimal $\varGamma _m$ . The non-dimensional group velocity is given by

(5.4) \begin{equation} \hat {c}_{g,n} := \frac {N}{f}\frac {c_{g,n}}{fL_x} = -\frac {c_m^3}{f^3L_x^3}\frac {1}{\tilde {\omega }_n}\frac {\big \langle {\mathcal{X}_n^\prime }^2\big \rangle }{\big \langle {\mathcal{X}_n}^2\big \rangle } = -|\textit{Ro}_{ac}|^{3/2}\tilde {\omega }_n^{-1}\varGamma _m^{3/2}\frac {\big \langle {\mathcal{X}_n^\prime }^2\big \rangle }{\big \langle {\mathcal{X}_n}^2\big \rangle }.\end{equation}

In this expression, there is explicit dependence on both $|\textit{Ro}_{ac}|$ and $\varGamma _m$ . The contribution from the horizontal structure

(5.5) \begin{equation} \mathcal{S}_n := \frac {\big \langle {\mathcal{X}_n^\prime }^2\big \rangle }{\big \langle {\mathcal{X}_n}^2\big \rangle } \end{equation}

depends only on $\varGamma _m$ and $\xi$ , but not $|\textit{Ro}_{ac}|$ . The frequency term depends on all three dimensionless parameters. However, under the YBJ approximation the frequency dependence drops out,

(5.6) \begin{equation} \hat {c}_{g,n}^{(\textit{YBJ}\kern1pt)} := -\frac {N}{f}\frac {c_{g,n}^{(\textit{YBJ}\kern1pt)}}{fL_x} = -|\textit{Ro}_{ac}|^{3/2}\varGamma _m^{3/2}\frac {\big \langle {\mathcal{X}_n^\prime }^2\big \rangle }{\big \langle {\mathcal{X}_n}^2\big \rangle }. \end{equation}

Ignoring the frequency term for the moment, we see that $\varGamma _m$ has two roles in determining the group velocity. First, there is the explicit dependence on $\varGamma _m^{3/2}$ and second, it determines the structural contribution $\mathcal{S}_n$ .

5.1.2. Delta-function filament

For the delta-function filament, we can compute $\mathcal{S}_{2n}$ exactly using (4.24a ),

(5.7) \begin{equation} \mathcal{S}_{2n} := \frac {\big \langle {\mathcal{X}_{2n}^\prime }^2\big \rangle }{\left \langle {\mathcal{X}_{2n}} ^2\right \rangle } = K_{2n}^2\frac {\int _0^1\sin ^2{K_{2n}(1 - \hat {x})}\,\mathrm{d}\hat {x}}{\int _0^1\cos ^2{K_{2n}(1 - \hat {x})}\,\mathrm{d}\hat {x}} = K_{2n}^2\frac {1 - \frac {1}{2K_{2n}}\sin {2K_{2n}}}{1 + \frac {1}{2K_{2n}}\sin {2K_{2n}}}. \end{equation}

Furthermore, we can use the asymptotic expansions for $K_{2n}$ and $\tilde {\omega }_n$ , in the limit of small and large $\varGamma _m$ , derived in § 4.4 to compute $\mathcal{S}_{2n}$ and $\hat {c}_{g,2n}$ in these limits.

However, in all but one case the leading-order term in the expansion for $K_{2n}$ is a constant. Consequently, $\mathcal{S}_{2n}$ also tends to a constant. For example, as $\varGamma _m \to 0$ , $K_0 \to \pi /2$ and, thus, $\mathcal{S}_0 \to \pi ^2/4$ . In general, as $\varGamma _m \to 0$ , $\mathcal{S}_{2n} \to (n + ({1}/{2}))^2\pi ^2$ . Moreover, $\tilde {\omega }_n \to \sqrt {1 - |\textit{Ro}_{ac}|}$ and, hence, $\hat {c}_{g,n} \sim \varGamma _m^{3/2}$ (figure 10 a). In the limit of small $\varGamma _m$ we have recovered the familiar $k_z^{-3}$ vertical wavenumber dependence of ray-tracing.

Figure 10. (a) Non-dimensionalised group velocity of the first five even modes for the Gaussian filament with $\textit{Ro}_{ac} = -0.5$ and $\xi = 10$ . (b) Horizontal structure contribution $\mathcal{S}_0$ to the group velocity of the zeroth mode for $\xi = 10, 25$ with asymptotic scalings for the delta-function filament in red and blue. (c) Non-dimensionalised group velocity of the zeroth mode for $\xi = 10$ and four values of $\textit{Ro}_{ac}$ . Both the Klein–Gordon (5.2, black) group velocity and the YBJ approximation (5.3, blue) are plotted. The maxima of the group velocities are denoted by stars. (d) Optimal value of $\varGamma _m$ , $\varGamma _m^\star$ , for radiating NIW energy as a function of $\xi ^{-1}$ for four values of $\textit{Ro}_{ac}$ . Under the YBJ approximation the optimal value is independent of $\textit{Ro}_{ac}$ . Green lines in (b) and (d) denote $\varGamma _m = 0.956$ , the location of the maximum of $\varGamma _m^{3/2}\mathcal{S}_0$ for the delta-function filament.

The limit $\varGamma _m \to \infty$ is more interesting. For $n \geqslant 1$ (4.26a ), $K_{2n} \to n\pi$ and, thus, $\mathcal{S}_{2n} \to n^2\pi ^2$ . However, in this case the leading-order contribution to the frequency also depends on $\varGamma _m$ , i.e.

(5.8a) \begin{equation} \tilde {\omega }_n \to n\pi |\textit{Ro}_{ac}|^{1/2} \varGamma _m ^ {1/2} \end{equation}

and, thus, the group velocity is only linear in $\varGamma _m \sim k_z^{-2}$ (figure 10 a),

(5.8b) \begin{equation} \hat {c}_{g,n} \to - n\pi |\textit{Ro}_{ac}|\varGamma _m \equiv -n\pi \frac {c_m^2}{f^2L_x^2}. \end{equation}

Finally, consider the zeroth mode in the strong dispersion regime $\varGamma _m \to \infty$ . In this limit, $K_0 \to \varGamma _m^{-1/2}\ \implies \mathcal{S}_0 \to (1/3)K_0^4 \to (1/3)\varGamma _m^{-2}$ (figure 10 b). Since $\tilde {\omega }_0 \to 1$ , we have

(5.9) \begin{equation} \hat {c}_{g,0} \to -\frac {1}{3}|\textit{Ro}_{ac}|^{3/2}\varGamma _m^{-1/2}. \end{equation}

In the strong dispersion regime the contribution from the horizontal structure is very sensitive to the vertical wavenumber to the extent that the group velocity is proportional to $k_z$ .

This result for the zeroth mode in the strong dispersion regime generalises to arbitrary periodic background flows through the strong dispersion approximation. From (4.22a ) we have $\mathcal{S}_0 \to \varGamma _m^2\langle \hat {V}^2\rangle$ and, thus,

(5.10a) \begin{equation} \hat {c}_{g,0} \to -|\textit{Ro}_{ac}|^{3/2}\langle \hat {V}^2 \rangle \varGamma _m^{-1/2}. \end{equation}

Redimensionalising,

(5.10b) \begin{equation} c_{g,0} \to -\frac {f\langle \hat {V}^2 \rangle }{Nc_m} \equiv -\frac {f\langle \hat {V}^2 \rangle }{N^2}k_z ,\end{equation}

consistent with Young & Ben Jelloul (Reference Young and Jelloul1997). We again emphasise that under the strong dispersion approximation the dynamics are independent of the length scales of the background flow and that results derived using the YBJ approximation extend to the Klein–Gordon equation as the zeroth mode frequency tends to $f$ .

5.1.3. Optimal radiation

The special behaviour of the zeroth mode, i.e. that the group velocity magnitude becomes a decreasing function of $\varGamma _m$ for large $\varGamma _m$ , means that there exists an optimal value of $\varGamma _m$ for radiating energy, generalising the findings of Balmforth et al. (Reference Balmforth, Llewellyn and Young1998). We call this value $\varGamma _m^\star$ . Thus, for a given background flow, i.e. fixed values of $\Delta \kern-1pt V$ , $L_f$ , $L_x$ , $N$ and $f$ , the optimal vertical wavenumber for radiating NIW energy is

(5.11) \begin{equation} k_z^\star := \sqrt {\frac {2 N^2}{\Delta \kern-1pt V f L_x}\frac {1}{\varGamma _m^\star }}. \end{equation}

Under the YBJ approximation, the group velocity (5.3) factorises into a term depending on $|\textit{Ro}_{ac}|$ and a term depending on $\varGamma _m$ (and $\xi$ ) but not $|\textit{Ro}_{ac}|$ . As a result, $\varGamma _m^\star$ is independent of $|\textit{Ro}_{ac}|$ (figures 10 c and 10 d) and occurs at the maximum of $\varGamma _m^{3/2}\mathcal{S}_0$ . For the delta-function filament, this maximum occurs at $\varGamma _m^\star = 0.956$ with value $\varGamma _m^{3/2}\mathcal{S}_0 = 0.20$ (figure 10 b).

Using the Klein-Gordon group velocity expression with frequency dependence (5.2) not only increases the group velocity of the zeroth mode but also shifts the location of the maximum to smaller $\varGamma _m$ (figures 10 c and 10 d). Furthermore, at finite $\xi$ , $\varGamma _m^\star$ is smaller than the delta-function value (figure 10 d). For the Gaussian filament with $\xi = 10$ and $\textit{Ro}_{ac} = -1$ , $\varGamma _m^\star$ is approximately half the delta-function value of $0.956$ .

The non-monotonic $\varGamma _m$ dependence of the zeroth mode group velocity and the subsequent existence of an optimal $\varGamma _m$ for vertically radiating NIW energy mirrors the $\gamma _m$ dependence of the decay time scale of the lateral radiation problem in § 3. Both problems highlight the importance of $O(1)$ parameter regimes that are not easily tackled by asymptotic approaches. Notably, this non-monotonic behaviour only occurs for the zeroth horizontal mode. However, as we argued in § 4.5, for uniform initial conditions, this mode dominates the energetics even for small $\varGamma _m$ . Moreover, the higher horizontal modes have larger group velocities since $\mathcal{S}_n$ is larger and so they propagate away more quickly leaving the zeroth mode behind. We now demonstrate this explicitly with 2-D linear simulations.

5.2. Linear simulations

We run 2-D linear simulations studying the evolution of a slab initial condition in the across-filament velocity,

(5.12) \begin{equation} u_i := \frac {1}{2}\left (1 + \mathrm{erf}\left (3 + \frac {2z}{H}\right ) \! \right )\!. \end{equation}

We project $u_i$ onto cosine modes and evolve each vertical mode, with the exception of the barotropic mode that we discard, according to the Klein–Gordon equation. For the background flow, we use the Gaussian filament with $\textit{Ro}_{ac} = -0.5$ and $\xi = 10$ . We take the depth of the domain to be $L_z = 0.1L_x$ and $H = 2\times 10^{-3}L_x$ . We use $512$ vertical modes and $512$ horizontal grid points. With $u$ expressed in cosine modes, it is simple to compute $w$ from continuity and then $v$ and $b$ from (2.1b , d ).

In the initial condition the energy content of each vertical mode is a decreasing function of $m$ , but the shear is maximal at $k_z \approx 1/H$ (figure 11 a). We define bulk parameters $\varGamma _H$ and $\gamma _H$ , analogous to $\varGamma _m$ and $\gamma _m$ , using the shear scale $H$ by

(5.13a,b) \begin{align} \varGamma _H := \frac {1}{|\textit{Ro}_{ac}|}\frac {N^2H^2}{f^2L_x^2}, \quad \gamma _H := \frac {1}{2}\xi \varGamma _H. \end{align}

These bulk parameters best describe the shear containing modes but most of the energy is contained in lower modes with $\varGamma _m \gt \varGamma _H$ and $\gamma _m \gt \gamma _H$ (figure 11 a). Having fixed the value of $H/L_x$ and $|\textit{Ro}_{ac}|$ , we vary $\varGamma _H$ by changing the strength of the stratification, i.e. $N/f$ . We focus on three cases, $N/f = 100$ , $250$ and $500$ , which respectively give $\varGamma _H = 0.08$ , $0.5$ and $2$ (table 6). These particular choices of $N/f$ are inspired by the conditions in the northern Gulf of Mexico and are thus towards the large end of the range of realistic oceanic values, at least away from the low latitudes. However, the ratio $N/f$ has a limited independent value outside of the parameters $\varGamma _H$ and $\gamma _H$ . We could achieve the same $\varGamma _H$ and $\gamma _H$ values for different stratifications by adjusting the value of $H/L_x$ . We discuss the range of realistic values of $\varGamma _H$ and $\gamma _H$ in § 6. That being said, $N/f$ does have some independent value in that it is used to non-dimensionalise the group velocity (5.4). Consequently, if we consider two cases with the same $\varGamma _H$ but different stratifications then the NIW energy will propagate more rapidly in the less stratified case.

Table 6. Two-dimensional linear simulation parameters.

Figure 11. (a) Power spectrum of the initial across-filament velocity (5.12) and vertical shear as a function of the vertical mode number. The vertical grey line indicates $k_z = 1/H$ . (b) Theoretical prediction of the group velocity (5.2) of the zeroth horizontal mode as a function of the vertical mode number for the three different values of $N/f$ simulated. (c–e) Across-filament vertical shear at the centre of the filament. In (c–e) the $x$ axes are in inertial periods and the $y$ axes are $\hat {z} := z/L_x$ . White lines indicate a ray moving at the theoretical maximum group velocity, i.e. the maxima of (b).

An animation of the evolution of $u$ in the three simulations (see supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.10637) demonstrates the two major effects of varying the stratification. Firstly, in the $N/f = 100$ case we observe beams radiating down and out from the centre of the filament. This classical behaviour of inertia-gravity waves is indicative of high frequency waves radiating along rays and occurs because at low $\gamma _m$ , where the zeroth horizontal mode is strongly evanescent in the filament, the initial condition projects broadly onto horizontal modes. After the beams have propagated away, we are left with the slower propagating low horizontal modes – primarily the zeroth mode. The beams are not observed in the higher stratification cases where the projection onto higher horizontal modes is much weaker, although a weak signal from the other low horizontal modes, mostly mode 2, can be observed. The stronger tunnelling in the higher stratification cases is also evidenced by the greater shear amplitude in the slowly propagating wavepacket at the centre of the filament (figure 11 c–e). A higher frequency wavepacket associated with the second horizontal mode is also present in figures 11(d) and 11(e).

The second observed effect is the slower propagation of the zeroth mode as the stratification increases. Plotting the group velocity (5.2) as a function of vertical mode number, which is proportional to $k_z$ , we see that the location of the maximum group velocity increases with $N/f$ since the dispersivity is proportional to $(N/fk_z)^2$ (figure 11 b). Furthermore, the magnitude of the maximum group velocity decreases, which follows from the explicit dependence of the group velocity on $N/f$ . In other words, increasing $N/f$ is equivalent to stretching the group velocity–mode number curve along the $x$ axis while compressing it down the $y$ axis. The result is that, for the energy containing low vertical modes, the magnitude of the group velocity decreases with increasing $N/f$ . This counter-intuitive dependence on stratification in the strong dispersion regime (small $\varGamma _m$ ), first noted by Balmforth et al. (Reference Balmforth, Llewellyn and Young1998), is the opposite of the behaviour one would predict by naively applying ray-tracing. If the energy were contained in the high modes then we could apply ray-tracing and we would observe the opposite behaviour – at high vertical wavenumbers the group velocity is an increasing function of $N/f$ (figure 11 b). In addition, the slices of the across-filament shear at the filament centre (figure 11 c–e) demonstrate vertically coherent wavepackets radiating at the maximum group velocity. The maximum is a stationary point of the group velocity and, hence, wavepackets propagating at the maximum group velocity are weakly dispersive and able to maintain a coherent spatial structure.