2.1. Decomposition into horizontal modes

We begin by non-dimensionalising the YBJ equation. Given the scalings $x,y\sim L$, $\psi \sim \varPsi$ and $t\sim L^2/\varPsi$, we obtain the following non-dimensional form of the YBJ equation:

(2.1)\begin{equation} \frac{\partial{\phi}}{\partial{t}} + \operatorname{\textit{J}}(\psi,\phi) + \frac{{\rm i}\zeta}{2}\phi - \frac{{\rm i}\varepsilon^2}{2}\nabla^2\phi = 0. \end{equation}

Here $\varepsilon ^2 = f\lambda ^2/\varPsi$ is the wave dispersiveness (assuming $f>0$). For readers familiar with Young & Ben Jelloul (Reference Young and Ben Jelloul1997), our $\varepsilon ^2$ is equivalent to their $\varUpsilon ^{-1}$. We remind the reader that we have assumed a single baroclinic mode, but $\varepsilon$ does vary among baroclinic modes through $\lambda$. The wave dispersiveness also varies spatially throughout the ocean (figure 1). We calculate $\varepsilon$ for the first four baroclinic modes from observations as described in Appendix A. Except for the high latitudes, the first and second baroclinic modes are almost entirely in the strongly dispersive regime ($\varepsilon \gg 1$). Higher baroclinic modes are to be more weakly dispersive, with $\varepsilon <1$ almost everywhere for mode 4. For a given baroclinic mode, low-latitude regions are more strongly dispersive, while higher latitudes and western boundary currents are more weakly dispersive.

Figure 1. Wave dispersiveness $\varepsilon ^2 = |\,f| \lambda ^2 / \varPsi$ plotted throughout the ocean for the first four baroclinic modes, with the deformation radius $\lambda$ estimated from hydrography and the streamfunction magnitude $\varPsi$ from altimetry. The equatorial band is blocked out because the mean flow amplitude cannot be estimated with confidence there.

Note that (2.1) is a Schrödinger equation. This parallel is made clear if we write (2.1) as

(2.2a,b)\begin{equation} {\rm i}\frac{\partial{\phi}}{\partial{t}} = H\phi, \quad H ={-}\frac{\varepsilon^2}{2}\nabla^2 - {\rm i} \operatorname{\textit{J}}(\psi, {\cdot} ) + \frac{\zeta}{2} . \end{equation}

The operator $H$ is known as the Hamiltonian operator. While the presence of first derivatives in the Hamiltonian stemming from advection may be unfamiliar to some, such terms arise in quantum mechanics when describing a charged particle in a magnetic field. This analogy to quantum mechanics was pointed out by Balmforth, Llewellyn Smith & Young (Reference Balmforth, Llewellyn Smith and Young1998) and we will exploit it here extensively. Rocha et al. (Reference Rocha, Wagner and Young2018) also used this analogy to derive the equivalent of Ehrenfest's theorem for NIWs, while Danioux, Vanneste & Bühler (Reference Danioux, Vanneste and Bühler2015) explained the concentration of NIWs into anticyclones via the analogue of quantum conservation laws. While we are setting $\beta =0$ in this paper, we note that the quantum analogue to $\beta \neq 0$ is known as the ‘Wannier–Stark ladder’, where the potential due to the mesoscale vorticity modulates a linear ramp due to $\beta$ (Balmforth & Young Reference Balmforth and Young1999).

The operator $H$ is Hermitian, i.e.

(2.3)\begin{equation} \int \varphi^* H \phi \, \mathrm{d}^2 \boldsymbol{x} = \int (H \varphi)^* \phi \, \mathrm{d}^2 \boldsymbol{x} \end{equation}

for sufficiently regular functions $\varphi$ and $\phi$ so it has real eigenvalues. We also assume $H$ is compact so that the eigenmodes form a complete orthonormal basis. Let $\boldsymbol {\mu }$ label the eigenmodes $\hat {\phi }_{\boldsymbol {\mu }}(x,y)$ and associated eigenvalues $\omega _{\boldsymbol {\mu }}$ of the operator $H$,

(2.4)\begin{equation} H \hat{\phi}_{\boldsymbol{\mu}} = \omega_{\boldsymbol{\mu}} \hat{\phi}_{\boldsymbol{\mu}}. \end{equation}

The field $\phi$ can then be expanded in the eigenmode basis as

(2.5)\begin{equation} \phi(x,y,t) = \sum_{\boldsymbol{\mu}} a_{\boldsymbol{\mu}}(t)\hat{\phi}_{\boldsymbol{\mu}}(x,y), \end{equation}

where $a_{\boldsymbol {\mu }}(t)$ is the projection of $\phi$ onto the eigenmode $\hat \phi _{\boldsymbol {\mu }}$. The coefficients $a_{\boldsymbol {\mu }}(t)$ then evolve according to

(2.6)\begin{equation} \frac{\mathrm{d}{a_{\boldsymbol{\mu}}}}{\mathrm{d}{t}} ={-}{\rm i}\omega_{\boldsymbol{\mu}}a_{\boldsymbol{\mu}}, \quad \text{so} \quad a_{\boldsymbol{\mu}}(t)=a_{\boldsymbol{\mu}}(0)\exp({-{\rm i}\omega_{\boldsymbol{\mu}}t}). \end{equation}

Therefore, the eigenvalue $\omega _{\boldsymbol {\mu }}$ represents the frequency shift of the mode away from $f$. The total dimensional NIW frequency is hence given by $f(1+\textit {Ro}\,\omega _{\boldsymbol {\mu }})$, where $\textit {Ro}=\varPsi /fL^2$ is the Rossby number. Furthermore, because the eigenvalues are real and the modes are orthogonal, the kinetic energy of the waves is conserved.

We consider this problem on a doubly periodic domain with size $2{\rm \pi} \times 2{\rm \pi}$. This is intended to represent a local view of an ocean that is filled with a random sea of eddies. The solutions we calculate are perfectly periodic and extend across all eddies. In reality, of course, the background field is not perfectly periodic and this causes the solutions to become localised in certain regions. Therefore, the solutions we calculate on the $2{\rm \pi} \times 2{\rm \pi}$ should be thought of similarly.

As a key example in this paper, we consider a $2{\rm \pi} \times 2{\rm \pi}$ domain that contains a dipole vortex given by (figure 2; cf. Asselin et al. Reference Asselin, Thomas, Young and Rainville2020)

(2.7)\begin{equation} \psi = \tfrac{1}{2}\left(\sin{x}-\sin{y}\right). \end{equation}

The analysis below is general, however, and can be applied to more general background flows.

Figure 2. Dipole vorticity with an anticyclone in the upper left corner and a cyclone in the lower right corner. The contours depict the streamfunction with positive values denoted by solid lines and negative values denoted by dashed lines.

2.2. Numerical calculation of eigenvalues and eigenmodes

For most choices of the background flow $\psi$, analytical solutions for the eigenfunctions of $H$ do not exist and numerical solutions are required. Solving the eigenvalue equation numerically requires us to discretise the operator $H$. The discrete eigenfunction is expressed as a vector, and the problem reduces to finding the eigenvalues of a finite matrix. The operator $H$ is Hermitian and so it is desirable for any discrete representation of $H$ to also be Hermitian. A fourth-order central finite difference scheme for the Laplacian term preserves this property. More care is required for the advection operator, for which we use the enstrophy-conserving scheme from Arakawa (Reference Arakawa1966) to preserve the Hermitian nature of the operator and guarantee that the eigenvalues of the matrix are real. Having real eigenvalues ensures that the conservation of NIW kinetic energy is respected in the discrete system. The exact method of numerically solving the eigenvalue problem is detailed in Appendix B.

3.1. The leading-order mode

Near-inertial waves are forced by atmospheric storms, which have a much larger horizontal scale than mesoscale eddies and can be idealised as a uniform forcing. We assume that the result of this forcing is the excitation of a constant non-zero $\phi$. The projection of this initial condition onto a given mode can thus be found by integrating that mode across the domain. For plane waves, a domain integral will vanish unless $\boldsymbol {\mu } = 0$, such that a uniform forcing will only project onto the $\boldsymbol {\mu } = 0$ mode in the unperturbed case. We begin by focusing on that case to obtain expressions for the perturbations to its spatial structure as well as its frequency shift. A small part of the forcing, however, projects onto modes with $\boldsymbol {\mu } \neq 0$ and we will return to these higher modes below.

The leading-order solution for $\boldsymbol {\mu } = 0$ is $\hat {\phi }_0^{(0)}=1$ and $\omega ^{(0)}_0 = 0$, and there is no modulation of the waves by the mesoscale eddy field. To obtain this modulation, we must go to higher order. At $O(\varepsilon ^0)$, the eigenvalue problem is

(3.5)\begin{equation} H^{(0)}\hat{\phi}^{(1)}_0 + H^{(1)}\hat{\phi}^{(0)}_0 = \omega^{(0)}_0\hat{\phi}^{(1)}_0 + \omega^{(1)}_0\hat{\phi}^{(0)}_0.\end{equation}

With $\omega _0^{(0)}=0$ and the advection term in $H^{(1)}$ vanishing when acting on $\hat {\phi }_0^{(0)} = 1$, this reduces to

(3.6)\begin{equation} {-}\frac{1}{2}\nabla^2\hat{\phi}_{0}^{(1)} + \frac{\zeta}{2} = \omega_0^{(1)}. \end{equation}

The two terms on the left vanish when integrated over the doubly periodic domain, so we conclude that $\omega _0^{(1)} = 0$.

There is, however, a correction to the eigenfunction at this order, determined by

(3.7)\begin{equation} \nabla^2\hat{\phi}_0^{(1)} = \nabla^2\psi. \end{equation}

With periodic boundary conditions, the solution to this is

(3.8)\begin{equation} \hat{\phi}_0^{(1)} = \psi,\end{equation}

where we have assumed that $\psi$ is defined such that it has zero domain average. This recovers the expression for $\hat {\phi }$ from YBJ. The structure of the mesoscale eddy field is imprinted onto the waves by the $\varepsilon ^{-2}\hat {\phi }^{(1)}_0$ term. Because the modulation is by the real streamfunction $\psi$, only the NIW amplitude is modulated by mesoscale eddies. The NIW field remains in phase across the domain.

We now also seek the leading non-zero correction to the eigenvalue, for which we go up another order. The eigenvalue equation at $O(\varepsilon ^{-2})$ is

(3.9)\begin{equation} H^{(0)}\hat{\phi}_0^{(2)} + H^{(1)}\hat{\phi}_0^{(1)} = \omega_0^{(0)}\hat{\phi}_0^{(2)} + \omega_0^{(1)}\hat{\phi}_0^{(1)} + \omega_0^{(2)}\hat{\phi}_0^{(0)}. \end{equation}

With $\omega _0^{(0)}=\omega _0^{(1)}=0$ and $\operatorname {\textit {J}}(\psi,\psi )=0$, this simplifies to

(3.10)\begin{equation} {-}\tfrac{1}{2}\nabla^2\hat{\phi}_0^{(2)} + \tfrac{1}{2} \psi \nabla^2 \psi = \omega_0^{(2)}. \end{equation}

The first term on the left vanishes under domain integration. Integrating the second term on the left by parts yields

(3.11)\begin{equation} \omega_0^{(2)} ={-}\frac{1}{2}\frac{\displaystyle\int |\boldsymbol{\nabla}\psi|^2 \, \mathrm{d}^2 \boldsymbol{x}}{\displaystyle\int \mathrm{d}^2 \boldsymbol{x}}. \end{equation}

The leading-order frequency shift is $\varepsilon ^{-2} \omega _0^{(2)}$. Given that $\varepsilon ^{-2} \ll 1$, the frequency shift away from $f$ is suppressed substantially, even compared with the small frequency shift assumed from the outset. Re-dimensionalising the expression results in

(3.12)\begin{equation} \omega_0^{(2)} ={-}\frac{1}{2f_0\lambda^2}\frac{\displaystyle\int |\boldsymbol{\nabla}\psi|^2 \, \mathrm{d}^2 \boldsymbol{x}}{\displaystyle\int \mathrm{d}^2 \boldsymbol{x}}. \end{equation}

This agrees with the YBJ result for the dispersion relation in the strong-dispersion regime, indicating that the frequency shift is proportional to the average kinetic energy of the eddy field.

3.2. Higher-order modes

We now return to the higher modes with $\boldsymbol {\mu }\neq 0$. These modes are degenerate to leading order. For example, the modes $(1,0),(-1,0),(0,1)$ and $(0,-1)$ all have $\omega ^{(0)}_{\boldsymbol {\mu }}=\frac {1}{2}$. Degenerate perturbation theory is necessary to calculate the first-order corrections to the eigenvalues and eigenfunctions (e.g. Sakurai & Napolitano Reference Sakurai and Napolitano2020). To obtain these corrections, we proceed naively with the calculation. We will run into a contradiction that motivates us to choose a different basis set than was chosen in (3.4a,b). To those familiar with degenerate perturbation theory, this may seem unnecessary, but we believe it to be more pedagogical.

We again start from the $O(\varepsilon ^0)$ equation, which now reads

(3.13)\begin{equation} H_0\hat{\phi}_{\boldsymbol{\mu}}^{(1)}+H_1\hat{\phi}_{\boldsymbol{\mu}}^{(0)}=\omega_{\boldsymbol{\mu}}^{(0)}\hat{\phi}_{\boldsymbol{\mu}}^{(1)} + \omega_{\boldsymbol{\mu}}^{(1)}\hat{\phi}_{\boldsymbol{\mu}}^{(0)}. \end{equation}

Multiplying this equation by $\hat {\phi }_{\boldsymbol {\nu }}^{(0)*}$, with both $\boldsymbol {\nu }$ and $\boldsymbol {\mu }$ labelling one of the modes in the degenerate group, and integrating over the domain results in

(3.14)\begin{equation} \int\hat{\phi}_{\boldsymbol{\nu}}^{(0)*}\left(H_0-\omega_{\boldsymbol{\mu}}^{(0)}\right)\hat{\phi}_{\boldsymbol{\mu}}^{(1)} \, \mathrm{d}^2 \boldsymbol{x} = \omega_{\boldsymbol{\mu}}^{(1)}\int \hat{\phi}_{\boldsymbol{\nu}}^{(0)*}\hat{\phi}_{\boldsymbol{\mu}}^{(0)} \, \mathrm{d}^2 \boldsymbol{x} - \int \hat{\phi}_{\boldsymbol{\nu}}^{(0)*} H_1 \hat{\phi}_{\boldsymbol{\mu}}^{(0)} \, \mathrm{d}^2 \boldsymbol{x}. \end{equation}

Using integration by parts, the $H_0$ on the left can be swapped for $\omega _{\boldsymbol {\nu }}^{(0)}$. Because the modes are degenerate to this order, the left-hand side vanishes. Furthermore, using the orthonormality of the eigenfunctions, the corrections to the eigenvalues are determined by

(3.15)\begin{equation} \omega_{\boldsymbol{\mu}}^{(1)} \delta_{\boldsymbol{\nu}\boldsymbol{\mu}} = \frac{1}{4{\rm \pi}^2} \int \hat{\phi}_{\boldsymbol{\nu}}^{(0)*} H_1 \hat{\phi}_{\boldsymbol{\mu}}^{(0)} \, \mathrm{d}^2 \boldsymbol{x}. \end{equation}

We have now arrived at our contradiction. The left-hand side of this equation is diagonal, whereas the right-hand side is not necessarily so. In (3.4a,b) we chose a basis for the unperturbed eigenfunctions: $\{{\rm e}^{{\rm i}x},{\rm e}^{-{\rm i}x},{\rm e}^{{\rm i}y},{\rm e}^{-{\rm i}y}\}$ for $|\boldsymbol {\mu }| = 1$. The key to degenerate perturbation theory is to choose a basis of the degenerate space to avoid this contradiction. It is clear that the correct basis must diagonalise $H_1$, which our original choice does not.

To proceed, we calculate the right-hand side of (3.15) in the original basis. This results in a $4 \times 4$ matrix. We diagonalise this matrix and find the corresponding linear combination of the original basis functions that diagonalises $H_1$. The corresponding eigenfunction corrections can be found by solving the screened Poisson equation obtained from the first-order equation (3.13). Rearranging, we have

(3.16)\begin{equation} \left( H_0-\omega_{\boldsymbol{\mu}}^{(0)} \right) \hat{\phi}_{\boldsymbol{\mu}}^{(1)} = \left(\omega_{\boldsymbol{\mu}}^{(1)}-H_1\right)\hat{\phi}_{\boldsymbol{\mu}}^{(0)}, \end{equation}

where the $\hat {\phi }_{\boldsymbol {\mu }}^{(0)}$ should be in the basis diagonalising $H_1$. If $H_1$ identically vanishes in this subspace, the degeneracy must be lifted at the next order, as in the example below.

3.3. Dipole flow solutions

We now consider the specific example of the dipole flow (2.7). Numerical solutions for $\varepsilon =2$ show that a uniform initial condition projects strongly (98.5 % of the energy) onto the $\hat {\phi }_0$ mode (figure 3). There is a small but negative frequency shift of $\omega _0 = -0.03104$. This agrees excellently with the predicted frequency shift from (3.11) of $\varepsilon ^{-2} \omega _0^{(2)} = \frac {1}{32} = -0.03125$. Additionally, there is weak horizontal structure that aligns with the streamfunction as expected. The root-mean-squared error between the numerical eigenmode and the analytical eigenmode $\hat {\phi }_0^{(0)} + \varepsilon ^{-2} \hat {\phi }_0^{(1)} = 1 + \varepsilon ^{-2} \psi$ is 1 %. The agreement is excellent despite $\varepsilon$ not being particularly large.

Figure 3. Numerical solution to the eigenvalue problem (2.4) with $\varepsilon = 2$ for the dipole flow. A uniform forcing primarily projects onto the mode shown on the left, with the mode with the second-highest projection shown on the right only making up less than 2 % of the energy. The eigenvalues $\omega$ and projection fractions (of energy) are shown in the panel titles. Vectors show the corresponding NIW velocities.

For the dipole flow, the right-hand side of (3.15) is zero for all combinations of basis functions of the $\omega _{\boldsymbol {\mu }}^{(0)} = \frac {1}{2}$ subspace. Therefore, there are no first-order frequency shifts, $\omega _{\boldsymbol {\mu }}^{(1)} = 0$, and the degeneracy is not lifted at this order. Performing the same procedure that led to (3.15) on the second-order equation yields

(3.17)\begin{equation} \omega_{\boldsymbol{\mu}}^{(2)} \delta_{\boldsymbol{\nu}\boldsymbol{\mu}} = \frac{1}{4{\rm \pi}^2} \int \hat{\phi}_{\boldsymbol{\nu}}^{(0)*} H_1 \hat{\phi}_{\boldsymbol{\mu}}^{(1)} \, \mathrm{d}^2 \boldsymbol{x}. \end{equation}

For our trial basis consisting of the four plane waves, we solve the screened Poisson equation (3.16) for the corresponding $\hat {\phi }_{\boldsymbol {\mu }}^{(1)}$. This is tedious but doable because the right-hand side is just a sum of sines and cosines. The equation for the second-order frequency shift can be diagonalised, and this time the eigenvalues are not zero and the degeneracy is lifted. We find for $\omega _{\boldsymbol {\mu }}^{(2)}$ the values $-\frac {1}{96}$, $-\frac {7}{96}$, $-\frac {49}{96}$ and $-\frac {55}{96}$, only the first of which corresponds to an eigenfunction that the forcing projects onto at this order. The leading-order eigenfunction of that mode is $\hat {\phi }_{\boldsymbol {\mu }}^{(0)} = -\psi$ (figure 3). The eigenvalue $\varepsilon ^2 \omega _{\mu }^{(0)} + \varepsilon ^{-2} \omega _{\boldsymbol {\mu }}^{(2)} = 1.99739$ is again in excellent agreement with the numerical eigenvalue of $1.99729$.

In this regime, the horizontal structure in the waves primarily arises due to $\hat {\phi }_0^{(1)}$, which is suppressed by $O(\varepsilon ^{-2})$. There is also horizontal structure due to modes with $\boldsymbol {\mu }\neq 0$, but these are projected onto weakly; the fraction of the variance accounted for by such a mode is $O(\varepsilon ^{-4})$ (Sakurai & Napolitano Reference Sakurai and Napolitano2020). As such, the wave potential energy, which depends on horizontal gradients in the wave field, is also suppressed. Xie & Vanneste (Reference Xie and Vanneste2015) associated the generation of wave potential energy with a sink of the background eddy kinetic energy in a process known as stimulated generation. Given the weak generation of horizontal structure, stimulated generation is weak in the strong-dispersion regime.

4.1. Parallel shear flow

We begin with an example of a parallel shear flow in which the streamfunction $\psi$ is a function of $x$ only. The symmetry in $y$ means the problem reduces to a one-dimensional eigenvalue problem. Balmforth et al. (Reference Balmforth, Llewellyn Smith and Young1998) considered this problem for a specific example of a shear flow that can be solved in closed form. Zhang & Xie (Reference Zhang and Xie2023) considered the limits of strong and weak dispersion for the same mean flow. Here, we address how the weak-dispersion limit can be analysed for a general parallel shear flow and apply the procedure to the example flow from Balmforth et al. (Reference Balmforth, Llewellyn Smith and Young1998). Our goal is to calculate the structure of the eigenmodes and their corresponding eigenvalues. We begin by introducing the WKB method from which the two scalings arise. For the anisotropic scaling, the key result is (4.8); for the isotropic scaling, the equivalent result is (4.12).

We assume that the streamfunction $\psi (x)$ is periodic on the domain $[-{\rm \pi},{\rm \pi} ]$. The eigenvalue problem (2.4) reduces to

(4.1)\begin{equation} {-}\frac{\varepsilon^2}{2}\nabla^2\hat{\phi} - {\rm i}v\frac{\partial{\hat{\phi}}}{\partial{y}} + \frac{\zeta}{2}\hat{\phi}= \omega\hat{\phi}, \end{equation}

where $\zeta = \nabla ^2 \psi$ and $v = \partial _x \psi$ are both functions of $x$ only, and we have suppressed the label on the eigenmode. The coefficients are independent of $y$, which motivates the ansatz $\hat {\phi } = \varPhi (x){\rm e}^{{\rm i}my}$. Given that the domain has width $2{\rm \pi}$ in $y$, the wavenumber $m$ must be an integer. With this ansatz, we are left with the one-dimensional eigenvalue problem

(4.2)\begin{equation} {-}\frac{\varepsilon^2}{2}\frac{\mathrm{d}^2 \varPhi}{\mathrm{d}\kern0.7pt x^2} + \left( \frac{\varepsilon^2 m^2}{2} + m v + \frac{\zeta}{2}\right) \varPhi = \omega\varPhi . \end{equation}

This is the Schrödinger equation of a particle in one-dimensional potential, with the bracketed term playing the role of the potential (Balmforth et al. Reference Balmforth, Llewellyn Smith and Young1998).

As $\varepsilon$ is small, WKB analysis can be used to find approximations to the eigenvalues and eigenfunctions (e.g. Bender & Orszag Reference Bender and Orszag1999). In WKB theory the field $\varPhi$ is expanded as

(4.3)\begin{equation} \varPhi(x) = \exp \frac{1}{\delta} \sum_{j = 0}^\infty \delta^j S_j(x), \end{equation}

where $\delta \ll 1$ is a scaling parameter that we are yet to determine. Substituting this into (4.2) yields

(4.4)\begin{equation} {-}\frac{\varepsilon^2}{2}\left[\frac{1}{\delta^2}\left(\sum_{j = 0}^\infty \delta^j\frac{\mathrm{d}{S_j}}{\mathrm{d}\kern0.7pt{x}}\right)^2+\frac{1}{\delta} \sum_{j = 0}^\infty \delta^j\frac{\mathrm{d}{^2 S_j}}{\mathrm{d}\kern0.7pt{x^2}}\right]+\frac{\varepsilon^2m^2}{2}+mv+\frac{\zeta}{2}=\omega. \end{equation}

If we assume that $m \sim O(1)$, both the refraction term and the advection terms are $O(1)$, and they must be balanced by a dispersion term of the same order. Requiring the lowest-order dispersion term to be $O(1)$ implies that $\delta = \varepsilon$, and the $O(1)$ equation becomes

(4.5)\begin{equation} {-}\frac{1}{2}\left(\frac{\mathrm{d}{S_0}}{\mathrm{d}\kern0.7pt{x}}\right)^2 + mv + \frac{\zeta}{2} = \omega.\end{equation}

By writing $\varepsilon ^{-1} \,\mathrm {d} S_0/\mathrm {d}\kern0.7pt x = {\rm i} k$, this equation is analogous to the dispersion relation (1.1a{–}c) specialised to this parallel shear flow. The function $S_0$ is found to be

(4.6)\begin{equation} S_0(x) ={\pm} \sqrt{2}{\rm i} \int^x \sqrt{\omega- m v(x') - \frac{\zeta(x')}{2}} \, \mathrm{d}\kern0.7pt x' \end{equation}

and determines the leading-order phase variations of the solution. One can additionally show (Bender & Orszag Reference Bender and Orszag1999, (10.1.12)) that the next-order solution is

(4.7)\begin{equation} S_1(x)={-}\frac{1}{4}\ln{\left(\omega-m v-\frac{\zeta}{2}\right)}, \end{equation}

which determines the leading-order amplitude modulation of the solution.

This asymptotic expansion is valid away from regions where the integrand above is zero. These are known as turning points of the problem and exist if $\omega < \max ( mv + \zeta /2)$. The associated eigenfunctions are referred to as bound states. Near turning points, $\omega -mv-\zeta /2$ can be approximated by a linear function of $x$, and solutions to (4.5) are given by Airy functions. The Airy function solutions must be asymptotically matched to the solutions away from the turning points. This yields an integral constraint from which the eigenvalues $\omega$ can be determined. The problem as formulated above is the classic two-turning point problem, and the asymptotic matching procedure is well documented (e.g. Bender & Orszag Reference Bender and Orszag1999, (10.5.6)). The resulting condition for $\omega$, often referred to as a quantisation condition, is

(4.8)\begin{equation} \frac{\sqrt{2}}{\varepsilon} \int_{x_0}^{x_1} \sqrt{\omega - mv(x) - \frac{\zeta(x)}{2}} \, \mathrm{d}\kern0.7pt x = \left(n+\frac{1}{2}\right){\rm \pi}, \quad \text{with} \ n = 0, 1, \dots,\end{equation}

where $x_0$ and $x_1$ are the turning points. The projection of a uniform forcing onto these modes can also be calculated asymptotically. The domain integral of a mode is dominated by contributions from the turning points (e.g. Bender & Orszag Reference Bender and Orszag1999, (10.4.24)).

If $\omega > \max (\zeta /2 + mv)$ then there are no turning points. The corresponding eigenmodes are referred to as free states, and the quantisation condition is replaced by

(4.9)\begin{equation} \frac{\sqrt{2}}{\varepsilon} \int_{-{\rm \pi}}^{\rm \pi} \sqrt{\omega - mv(x) - \frac{\zeta(x)}{2}} \, \mathrm{d}\kern0.7pt x = 2n{\rm \pi}, \quad \text{with} \ n = 0, 1\ldots\,. \end{equation}

Note the lack of a half-integer shift that for bound states arises from the Airy behaviour near turning points. The lack of turning points in the free states also means (4.6) is valid across the entire domain. Because the eigenfunctions of these free states are oscillatory in the entire domain, a uniform forcing projects only weakly onto them, and we do not discuss them any further. We also note that the discretisation of the free states is due to the periodic domain; they would be replaced by a continuum of modes in an infinite domain, while the bound states would remain discrete.

Under this scaling, the WKB modes are anisotropic. We assumed $m \sim {O}(1)$, which means that the modes’ phase varies in $y$ on a length scale $O(1)$. In contrast, the leading-order phase variations in $x$ come from $\varepsilon ^{-1} S_0$ and, therefore, occur on a scale $O(\varepsilon )$. The phase varies slowly along streamlines and rapidly across streamlines. This makes refraction and advection come in at the same order as cross-streamline dispersion.

An alternative would be to choose the scaling $m \sim O(\varepsilon ^{-2}).$ Repeating the WKB ansatz requires a choice of $\delta =\varepsilon ^2$ and $\omega \sim O(\varepsilon ^{-2})$ in order to end up with an equation of a similar form to (4.5):

(4.10)\begin{equation} {-}\frac{1}{2}\left(\frac{\mathrm{d}{S_0}}{\mathrm{d}\kern0.7pt{x}}\right)^2 + \frac{\varepsilon^4 m^2}{2} + \varepsilon^2 m v = \varepsilon^2 \omega. \end{equation}

With the scaling given above, each term is $O(1)$. We can solve for $S_0$ and the corresponding quantisation condition for bound modes:

(4.11)$$\begin{gather} S_0(x) ={\pm} \sqrt{2}{\rm i} \varepsilon \int^x\sqrt{\omega-\frac{\varepsilon^2m^2}{2}-mv(x')} \, \mathrm{d}\kern0.7pt x', \end{gather}$$

(4.12)$$\begin{gather}\frac{\sqrt{2}}{\varepsilon}\int_{x_0}^{x_1}\sqrt{\omega-\frac{\varepsilon^2m^2}{2}-mv(x)} \, \mathrm{d}\kern0.7pt x = \left(n+\frac{1}{2}\right){\rm \pi}. \end{gather}$$

These modes are isotropic. The phase variations in $y$ occur on a scale $O(\varepsilon ^2)$, which is the same as in $x$ because phase variations in $x$ now come from $\varepsilon ^{-2} S_0$. This makes advection and along-streamline dispersion come in at the same order, and it makes refraction negligible.

Despite the different characteristics of the two scalings, they lead to similar quantisation conditions that differ only by what terms are included. We can combine them into a uniformly valid quantisation condition:

(4.13)\begin{equation} \frac{\sqrt{2}}{\varepsilon}\int_{x_0}^{x_1}\sqrt{\omega- \frac{\varepsilon^2m^2}{2} - mv(x) - \frac{\zeta(x)}{2}} \, \mathrm{d}\kern0.7pt x = \left(n+\frac{1}{2}\right){\rm \pi}. \end{equation}

The ‘potential’ governing the wave evolution is therefore

(4.14)\begin{equation} V(x) = \frac{\varepsilon^2m^2}{2} + mv(x) + \frac{\zeta(x)}{2}. \end{equation}

Under the anisotropic scaling $m \sim O(1)$, the along-streamline dispersion term is suppressed by a factor $\varepsilon ^{2}$, leaving the $O(1)$ refraction and advection terms to dominate. Under the isotropic scaling $m \sim O(\varepsilon ^{-2})$, the advection and along-streamline dispersion terms are enhanced by a factor $\varepsilon ^{-2}$ and dominate over a now negligible refraction term. In both cases, the general equation is obtained by retaining a term that is of higher order, which is allowed in an asymptotic theory. A uniform forcing only projects onto modes with $m = 0$, so all the modes projected onto are of the anisotropic variety.

We now consider a specific example of a parallel shear flow that varies sinusoidally in $x$:

(4.15)\begin{equation} \psi = \cos{x}. \end{equation}

This shear flow has a region of anticyclonic vorticity at the centre of the domain and cyclonic vorticity centred on $x = \pm {\rm \pi}$ (figure 4a,b). This is a rare example in which the eigenvalue problem (4.2) can be solved exactly using Mathieu functions (Balmforth et al. Reference Balmforth, Llewellyn Smith and Young1998). The generally applicable WKB theory described above accurately predicts the eigenvalues, even for a modestly small $\varepsilon = \frac {1}{4}$ (figure 4c). We provide the analytical solutions to the WKB integrals in Appendix C. We also note that the symmetry of the problem means that a uniform wind forcing only projects onto modes with even $n$.

Figure 4. (a) Streamfunction and flow vectors for the shear flow example. (b) Vorticity, showing the anticyclonic vorticity concentrated in the centre of the domain and cyclonic vorticity on the outside. (c) Eigenvalues $\omega$ as a function of the integer wavenumber $m$ for $\varepsilon = \frac {1}{4}$. The results from WKB theory (orange crosses) are shown along with the exact eigenvalues found from numerical solutions (black circles). The WKB results agree remarkably well with the numerical results, although there are some spurious eigenvalues near the boundary between free and bound modes. The purple shading indicates the region where free modes exist, which are not shown here.

For $m = 0$, the eigenmodes are shaped by the potential $V = {\zeta }/{2}$ (figure 5a). Where $\omega > V$, $S_0$ is imaginary and the solutions are oscillatory; where $\omega < V$, $S_0$ is real and the solutions are decaying (figure 5a). Near the anticyclonic centre of the flow, the potential is at its lowest and all the modes are oscillatory. Moving further out into the cyclonic region, more and more of the modes become evanescent.

Figure 5. (a) The potential $V = {\zeta }/{2}$ (black line) of the parallel shear flow for the $m=0$ mode. The dashed lines show the level of each eigenvalue $\omega$ for $\varepsilon = \frac {1}{4}$. The solid coloured lines represent a scaled representation of each eigenfunction corresponding to a given eigenvalue, as identified by the colours. (b) Time evolution of the NIW amplitude $|\phi |$ for the parallel shear flow example with $\varepsilon =\frac {1}{4}$, starting by a uniform field.

The dependence of the potential on the vorticity $\zeta$ leads to trapping of NIWs in anticyclones. The trapping arises from the dephasing of the modes that make up the initial condition. This is analogous to the argument in Gill (Reference Gill1984) regarding the vertical propagation of NIWs due to the $\beta$ effect. The evolution occurs in three phases (figure 5b). First, refraction imprints the mesoscale vorticity onto the initially uniform wave phase, leaving the amplitude unchanged (cf. Asselin et al. Reference Asselin, Thomas, Young and Rainville2020). Second, once these phase gradients are sufficiently pronounced, cross-streamline dispersion becomes important and concentrates the wave energy into the centre of the anticyclone, and the wave field assumes a spatial scale $O(\varepsilon )$. Third, the amplitude remains elevated on average within the anticyclonic region but is more spread out than during the initial concentration. It is this long-time behaviour that corresponds to the fully dephased eigenmodes. The time it takes for this dephasing to occur depends inversely on the spacing of the eigenvalues $\omega$. As $\varepsilon$ decreases, the eigenvalues become more finely spaced. The $m = 0$ modes have a spacing $O(\varepsilon )$, so it takes $t \sim O(\varepsilon ^{-1})$ for them to dephase. Another way to think of this is that as $\varepsilon$ decreases, dispersion becomes weaker and it takes longer for phase gradients to build up to a level where dispersion is important.

Finally, we note that Asselin et al. (Reference Asselin, Thomas, Young and Rainville2020) discussed solutions to the YBJ equation for which phase lines are aligned with streamlines and straining is ineffective in driving a decrease in the spatial scale of the wave. That analysis sets the dispersion term to zero, however, and so only captures the initial phase in which refraction dominates. The WKB theory presented above shows that cross-streamline dispersion is of leading order and should not be dropped if the long-term evolution is of interest (cf. figure 5b). Our $m=0$ anisotropic modes can thus be understood as a generalisation of Asselin et al.'s (Reference Asselin, Thomas, Young and Rainville2020) solution. The ineffectiveness of straining due to the alignment of the wave phase with streamlines remains apparent, but cross-streamline dispersion is now taken into account such that the solution remains valid at late times. We further note that our anisotropic modes also allow for slow variations of the wave field along streamlines, such that advection assumes the same importance as refraction and cross-streamline dispersion. These anisotropic modes with $m>0$ may be excited by a non-uniform forcing, such as a passing atmospheric front (cf. Thomas, Smith & Bühler Reference Thomas, Smith and Bühler2017).

4.2. Axisymmetric flow

We now consider a streamfunction with axial symmetry, such that $\psi = \psi (r)$, where $r$ is the radial distance from the origin. Llewellyn Smith (Reference Llewellyn Smith1999) studied NIWs with azimuthal wavenumber zero in an axisymmetric vortex and provided asymptotic expressions for the frequency of the lowest radial mode. Kafiabad, Vanneste & Young (Reference Kafiabad, Vanneste and Young2021) studied a similar case but also considered the impact of NIWs back on the vortex. Using WKB theory, we consider NIWs with an arbitrary azimuthal wavenumber and provide a transcendental equation that can be solved for their frequency as for the parallel shear flows above.

We make the ansatz $\hat {\phi }=A(r){\rm e}^{{\rm i}m\theta }$, where $\theta$ is the azimuthal angle, and again we drop the mode label. In polar coordinates, (2.4) then reduces to

(4.16)\begin{equation} {-}\frac{\varepsilon^2}{2} \left( \frac{\mathrm{d}{^2A}}{\mathrm{d}{r^2}} + \frac{1}{r}\frac{\mathrm{d}{A}}{\mathrm{d}{r}} \right) + \left(\frac{\varepsilon^2 m^2}{2r^2} + \frac{mv}{r} + \frac{\zeta}{2} \right) A = \omega A, \end{equation}

where $v = \partial _r \psi$ denotes the azimuthal velocity. There are some subtleties involved in applying WKB theory to this equation. For modes with $m>0$, the potential diverges at the origin. This issue has long been noted in the quantum mechanics literature and can be addressed by performing a so-called Langer transform on the equation. For $m=0$, there is no divergence of the potential, but there is a phase shift at the origin. As pointed out by Berry & Ozorio de Almeida (Reference Berry and Ozorio de Almeida1973), both cases turn out to give the same quantisation condition:

(4.17)\begin{equation} \frac{\sqrt{2}}{\varepsilon}\int_{r_0}^{r_1}\sqrt{\omega-V(r)} \, \mathrm{d} r = \left(n+\frac{1}{2}\right){\rm \pi}, \quad \text{with} \ n = 0, 1, 2, \dots. \end{equation}

Here the potential is

(4.18)\begin{equation} V(r) = \frac{\varepsilon^2m^2}{2r^2}+ \frac{mv}{r} + \frac{\zeta}{2}. \end{equation}

If $m > 0$, the integration bounds $r_0$ and $r_1$ are the two zeros of the integrand; if $m = 0$, $r_0 = 0$ and $r_1$ is the one zero of the integrand. As in the case of a parallel shear flow, this expression is uniformly valid in the sense that it works for both $m \sim O(1)$ and $m \sim O(\varepsilon ^{-2})$. These again correspond to anisotropic and isotropic modes, respectively, with refraction, advection and dispersion along and across streamlines playing the same roles as before. The only difference is that the streamlines are now circular.

We consider the concrete example of an isolated Gaussian vortex on an infinite domain:

(4.19)\begin{equation} \psi(r)={\rm e}^{-{r^2}/{4}}. \end{equation}

This corresponds to an anticyclone in the centre of the domain that is surrounded by a halo of cyclonic vorticity (figure 6a,b). Again, the WKB calculation for $\varepsilon = \frac {1}{4}$ yields eigenvalues that agree extremely well with the exact eigenvalues (figure 6c). The structure of the first few modes is shown in figure 7. For $m=0$, the modes are concentrated in the anticyclone. For $m>0$, there is a repulsion from the very centre of the anticyclone due to the advection and dispersion terms in $V(r)$. This repulsion increases with $m$, but the modes remain primarily concentrated in the region of anticyclonic vorticity. The modes become more isotropic as $m$ is increased. Note that again, a uniform forcing only projects onto the $m=0$ mode due to the symmetry of the vorticity field. We also note that only the bound states form a discrete spectrum in an infinite domain, there will also be a continuum of free modes.

Figure 6. (a) Streamfunction and flow vectors for the axisymmetric flow example. (b) Vorticity field showing the anticyclonic vorticity concentrated in the centre of the domain, which is flanked by a halo of cyclonic vorticity. (c) Eigenvalues $\omega$ as a function of azimuthal wavenumber $m$ for $\varepsilon = \frac {1}{4}$. The results from WKB theory (orange crosses) are shown along with the exact eigenvalues found from numerical solutions (black circles). The WKB approximation agrees remarkably well with the numerical results.

Figure 7. Real part of the eigenfunctions for the axisymmetric flow example with $\varepsilon = \frac {1}{4}$. The radial wavenumber $n$ increases from left to right and corresponds to an increasing number of nodes in the radial direction. The azimuthal wavenumber increases from top to bottom and corresponds to an increasing number of nodes in the azimuthal direction.

4.3. General case

Based on the intuition gained above, we wish to construct a uniformly valid asymptotic expansion for a general two-dimensional background flow. We again make a WKB ansatz that leads to (4.23). This equation can be solved by the method of characteristics and recovers Kunze's ray tracing. The quantisation condition for the general case is (4.27a,b).

In analogy with the isotropic scaling, we begin by assuming a solution of the form

(4.20)\begin{equation} \hat{\phi}(x,y) = \exp\left[{\frac{1}{\varepsilon^2}\sum_{j=0}^\infty \varepsilon^{2j}S_{j}(x,y)}\right], \end{equation}

again dropping the mode label. Substituting this into (2.4) yields

(4.21)\begin{equation} {-}\frac{1}{2\varepsilon^2}\left|\sum_{j=0}^\infty \varepsilon^{2j}\boldsymbol{\nabla} S_j\right|^2 - \frac{1}{2}\sum_{j=0}^\infty\varepsilon^{2j}\nabla^2S_j - \frac{{\rm i}}{\varepsilon^2}\sum_{j=0}^\infty\varepsilon^{2j}\operatorname{\textit{J}}(\psi,S_j) + \frac{\zeta}{2} = \omega. \end{equation}

Assuming $\omega \sim O(\varepsilon ^{-2})$ and collecting leading-order terms, we obtain

(4.22)\begin{equation} {-}\tfrac{1}{2}|\boldsymbol{\nabla} S_0|^2 - {\rm i}\operatorname{\textit{J}}(\psi,S_0) = \varepsilon^2 \omega. \end{equation}

In the simple examples discussed above, we obtained a uniformly valid approximation by retaining the higher-order refraction term in the leading-order equation arising from an isotropic scaling. We do so again here:

(4.23)\begin{equation} {-}\frac{1}{2}|\boldsymbol{\nabla} S_0|^2 - {\rm i}\operatorname{\textit{J}}(\psi,S_0) +\frac{\varepsilon^2\zeta}{2} = \varepsilon^2 \omega. \end{equation}

We anticipate that the order of these terms again changes for anisotropic modes. If the phase varies slowly along streamlines, the advection term is reduced by a factor $O(\varepsilon ^2)$ and cross-streamline dispersion, acting on spatial variations on a scale of $O(\varepsilon )$ rather than $O(\varepsilon ^2)$, will attain the same order, whereas along-streamline dispersion becomes negligible. Equation (4.23) can therefore capture both isotropic and anisotropic modes.

We now introduce the wavenumber vector $\boldsymbol {k}$ by writing $\varepsilon ^{-2} \partial S_0 / \partial \boldsymbol {x} = {\rm i} \boldsymbol {k}$. Equation (4.23) can be solved using the method of characteristics:

(4.24a{–}c)\begin{equation} \frac{\mathrm{d}\kern0.7pt{\boldsymbol{x}}}{\mathrm{d}{\tau}} = \varepsilon^2 \boldsymbol{k} + \boldsymbol{u}, \quad \frac{\mathrm{d}{\boldsymbol{k}}}{\mathrm{d}{\tau}} ={-}\frac{\partial{}}{\partial{\boldsymbol{x}}} \left( \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{k} + \frac{\zeta}{2} \right), \quad \omega = \frac{\varepsilon^2 |\boldsymbol{k}|^2}{2} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{k} + \frac{\zeta}{2}. \end{equation}

These are the non-dimensionalised ray tracing equations of Kunze (Reference Kunze1985). We further elaborate on this connection between YBJ and Kunze's ray tracing below.

Numerical solutions for the dipole flow show that the majority of a uniform forcing projects onto anisotropic modes that show little structure along streamlines and vary more rapidly across streamlines (figure 8). With $\varepsilon = \frac {1}{4}$ there is also some projection onto modes that show more characteristics of isotropic phase variations. The variations are more rapid, as emerges from the isotropic scaling discussed above.

Figure 8. Real part of the eigenfunctions of the dipole flow with $\varepsilon = \frac {1}{4}$. Together, these eight eigenfunctions represent over 97 % of the energy excited by a uniform forcing. They are the eight modes with the strongest projection and are then ordered by eigenvalue $\omega$. The eigenvalues are shown in the top left corner and the projections of a uniform forcing onto the eigenfunction (energy fraction) are shown in the top right corner.

Finally, we show how approximations to the eigenvalues can be obtained in the weak-dispersion limit when the flow problem is not separable, as it was in the cases of a parallel shear flow or axisymmetric flow. To this end, we utilise results from the quantum mechanics literature. Recall that the YBJ equation is equivalent to the Schrödinger equation, with the YBJ operator

(4.25)\begin{equation} H ={-}\frac{\varepsilon^2}{2} \nabla^2 - {\rm i} \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} + \frac{\zeta}{2} \end{equation}

playing the role of the Hamiltonian. The weak-dispersion limit corresponds to the classical limit of the equivalent quantum system, and the ray tracing equations are the analogue of the classical Hamiltonian dynamics. The classical Hamiltonian is obtained from $H$ by making the substitution $\boldsymbol {\nabla } \mapsto {\rm i}\boldsymbol {k}$, yielding the dispersion relation in (4.24a{–}c). The Hamiltonian dynamics is then

(4.26a,b)\begin{equation} \frac{\mathrm{d}\kern0.7pt{\boldsymbol{x}}}{\mathrm{d}{\tau}} = \frac{\partial{\omega}}{\partial{\boldsymbol{k}}} \quad \text{and} \quad \frac{\mathrm{d}{\boldsymbol{k}}}{\mathrm{d}{\tau}} ={-}\frac{\partial{\omega}}{\partial{\boldsymbol{x}}}, \end{equation}

the ray tracing equations stated in (4.24a{–}c). The connection with the Schrödinger equation is most easily seen in the Hamilton–Jacobi description of classical mechanics (e.g. Bühler Reference Bühler2006; Sakurai & Napolitano Reference Sakurai and Napolitano2020).

The quantisation conditions derived above for separable problems, from which we obtained good approximations of the frequency shifts $\omega$, can be generalised to some extent to non-separable problems like the dipole flow (figure 2). This semi-classical analysis of a quantum system was developed by Einstein (Reference Einstein1917), Brillouin (Reference Brillouin1926) and Keller (Reference Keller1958), extending the Bohr–Sommerfeld quantum theory. The resulting approach is referred to as the Einstein–Brillouin–Keller (EBK) method (see also Berry & Mount Reference Berry and Mount1972; Percival Reference Percival1977; Keller Reference Keller1985). The starting point is that the rays (classical trajectories in the quantum problem), being constrained by the invariant $\omega$ (energy in the quantum problem), trace out invariant tori in the phase space spanned by $\boldsymbol {x}$ and $\boldsymbol {k}$. A ray starting on such a torus will remain on it forever. The quantisation condition selects invariant tori that correspond to allowed bound states by insisting that phase increments along closed loops on the invariant torus integrate to multiples of $2{\rm \pi}$. Recalling that $\varepsilon ^{-2} S_0 = {\rm i} \boldsymbol {k}$, so $\boldsymbol {k}$ is the spatial gradient of the phase and $\boldsymbol {k} \boldsymbol {\cdot } \mathrm {d}\kern0.7pt\boldsymbol {x}$ is a phase increment, the quantisation conditions read

(4.27a,b)\begin{equation} \oint_{\mathcal{C}_1}\boldsymbol{k}\boldsymbol{\cdot} \, \mathrm{d}\kern0.7pt\boldsymbol{x} = 2{\rm \pi}\left(n+\frac{1}{2}\right), \quad \oint_{\mathcal{C}_2}\boldsymbol{k}\boldsymbol{\cdot} \, \mathrm{d}\kern0.7pt\boldsymbol{x} = 2{\rm \pi} m, \end{equation}

where $n$ and $m$ are integers. The contours $\mathcal {C}_1$ and $\mathcal {C}_2$ are topologically independent closed curves on the invariant torus (figure 9a). In our example, the curve $\mathcal {C}_1$ passes through the hole of the phase space torus, whereas the curve $\mathcal {C}_2$ goes around the hole. The two curves are independent in the sense that neither one can be continuously deformed into the other. There is a half-integer phase shift in the quantisation condition arising from the integral along the curve $\mathcal {C}_1$ because this curve passes through two caustics, the generalisation of a turning point, where additional phase shifts are incurred (Brillouin Reference Brillouin1926; Keller Reference Keller1958; Maslov Reference Maslov1972). The curve $\mathcal {C}_2$ encounters no caustics. The integer wavenumbers $n$ and $m$ correspond to the cross- and along-streamline variations, respectively. These EBK quantisation conditions are entirely analogous to the WKB quantisation conditions derived above for the separable parallel shear flow and axisymmetric flow.

Figure 9. (a) Example of a trajectory tracing out an invariant torus for the dipole case. This torus corresponds to $n=2$, $m=0$ for $\varepsilon = \frac {1}{4}$. The background colours show the vorticity field. The black line shows a finite-time trajectory on the torus. The green and magenta lines represent a choice for the two invariant curves on the torus. They are independent because no continuous deformation of one can transform it into the other. (b) Different initial conditions result in different trajectories. This example is not bound to an invariant torus but is instead an example of a chaotic trajectory.

We apply the EBK quantisation to the dipole flow with $\varepsilon = \frac {1}{4}$. Our procedure closely follows Percival & Pomphrey (Reference Percival and Pomphrey1976): we find the invariant tori satisfying the quantisation condition by writing the Hamiltonian equations in action–angle variables and employing Newton's method; see Appendix D for details. All eigenvalues calculated by this EBK method show excellent agreement with the numerical values (figure 10).

Figure 10. Numerical eigenvalues (black circles) and EBK eigenvalues (orange crosses) calculated for the dipole flow with $\varepsilon = \frac {1}{4}$. Here EBK calculations are only shown for the sufficiently confined modes where the invariant tori are easy to calculate. The EBK values agree with the numerical values to $O(10^{-3})$.

As foretold by Einstein (Reference Einstein1917), not all modes are accessible by the EBK approach. If the system is non-integrable, trajectories in phase space can become chaotic instead of tracing out an invariant torus (figure 9b). States corresponding to such chaotic trajectories are not amenable to the EBK method. This ‘quantum chaos’ has received much attention in the physics literature and has connections to random matrix theory (e.g. Gutzwiller Reference Gutzwiller1992; Stone Reference Stone2005). Methods exist to estimate eigenvalues as well as their statistics (e.g. Edelman & Rao Reference Edelman and Rao2005; Edelman & Sutton Reference Edelman and Sutton2007). We do not pursue these issues any further here, in part because a uniform forcing projects most strongly onto the regular modes accessible with the EBK method (figure 8).