2.1. Wave dynamics in a shallow upper layer

Consider the set-up sketched in figure 1, namely a two-layer model with upper-layer density $\rho _1$ and lower-layer density $\rho _2\gt \rho _1$ , in a frame rotating at a positive rate $f/2$ around the vertical axis. Denoting time as $\tau$ , the upper layer has depth $h(x,y,\tau )$ , with $h=H_0$ in the rest state. The lower layer is infinitely deep. The base state consists of a steady, vertically invariant background flow $\boldsymbol {U}_{\!g}(x,y)$ spanning both layers. The background flow is in geostrophic balance with a vertically invariant lateral pressure gradient. It stems from a streamfunction $\psi (x,y)$ , that is, $\boldsymbol { U}_{\!g}=[U_{\!g}(x,y),V_{\!g}(x,y),0]=-\boldsymbol{\nabla } \times (\psi \, \boldsymbol {e}_z)$ . There is no deformation of the interface associated with such a vertically invariant balanced flow. Following DVB, we consider the rotating shallow water equations in the upper layer, linearised around the background balanced flow:

(2.1) \begin{align} \partial _\tau u +J(\psi ,u) +u \partial _x U_{\!g} + v \partial _y U_{\!g} - f v & = -g' \partial _x h , \\[-9pt] \nonumber \end{align}

(2.2) \begin{align} \partial _\tau v + J(\psi ,v) + u \partial _x V_{\!g} + v \partial _y V_{\!g} + f u & = -g' \partial _y h , \\[-9pt] \nonumber \end{align}

(2.3) \begin{align} \partial _\tau h + J(\psi ,h) + H_0 \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol {u} & = 0 , \\[12pt] \nonumber \end{align}

where $\boldsymbol {u}(x,y,\tau )=(u,v)$ denotes the horizontal velocity in the upper layer, $\boldsymbol{\nabla }=(\partial _x,\partial _y)$ , the Jacobian operator is $J(s,q)=(\partial _x s) (\partial _y q) - (\partial _x q) (\partial _y s)$ and the reduced gravity is $g'=g\,(\rho _2-\rho _1)/\rho _1$ with $g$ the acceleration of gravity.

Figure 1. A two-layer model with an infinitely deep lower layer. The base state consists of a vertically invariant steady horizontal flow $\boldsymbol {U}_{\!g}(x,y)$ spanning both layers, together with a flat interface between the two layers. We consider perturbations $\boldsymbol {u}(x,y,t)$ to the horizontal velocity in the upper layer only, whose depth is then denoted as $h(x,y,t)$ . In line with the rigid-lid approximation, we neglect the fluctuations of the free surface as compared with $h$ .

In the absence of background flow, $U_{\!g}=V_{\!g}=0$ , (2.1)–(2.3) support interfacial waves of frequency $\omega =f \sqrt {1 + k^2 \lambda ^2}$ , where $k$ denotes the wavenumber and $\lambda =\sqrt {g' H_0}/f$ is the small Rossby deformation radius associated with the shallow mixed layer. NIWs correspond to interfacial waves with wavelength much greater than $\lambda$ , their frequency $\omega \simeq f (1 + k^2 \lambda ^2/2)$ being close to $f$ . The large-scale atmospheric forcing induces NIW in the upper ocean, and these waves remain near-inertial because the background geostrophic flow has a typical scale $L_\psi \gg \lambda$ .

We non-dimensionalise (2.1)–(2.3) in such a way that the dimensionless fields and variables are ${\mathcal{O}}(1)$ at leading order in the expansion to come. Time is non-dimensionalised with $1/f$ and horizontal scales with $L_\psi$ . The background-flow streamfunction is non-dimensionalised using its root-mean-square (r.m.s.) value $\psi _{\textit{rms}}$ , where the mean is performed over space. Denoting as $U_w$ the infinitesimal velocity scale of the wave field, the wavy displacement of the interface scales as $H_0 U_w/(f L_\psi )$ . With such scalings, the dimensionless fields and variables read

(2.4) \begin{align} & \tau =\frac {\tilde {\tau }}{f} \, , \quad {\boldsymbol{x}} = L_{\psi } \tilde {{\boldsymbol{x}}} \, , \quad (u,v)=U_w (\tilde {u},\tilde {v}) \, , \quad h = H_0\left ( 1 + \frac {U_w}{f L_{\psi }}\tilde {h} \right ) \! , \\[-12pt] \nonumber \end{align}

(2.5) \begin{align} & \psi =\psi _{\textit{rms}} \chi \, , \quad \boldsymbol {U}_{\!g}=\frac {\psi _{\textit{rms}}}{L_\psi } \tilde {\boldsymbol {U}} \quad \text{with} \quad \tilde {\boldsymbol {U}}=-\tilde {\boldsymbol{\nabla }} \times (\chi \boldsymbol {e}_z) , \\[9pt] \nonumber \end{align}

where tildes denote dimensionless quantities and derivatives with respect to dimensionless variables, and $\chi (x,y)$ is the dimensionless streamfunction. Denoting space average with angular brackets, the latter satisfies $ \langle \chi ^2 \rangle =1$ . Substituting (2.4) into (2.1)–(2.3) leads to the dimensionless equations

(2.6) \begin{align} \partial _{\tilde {\tau }} \tilde {u} +\textit{Ro}_\psi \big[\tilde {J}(\chi ,\tilde {u}) -\tilde {u} \chi _{\tilde {x}\tilde {y}} - \tilde {v} \chi _{\tilde {y}\tilde {y}} \big] - \tilde {v} & = - \epsilon \partial _{\tilde {x}} \tilde {h} , \\[-9pt] \nonumber \end{align}

(2.7) \begin{align} \partial _{\tilde {\tau }} \tilde {v} + \textit{Ro}_\psi \big[\tilde {J}(\chi ,\tilde {v}) + \tilde {u} \chi _{\tilde {x}\tilde {x}} + \tilde {v} \chi _{\tilde {x}\tilde {y}} \big] + \tilde {u} & = - \epsilon \partial _{\tilde {y}} \tilde {h} , \\[-9pt] \nonumber \end{align}

(2.8) \begin{align} \partial _{\tilde {\tau }} \tilde {h} + \textit{Ro}_\psi \tilde {J}(\chi ,\tilde {h}) + \boldsymbol{\tilde {\boldsymbol{\nabla }}} \boldsymbol{\cdot }{\tilde {\boldsymbol {u}}} & = 0 , \\[9pt] \nonumber \end{align}

where $\epsilon =(\lambda /L_\psi )^2$ denotes the Burger number and $\textit{Ro}_\psi =\psi _{\textit{rms}}/(f L_\psi ^2)$ denotes the Rossby number of the background flow.

2.2. The Young–Ben Jelloul (YBJ) equation

Young & Ben Jelloul (Reference Young and Ben Jelloul1997) consider the distinguished asymptotic regime $\epsilon \ll 1$ , $\textit{Ro}_\psi \ll 1$ , keeping a finite ratio $\gamma =\textit{Ro}_\psi /\epsilon ={\mathcal{O}}(\epsilon ^0)$ . Through an asymptotic expansion recalled in Appendix A, YBJ show that the horizontal velocity field $(u,v)$ consists of inertial oscillations whose complex amplitude slowly varies with time. They derive a reduced equation governing the modulation of this complex amplitude as a result of advection and refraction by the weak background flow, together with wave dispersion. For the present set-up, this procedure results in the following evolution equation for the demodulated complex velocity field $M({\boldsymbol{x}},t)=(u+iv)e^{i \tau }$ :

(2.9) \begin{align} \partial _t M + \underbrace {\vphantom {\frac {i}{2}} \gamma J(\chi , M)}_{\text{advection}} + \underbrace {\frac {i \gamma }{2}\left (\Delta \chi \right ) M}_{\text{refraction}} - \underbrace {\frac {i}{2}\Delta M}_{\text{dispersion}}= 0 , \end{align}

where $\Delta = \partial _{xx} + \partial _{yy}$ is the Laplace operator and $t=\epsilon \tau$ is a slow time variable (dropping the tildes to alleviate notation). Equation (2.9) is the simplest instance of the YBJ model. It involves the single dimensionless parameter $\gamma \geqslant 0$ characterising the strength of the background flow relative to wave dispersion. In terms of dimensional variables, the expression of $\gamma$ is

(2.10) \begin{align} \gamma = \frac {\psi _{\textit{rms}} f}{g' H_0} . \end{align}

Pure inertial oscillations with $M=\text{const.}$ are valid solutions to the YBJ equation (2.9) in the absence of background flow only, that is, for $\gamma =0$ . For $\gamma \neq 0$ , a uniform initial condition for $M$ evolves with time, developing some spatial structure as a result of refraction and advection by the background flow, and wave dispersion. We are interested in the fate of a uniform NIW field induced by a large-scale atmospheric storm. Equation (2.9) being linear and invariant to a uniform phase shift of the complex variable $M$ , we focus on the initial condition $M({\boldsymbol{x}},t=0)=1$ in the following.

3.1. Particle in a steady electromagnetic field

Consider a particle of mass $m$ and positive charge $q$ in a steady electromagnetic field whose potentials depend on $x$ and $y$ . Using a set of units such that $m=q=\hbar =1$ (that is, using a non-dimensionalisation based on $q$ , $m$ and $\hbar$ ), the dimensionless Hamiltonian reads

(3.2) \begin{align} H({\boldsymbol{x}},\boldsymbol {p}) = \frac {1}{2}[\boldsymbol {p}-\boldsymbol {A}(x,y)]^2 +V(x,y) , \end{align}

where $\boldsymbol {A}(x,y)$ denotes the vector potential and $V(x,y)$ denotes the electrostatic potential, both dimensionless. The quantum dynamics of the particle are governed by the Schrödinger equation, $i \partial _t \phi = H\{ \phi \}$ , where $\phi (x,y,t)$ denotes the wave function and the momentum $\boldsymbol {p}$ in the Hamiltonian (3.2) is replaced by the operator $-i \boldsymbol{\nabla }$ . Now, upon choosing dimensionless potentials that are related to the streamfunction of the NIW problem through

(3.3) \begin{align} \boldsymbol {A}= \gamma \boldsymbol{\nabla } \times (\chi \boldsymbol {e}_z) \quad \text{and} \quad V=\frac {1}{2} \big (\gamma {\Delta \chi }- \gamma ^2 {|\boldsymbol{\nabla } \chi |^2} \big ) , \end{align}

the Schrödinger equation for the wave function $\phi (x,y,t)$ reduces precisely to the YBJ equation (3.1) for the demodulated velocity $M(x,y,t)$ . We conclude that there is an exact analogy between the YBJ equation and the quantum dynamics of a charged particle in the electromagnetic field given by the potentials (3.3), the demodulated velocity $M(x,y,t)$ playing the role of the wave function of the charged particle. Within this analogy, the vector potential $\boldsymbol {A}$ equals minus the background flow velocity (therefore, $\boldsymbol {A}$ satisfies the Coulomb’s gauge condition $\boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol {A}=0$ ), while the scalar potential $V$ equals one half the background flow vorticity, minus the background flow kinetic energy. Table 1 further lists analogous quantities between the two systems.

Table 1. Summary of the analogy between the Schrödinger equation for a charged particle (left column) and the YBJ equation (right column).

3.2. Conserved quantities

There are two ways of determining the conserved quantities of the YBJ equation. One can directly deduce them from the equation or one can readily infer them from the quantum analogy. Consider the YBJ equation inside a doubly periodic domain $(x,y)\in {\mathcal D}=[0,1]^2$ . Multiplying the YBJ equation (2.9) with $M^*$ before adding the complex conjugate and averaging over the domain $\mathcal D$ yields, after a few integrations by parts using the periodic boundary conditions:

(3.4) \begin{align} \frac {\mathrm{d} {\mathcal A}}{\mathrm{d}t}=0 \quad \text{with} \quad {\mathcal A}=\big \langle |M|^2 \big \rangle , \end{align}

where the angular brackets denote space average over the domain $\mathcal D$ . Alternatively, the conservation of $\mathcal A$ is readily inferred from the quantum analogy, as $\mathcal A$ corresponds to the conserved total probability of finding the particle somewhere inside the domain $\mathcal D$ . In the YBJ context, $\mathcal A$ corresponds to wave action, usually defined as the ratio of the wave energy to the wave frequency. The mechanical energy of NIWs is dominated by the kinetic energy $ \langle |M|^2 \rangle$ (omitting the prefactor $1/2$ ), while the frequency is equal to $f$ to lowest order. The conservation of wave action thus reduces to the conservation of the space-averaged kinetic energy $ \langle |M|^2 \rangle$ of the wave field.

The conservation of $\mathcal A$ is discussed in the original YBJ paper (Young & Ben Jelloul Reference Young and Ben Jelloul1997). Eighteen years later, a second independent conserved quantity was uncovered by DVB based on manipulations of the YBJ equation. Once again, this second invariant is readily inferred from the quantum analogy. Indeed, the Hamiltonian being time-independent, its expectation value is conserved over time: the mechanical energy of the charged particle is conserved. In the quantum context, this expectation value is $\left \langle \phi | H |\phi \right \rangle$ . In the YBJ context, this quantity becomes $\int _{\mathcal D} M^* H\{ M \}\,\mathrm{d}{\boldsymbol{x}}$ , where $H\{ M \}$ is given by the right-hand side of (3.1). The conserved quantity finally reads

(3.5) \begin{align} \nonumber E & = \left \langle M^* \left [ -\frac {1}{2} \Delta M -i \gamma J(\chi ,M) + \frac {\gamma }{2} (\Delta \chi )M \right ] \right \rangle \nonumber \\ & = \left \langle \underbrace {\frac {|\boldsymbol{\nabla } M|^2}{2}}_{\text{potential}} + \underbrace {\frac {\gamma \Delta \chi }{2} |M|^2 \vphantom {\frac {|\boldsymbol{\nabla } M|^2}{2}}}_{\text{refraction}} + \underbrace {i \gamma \chi J(M^*,M) \vphantom {\frac {|\boldsymbol{\nabla } M|^2}{2}}}_{\text{advection}} \right \rangle \!, \end{align}

where we have performed various integrations by parts using the periodic boundary conditions to obtain the second expression. We refer to (3.5) as the wave energy. The various terms on the right-hand side of (3.5) correspond to the potential energy, the contribution from the refractive term and the contribution from the advective term. In addition, the equation $\mathrm{d}E/\mathrm{d}t=0$ can be recast as an evolution equation for a single one of these energy terms, provided one substitutes the YBJ expression for $\partial _t M$ in the time derivative of the other forms of energy, see Rocha et al. (Reference Rocha, Wagner and Young2018).

Strictly speaking, the total mechanical energy of the waves consists of a leading-order kinetic energy term, proportional to $\mathcal A$ , and the weaker contributions gathered in $E$ mentioned above. In the present context, $\mathcal A$ and $E$ are conserved independently. In the absence of background flow, $\gamma =0$ , only the potential energy contribution $ \langle |\boldsymbol{\nabla } M|^2/2 \rangle$ remains in (3.5).

4.1. Numerical set-up

The goal of the present study is to characterise the organisation of the NIW field over a steady background flow, comparing the theoretical predictions to numerical simulations of the YBJ equation (2.9) in the doubly periodic domain $\mathcal D$ . The simulations are performed using standard pseudo-spectral methods on a GPU with dealiasing and RK4 time stepping. The time step is fixed for a given simulation. No hyperviscosity is required, as the spatial spectrum of $M$ naturally exhibits an emergent cutoff wavenumber within the resolved scales. The parameter values of all numerical simulations are reported in Appendix F. The initial condition is $M({\boldsymbol{x}},t=0)=1$ . For the steady background velocity field entering the equation, we use an instantaneous velocity field extracted from a simulation of the two-dimensional (2-D) Navier–Stokes (NS) equations, following the same forcing and dissipation protocols as described by Meunier & Gallet (Reference Meunier and Gallet2025), albeit in a different parameter regime. We low-pass filter this frozen-in-time velocity field to remove excessively small scales with wavenumber $k \gtrsim 150$ . We display the streamfunction, kinetic energy and vorticity of the resulting background flow in figure 2. We stress the fact that this flow is time-independent in the YBJ equation.

Figure 2. Steady background flow used in the numerical simulations of the YBJ equation: (a) background streamfunction $\chi (x,y)$ ; (b) kinetic energy $|\boldsymbol{\nabla } \chi |^2$ and (c) vorticity field $\Delta \chi$ . The normalisation is such that $\langle \chi ^2 \rangle =1$ (see text). In all panels, the black contours correspond to streamlines of the background flow.

4.2. Quantities of interest

In the following, we mainly discuss the time-averaged spatial distributions of various forms of energy in the system. Denoting time-average with an overbar, we consider the spatial distribution of wave action $\overline {|M|^2}(\boldsymbol {x})$ , which also corresponds to the spatial distribution wave kinetic energy. We also consider the spatial distribution of the mean squared gradient of $M$ , $\overline {|\boldsymbol{\nabla } M|^2}(\boldsymbol {x})$ . The latter being the dominant contribution to the NIW potential energy in both limits $\gamma \ll 1$ and $\gamma \gg 1$ , we simply refer to it as the NIW potential energy in the following. In the two-way coupled model derived by Xie & Vanneste (Reference Xie and Vanneste2015), $|\boldsymbol{\nabla } M|^2$ represents the NIW contribution to the total energy invariant, which makes it a quantity of interest for predictions. Together with these various forms of energy, we also discuss the time-averaged Stokes drift $\overline {\boldsymbol {u}_s}({\boldsymbol{x}})$ induced by the wave field. The precise definition of the Stokes drift is deferred to Appendix B, where we show that the time-averaged Stokes drift is related to the complex amplitude $M$ entering the YBJ equation through

(4.1) \begin{align} \overline {\boldsymbol {u}_s}({\boldsymbol{x}}) = \frac {1}{4} \big [ i(\overline {M \boldsymbol{\nabla } M^*} - \overline {M^* \boldsymbol{\nabla } M}) + \boldsymbol{\nabla } \times (\overline {|M|^2} \boldsymbol {e}_z) \big ]. \end{align}

The dimensionless Stokes drift appearing in (4.1) corresponds to the dimensional Stokes drift divided by $U_w^2/(f L_\psi )$ .

4.3. Two limiting regimes

Once the flow structure $\chi (x,y)$ and the initial condition $M=1$ have been fixed, the only dimensionless parameter entering the problem is the strength $\gamma$ of the background flow. Guided by the quantum analogy, in the following, we focus on two limiting situations of interest.

1. (i) $\gamma \ll 1$ : this is the ‘quantum’ or ‘strong-dispersion’ limit. The background flow is weak and the dispersive effects in the YBJ equation (2.9) are strong.

2. (ii) $\gamma \gg 1$ : this is the limit of ‘classical mechanics’. The YBJ equation is analogous to the dynamics of a quantum particle in the small- $\hbar$ limit.

6.1. Ergodic theory and microcanonical ensemble

Instead of Newton’s third law, the equilibrium statistical mechanics of a classical system starts from the classical version of the Hamiltonian (3.2). Hamilton’s equations govern the evolution of a narrow wave packet located at ${\boldsymbol{x}}(t)=[x(t),y(t)]$ with wavevector $\boldsymbol { p}(t)=[p_x(t),p_y(t)]$ (see sketch in figure 4):

(6.3) \begin{align} \frac {\mathrm{d} {\boldsymbol{x}}}{\mathrm{d}t}= \frac {\partial H}{\partial \boldsymbol {p}} \, , \qquad \frac {\mathrm{d} \boldsymbol {p}}{\mathrm{d}t}= - \frac {\partial H}{\partial {\boldsymbol{x}}} , \end{align}

where the equations are to be understood componentwise. Consider a cloud of initial conditions in the phase space $(x,y,p_x,p_y)$ . Liouville’s theorem states that, following the Hamiltonian evolution equations (6.3), the cloud will deform in phase space conserving its initial volume. In other words, the volume in phase space is conserved by the dynamics because (6.3) correspond to an incompressible flow in phase space.

Figure 4. A narrow wave packet with mean position ${\boldsymbol{x}}(t)$ and wavevector $\boldsymbol {p}(t)$ behaves like a charged classical particle in a steady 2-D electromagnetic field.

Consider now an ensemble of particles with the same initial energy $E_0$ . Because energy is conserved, these particles only have access to the hypersurface $H({\boldsymbol{x}},\boldsymbol {p})=E_0$ in phase space. Like a patch of dye getting homogenised by a chaotic flow and achieving uniform concentration in the long-time limit, we expect the Hamiltonian phase-space flow (6.3) to homogenise a cloud of initial conditions with initial energy $E_0$ over the hypersurface $H({\boldsymbol{x}},\boldsymbol { p})=E_0$ . Introducing a probability density ${\mathcal P}({\boldsymbol{x}},\boldsymbol {p})$ such that ${\mathcal P}({\boldsymbol{x}},\boldsymbol { p})\,\mathrm{d}{\boldsymbol{x}} \,\mathrm{d}\boldsymbol {p}$ is the probability for a particle to be in a phase-space volume $\mathrm{d}{\boldsymbol{x}}\, \mathrm{d}\boldsymbol {p}$ around the point $({\boldsymbol{x}},\boldsymbol {p})$ , this ergodic assumption translates into

(6.4) \begin{align} {\mathcal P}({\boldsymbol{x}},\boldsymbol {p}) = {\mathcal C} \, \delta [H({\boldsymbol{x}},\boldsymbol {p})-E_0] , \end{align}

where the constant $\mathcal C$ is a normalisation factor. The validity of the ergodic prescription (6.4) is a lively topic of research at the crossroads of mathematics and physics, known as ‘quantum chaos’ (Berry Reference Berry1977). Rigorous mathematical results have been obtained based on the asymptotic behaviour of the Wigner transform of the wavefunction in the classical limit (Voros Reference Voros1976). A detailed discussion of this topic goes beyond the scope of the present study. Instead, we motivate (6.4) based on the microcanonical ensemble prescription of equilibrium statistical mechanics (Bouchet & Venaille Reference Bouchet and Venaille2012), which is expected to correctly describe the statistics of the quantum system in the classical limit $\gamma \gg 1$ . In the following, we thus assume that the ergodic assumption holds and we replace long-time averages with averages in phase space using the probability density (6.4).

6.2. Distribution of NIW kinetic energy

As a first illustration, let us determine the time-averaged distribution of NIW kinetic energy (or wave action) $\overline {|M|^2}({\boldsymbol{x}})$ using an average in phase space. According to table 1, $\overline {|M|^2}({\boldsymbol{x}})$ corresponds to the time-averaged probability of finding the quantum particle at location $\boldsymbol{x}$ . Additionally, in the classical limit, this reduces to the probability of finding the classical particle at location $\boldsymbol{x}$ regardless of its momentum $\boldsymbol { p}$ . Using the microcanonical probability density (6.4), the latter probability is given by

(6.5) \begin{align} \overline {|M|^2}({\boldsymbol{x}}) & = \int _{{\boldsymbol{x}}'\in {\mathcal D} , \, \boldsymbol {p}\in \mathbb{R}^2} \underbrace {\delta ({\boldsymbol{x}}'-{\boldsymbol{x}})}_{\text{particle located at ${\boldsymbol{x}}$}} \, { {\mathcal P}}({\boldsymbol{x}}',\boldsymbol {p}) \, \mathrm{d}{\boldsymbol{x}}' \mathrm{d}\boldsymbol {p} \\[-12pt] \nonumber \end{align}

(6.6) \begin{align} & = {\mathcal C} \int _{\boldsymbol {p}\in \mathbb{R}^2} \,\delta \left [ \frac {1}{2}(\boldsymbol {p}+\gamma \boldsymbol { U}({\boldsymbol{x}}))^2 + V({\boldsymbol{x}}) -E_0 \right ] \, \mathrm{d}\boldsymbol {p} . \\[9pt] \nonumber \end{align}

Changing integration variable to $\boldsymbol {K}=\boldsymbol {p}+\gamma \boldsymbol {U}({\boldsymbol{x}})$ with norm $K=|\boldsymbol { K}|$ , the integral becomes

(6.7) \begin{align} \overline {|M|^2}({\boldsymbol{x}}) & = {\mathcal C} \int _{K \in \mathbb{R}^+} \,\delta \left [ \frac {1}{2}K^2 + V({\boldsymbol{x}}) -E_0 \right ] 2 \pi\! K \, \mathrm{d}K \\[-12pt] \nonumber \end{align}

(6.8) \begin{align} & = 2 \pi {\mathcal C} \int _0^\infty \,\delta \left [ s + V({\boldsymbol{x}}) -E_0 \right ] \, \mathrm{d}s , \\[9pt] \nonumber \end{align}

where we changed the integration variable again to $s=K^2/2$ . The resulting integral equals one if $V({\boldsymbol{x}}) -E_0\lt 0$ and zero if $V({\boldsymbol{x}}) -E_0\gt 0$ , that is,

(6.9) \begin{align} \overline {|M|^2}({\boldsymbol{x}}) & = 2 \pi {\mathcal C} \, {\mathcal H}[E_0-V({\boldsymbol{x}})] , \end{align}

where $\mathcal H$ denotes the Heavyside function.

The initial energy of the particles is estimated by inserting the initial condition $M(x,y,t=0)=1$ into (3.5) for the energy. Only the term $\gamma (\Delta \chi )|M|^2/2$ remains: the local initial energy is of the order of the local vorticity and therefore it scales as $\gamma$ . In contrast, the potential (6.1) has much greater magnitude, of order $\gamma ^2$ , and it is negative almost everywhere. We conclude that the initial energy is negligible as compared with the potential $V\lt 0$ in the limit $\gamma \gg 1$ of interest here: $E_0 \simeq 0$ (see Appendix C.3 for a refined calculation taking into account the narrow distribution of $E_0$ around zero). To a good approximation, ${\mathcal H}[E_0-V({\boldsymbol{x}})]=1$ almost everywhere and we thus predict a uniform distribution of NIW kinetic energy, $\overline {|M|^2}({\boldsymbol{x}}) = 2 \pi {\mathcal C}$ . Because of action conservation, the space average of $|M|^2$ is conserved and equal to one. We thus obtain ${\mathcal C}=1/(2 \pi )$ , the prediction for the time-averaged spatial distribution of kinetic energy being simply

(6.10) \begin{align} \overline {|M|^2}({\boldsymbol{x}}) & = 1 . \end{align}

Somewhat surprisingly, based on statistical mechanics, we predict a uniform distribution of NIW kinetic energy, despite the spatial structure of the potential (6.1). This long-time behaviour contrasts with the early-time behaviour of the system, as described e.g. by DVB and by Rocha et al. (Reference Rocha, Wagner and Young2018). At early time, the uniform initial condition for $M$ is affected predominantly by the refraction term, which induces strong phase gradients driving an action flux towards the centre of anticyclones. For subsequent times, however, the advection term comes into play and, for $\gamma \gg 1$ , DVB show that the dominant balance in the YBJ equation is eventually between advection and dispersion. Similarly, one can check that the refraction term plays a negligible role in the line of arguments leading to the uniform prediction (6.10) from ergodic theory. In § 6.6, we address the influence of the small refraction term in more details to refine the prediction (6.10).

6.3. Distribution of NIW potential energy

As a second illustration of the statistical mechanics approach, we consider the time-averaged spatial distribution of NIW potential energy, $\overline {|\boldsymbol{\nabla } M|^2}({\boldsymbol{x}})$ . The dimensionless momentum operator being $-i \boldsymbol{\nabla }$ according to the quantum analogy, the NIW potential energy is analogous to the expectation value of the squared momentum. Alternatively, based on the sketch in figure 4, one estimates $\boldsymbol{\nabla } M \simeq i \boldsymbol {p} M$ and $|\boldsymbol{\nabla } M|^2 \simeq p^2 |M|^2$ , in line with the standard WKB approximation. We thus seek the averaged squared momentum of the particles located at $\boldsymbol{x}$ . In phase space, this average reads

(6.11) \begin{align} \overline {|\boldsymbol{\nabla } M|^2}({\boldsymbol{x}}) & = \int _{{\boldsymbol{x}}'\in {\mathcal D}, \, \boldsymbol {p}\in \mathbb{R}^2} \underbrace {p^2}_{\text{squared momentum}} \delta ({\boldsymbol{x}}'-{\boldsymbol{x}}) \, {\mathcal P}({\boldsymbol{x}}',\boldsymbol {p}) \, \mathrm{d}{\boldsymbol{x}}' \mathrm{d} \boldsymbol {p}, \\[-12pt] \nonumber \end{align}

(6.12) \begin{align} & = 2 \gamma ^2 \boldsymbol {U}({\boldsymbol{x}})^2 , \\[9pt] \nonumber \end{align}

the integration in phase space being detailed in Appendix C. We conclude that the spatial distribution of NIW potential energy is given by the kinetic energy distribution of the background flow, with an accumulation of NIW potential energy in fast-flow regions.

6.4. Time-averaged Stokes drift

As a last illustration of the statistical mechanics approach, we consider the time-averaged Stokes drift (4.1). In the $\gamma \gg 1$ limit, substitution of the prediction (6.10) for $\overline {|M|^2}({\boldsymbol{x}})$ into (4.1) shows that the second term vanishes. Inserting again the estimate $\boldsymbol{\nabla } M \simeq i \boldsymbol {p} M$ , the time-averaged Stokes drift (4.1) reduces to

(6.13) \begin{align} \overline {\boldsymbol {u}_s}({\boldsymbol{x}}) = \frac {i}{4} (\overline {M \boldsymbol{\nabla } M^*} - \overline {M^* \boldsymbol{\nabla } M}) \simeq \frac {1}{2} \overline {\boldsymbol {p} |M|^2}({\boldsymbol{x}}) . \end{align}

Once again, we assume ergodicity to compute the average appearing on the right-hand side in phase space. The integration in phase space is deferred to Appendix C, the end result being

(6.14) \begin{align} \overline {\boldsymbol {u}_s}({\boldsymbol{x}}) =- \frac {\gamma }{2} \boldsymbol {U}({\boldsymbol{x}}) . \end{align}

It is remarkable that we obtain the same prediction for the time-averaged Stokes drift in the strong-dispersion limit and in the strong-advection limit, see (5.9) and (6.14), although these two predictions arise from different terms in the expression (4.1) of the Stokes drift.

6.5. Comparison to numerical simulations

In figure 5, we compare the time-averaged spatial distributions of potential energy and Stokes drift with the ergodic predictions (6.12) and (6.14). The agreement is very good in both cases, both at the qualitative and at the quantitative levels, showing that NIW packets indeed tend to behave in an ergodic fashion in the strong-advection limit. While the large-scale structure of both fields is accurately captured by the ergodic theory, one can notice some small-scale roughness or ‘granularity’ in the numerical fields in figure 5. This phenomenon does not seem to disappear over longer time average. Rather, it stems from some form of interference pattern, reminiscent of the patterns obtained in studies of quantum chaos (Voros Reference Voros1976; Berry Reference Berry1977; Nonnenmacher Reference Nonnenmacher2013). The latter studies suggest that the granularity arises at the de Broglie wavelength of the system and that classical statistics apply to quantities averaged over a few de Broglie wavelengths. Equations (6.13) and (6.14) point to the scaling $p\sim \gamma$ for the typical momentum of the particle, and, therefore, a de Broglie wavelength that scales as $1/p \sim 1/\gamma$ (remembering that $\hbar =1$ with our choice of units). The same estimate is readily obtained by DVB who show that, in the large- $\gamma$ regime, the characteristic scale of $M$ is set through a balance between advection and dispersion. Some slight and slow time dependence in the position of the vortices of the background flow, or an ensemble average over a family of slightly structured initial conditions for $M$ , may be sufficient to disrupt the precise phase relations producing the interference pattern, thus smoothing out the observed granularity.

Figure 5. Ergodic predictions for the (a) time-averaged NIW potential energy $\overline {|\boldsymbol{\nabla } M|^2}({\boldsymbol{x}})$ and (b) Stokes’ drift $\overline {\boldsymbol{u}_s}({\boldsymbol{x}})$ . (c,d) Same fields as panels (a) and (b) extracted from a numerical run with $\gamma =30$ . In panels (a) and (c), black contours indicate isovalues 1, 4, 9 and 16.

Consider now the spatial distribution of NIW kinetic energy, illustrated in figure 6. The naive ergodic prediction (6.10) corresponds to a uniform distribution of wave kinetic energy. There is reasonable agreement with this prediction: for large $\gamma$ , the fields in figure 6 feature an extended uniform white region with $\overline {|M|^2}({\boldsymbol{x}}) \simeq 1$ . However, locally, we observe some strong departures from this uniform background. These departures are located near the points of vanishing kinetic energy of the background flow, as illustrated by the contour of low $|\boldsymbol { U}|$ . More precisely, we distinguish between two types of regions.

1. (i) In regions of vanishing $\boldsymbol {U}^2({\boldsymbol{x}})$ with cyclonic background vorticity $\gamma \Delta \chi ({\boldsymbol{x}})\gt 0$ , we observe a deficit of NIW kinetic energy (blue regions).

2. (ii) In regions of vanishing $\boldsymbol {U}^2({\boldsymbol{x}})$ with anticyclonic background vorticity $\gamma \Delta \chi ({\boldsymbol{x}})\lt 0$ , we observe an excess of NIW kinetic energy, see the narrow red regions in the two strongest anticyclones in figure 6.

Figure 6. (a) Time-averaged NIW kinetic energy $\overline {|M|^2}$ from a numerical simulation with $\gamma =10$ . (b) Idem for $\gamma = 30$ . (c) Refined ergodic prediction (6.17) for $\gamma = 30$ , taking into account non-uniform initial energy and trapping in the two main anticyclones. The colourscale lightness varies linearly on $[0,1]$ and on $[1,5]$ to keep $\overline {|M|^2}=1$ in white. Black contours indicate the isovalue $1/3$ of the r.m.s. value of $|\boldsymbol{U}|$ . These contours show that deviations from the uniform distribution preferentially occur in regions of slow background flow.

6.6. Refined ergodic prediction and anticyclonic trapping

We now describe these regions in more detail, starting with the deficit regions. As mentioned previously, the deficit regions correspond to regions of vanishing kinetic energy of the background flow, together with positive vorticity. The full potential $V$ in (3.3) is positive in such regions. Few particles have sufficient initial energy to rise on top of such potential hills; hence, the deficit in NIW kinetic energy. In Appendix C.3, we derive the following refined statistical mechanics prediction for the distribution of NIW kinetic energy, taking into account the narrow distribution of the initial energy $E_0$ of the particles:

(6.15) \begin{align} \overline {|M|_{\textit{erg}}^2}({\boldsymbol{x}}) = \frac {1}{2} \left \{ 1- \text{erf} \left [\frac {\Delta \chi (\boldsymbol {x}) - \gamma |\boldsymbol{\nabla } \chi |^2 (\boldsymbol {x})}{\sqrt {2} (\Delta \chi )_{\textit{rms}}} \right ] \right \} \! , \end{align}

where $(\Delta \chi )_{\textit{rms}}$ denotes the r.m.s. vorticity of the background flow (r.m.s. value of $\Delta \chi$ ). As illustrated in figure 6, this refined prediction accurately captures the location and intensity of the deficit regions.

We now turn to the excess regions, which coincide with the vortical cores of the two main anticyclones of the background flow, whose centres are located at ${\boldsymbol{x}}_1=(0.761,0.830)$ and ${\boldsymbol{x}}_2=(0.663,0.107)$ . For large $\gamma$ , these excess regions occupy an arguably narrow fraction of the domain and they contain a tiny fraction of the overall NIW kinetic energy (in all the numerical simulations with $\gamma \ge 5$ , the overall surplus of energy contained in these two anticyclones is only approximately 0.5 % of the total energy). The NIW accumulation results from a breaking of ergodicity in the vortical cores of the anticyclones as NIW kinetic energy from the initial condition remains trapped there. Such trapping of NIW kinetic energy can be explained by the existence of localised eigenmodes in the vicinity of the anticyclonic vortex cores, see e.g. the eigenmodes computed by Llewellyn Smith (Reference Llewellyn Smith1999) for an axisymmetric vortex. To test this assumption, we ran a few simulations with a tailored initial condition that features very little NIW kinetic energy in the anticyclones. This initial condition prevents any trapping in the non-ergodic regions, and the regions of excess NIW kinetic energy are indeed absent from the resulting $\overline {|M|^2}({\boldsymbol{x}})$ (not shown).

To further describe NIW accumulation in the two anticyclones in the large- $\gamma$ limit, we propose an idealised model consisting of classical particles trapped in an axisymmetric vortex. We start from a uniform initial distribution of classical particles in the anticyclone, corresponding to the initial condition $M(\boldsymbol {x},t=0)=1$ . A consequence of axisymmetry is that the subsequent motion of the particles takes place in the radial direction only. In other words, the additional conservation of angular momentum prevents chaotic motion and ergodic statistics for the trapped particles, allowing only for oscillatory radial motion. Early evolution under the YBJ equation leads to $M({\boldsymbol{x}},t)\simeq 1-i\gamma \zeta t/2 + {\mathcal{O}}(t^2)$ , where $\zeta = \Delta \chi$ denotes the vorticity of the background flow (see also Asselin et al. Reference Asselin, Thomas, Young and Rainville2020). Denoting as $r$ the radial coordinate, the early-time radial velocity is thus $-i\partial _r M=-\gamma (\partial _r \zeta ) t /2$ . It points towards the centre of the anticyclone ( $\partial _r \zeta \gt 0$ ), indicating accumulation at early time. To estimate the radius of influence – or ‘trapping radius’ – of each anticyclone, we extract the distance between the vortex centre and the first point at which $-\boldsymbol{\nabla } \zeta$ points away from the vortex centre. This leads to the trapping radius $R_1 = 0.030$ (respectively $R_2 = 0.029$ ) for anticyclone 1 (respectively anticyclone 2).

In Appendix D, we describe the subsequent radial motion of the trapped classical particles within the (approximately) axisymmetric vortical core. First, we note that the NIW kinetic energy trapped in the two anticyclones represents only a tiny fraction of the total NIW kinetic energy of the system. This fraction is independent of $\gamma$ for large $\gamma$ , being equal to the initial kinetic energy contained within the disks of centres ${\boldsymbol{x}}_1$ and ${\boldsymbol{x}}_2$ , and radii $R_1$ and $R_2$ : $\pi (R_1^2+R_2^2) \simeq 0.005 \ll 1$ (the total NIW kinetic energy being equal to one with our non-dimensionalisation). Second, we compute the distribution of excess kinetic energy within the two disk-shaped regions, taking into account the vorticity profiles $\zeta _{1,2}(r)$ of the two anticyclones. We obtain

(6.16) \begin{align} \overline {|M|_{\textit{trap},i}^2}(r) = \frac {1}{r} \int _{r}^{R_i} \frac {r_0 \, \mathrm{d}r_0}{\sqrt {\zeta _i(r_0)-\zeta _i(r)} \times \int_{s=0}^{s=r_0} \dfrac {\mathrm{d}s}{\sqrt {\zeta _i(r_0)-\zeta _i(s)}}} \quad (i=1,2), \end{align}

with $r \leqslant R_i$ the distance to the centre of the anticyclone considered. The full prediction for the spatial distribution of NIW kinetic energy is obtained by adding the two excess distributions (6.16) to the ergodic contribution (6.15):

(6.17) \begin{align} \overline {|M|^2}({\boldsymbol{x}}) = \overline {|M|_{\textit{erg}}^2}({\boldsymbol{x}}) + \overline {|M|_{\textit{trap,1}}^2}(|{\boldsymbol{x}}-{\boldsymbol{x}}_1|) + \overline {|M|_{{trap,2}}^2}(|{\boldsymbol{x}}-{\boldsymbol{x}}_2|) . \end{align}

This prediction is displayed in figure 6. Beyond capturing the deficit regions, it reproduces the accumulation observed in the numerical simulations, at least qualitatively. In particular, the trapped kinetic energy gets redistributed approximately as $1/r$ , with $r$ the distance from the vortex centre. Hence, the narrow regions of intense $\overline {|M|^2}$ inside the two anticyclonic cores in figure 6.