<figure><img src="image_1.png" /></figure>

**Subject Areas:**

geophysics, fluid mechanics,
mathematical physics

**Keywords:**

waves, stratification, rotating flows,
triaxial ellipsoids

**Author for correspondence:**

Jérémie Vidal

e-mail:

jeremie.vidal@univ-grenoble-alpes.fr

<figure><img src="image_2.png" /></figure>

# Inertia-gravity waves in geophysical vortices

Jérémie Vidal<fnref n="1" /> & Yves Colin de Verdière<fnref n="2" />

<figure><img src="image_3.png" /></figure>

Pancake-like vortices are often generated by turbulence in geophysical flows. Here, we study the inertia-gravity oscillations that can exist within such geophysical vortices, due to the combined action of rotation and gravity. We consider a fluid enclosed within a triaxial ellipsoid, which is stratified in density with a constant Brunt-Väisälä frequency (using the Boussinesq approximation) and uniformly rotating along a (possibly) tilted axis with respect to gravity. The wave problem is then governed by a mixed hyperbolic-elliptic equation for the velocity. As in the rotating non-stratified case considered by Vantieghem (2014, Proc. R. Soc. A, **470**, 20140093, doi:10.1098/rspa.2014.0093), we find that the spectrum is pure point in ellipsoids (i.e. only consists of eigenvalues) with smooth polynomial eigenvectors. Then, we characterise the spectrum using numerical computations (obtained with a bespoke Galerkin method) and asymptotic spectral theory. Finally, the results are discussed in light of natural applications (e.g. for Mediterranean eddies or Jupiter's vortices).

## 1. Introduction

Geophysical flows are often influenced by the action of density stratification and rotation. In particular, buoyancy supports the propagation of internal (gravity) waves in stratified fluids [1]. Similarly, the Coriolis force sustains inertial waves in homogeneous rotating fluids [2]. Combining density stratification and global rotation then gives birth to a new wave family [3], usually called inertia-gravity waves (IGWs). Such various waves, which are ubiquitous in rotating and stratified environments, are often believed to be key for understanding the (turbulent) dynamics of geophysical flows [4,5].

<footnote n="1">¹ Université Grenoble Alpes, CNRS, ISTerre, 38000 Grenoble, France</footnote>

<footnote n="2">² Université Grenoble Alpes, CNRS, Institut Fourier, 38000 Grenoble, France</footnote>

<footnote>D JV, 0000-0002-3654-6633</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 1. (a) Cross-section of the temperature anomaly (acoustic tomography) in a Mediterranean eddy offshore in the Eastern Atlantic. Adapted from figure 3 in [6]. (b) Picture of Jupiter's Great Red Spot (GRS) taken on 21 April 2014 with Hubble. Credits: NASA, ESA and A. Simon (Goddard Space Flight Center). Jupiter's axis of rotation is $\Omega$.</figcaption></figure>

Another characteristic of rotating stratified fluids is that long-lived coherent vortices are often observed in geophysical conditions. This is a direct manifestation of the high-Reynolds-number turbulence of geophysical flows. Rotating turbulence can sustain vortices with a large vertical extent without stratification (due to a nearly geostrophic balance [7]). However, such vortices become strongly flattened when density stratification is strong enough, leading to a pancake-like ellipsoidal shape [8]. Notably, the elliptical shape strongly depends on the interplay between rotational effects, background shear, and the difference between the vortex's stratification and that of the ambient medium [9–11]. For instance, Mediterranean eddies are long-lived anticyclones of lenticular shape, with radii 10–100 km and thicknesses of less than 1 km (figure 1a), formed by the accumulation of warm salty water flowing from the Mediterranean Sea in the Atlantic Ocean [12]. Flattened vortices are also observed in planetary atmospheres. The most striking example is Jupiter's Great Red Spot (GRS, figure 1b), which has persisted for more than a century [13].

The origin of the long-term stability of such (nearly) isolated vortices has thus received much attention in geophysical fluid dynamics. For mathematical tractability, the quasi-geostrophic approximation is often employed to model a geophysical vortex by a fluid ellipsoid with spatially uniform potential vorticity [14]. It has also been reported that such quasi-geostrophic vortices are often quite stable over time [15,16]. However, various hydrodynamic instabilities could exist in geophysical vortices [17,18] Ultimately, these instabilities may sustain (small-scale) bulk turbulence, which would likely affect the long-term vortex stability (because of additional dissipation). In particular, the strong elliptical shape of geophysical vortices may trigger the so-called elliptical instability [19]. This is a parametric instability, due to nonlinear couplings between an elliptical flow and two waves (e.g. inertial or internal gravity waves). Yet, it remains unclear whether the elliptical instability could be triggered for realistic parameters in rotating stratified fluids [20–23]. Before assessing the relevance of this mechanism for geophysical vortices, a preliminary step is to understand the wave properties in rotating stratified fluids.

Although IGWs have been extensively studied in unbounded fluids [3], the properties of the global oscillations in bounded domains (e.g. the inertia-gravity modes, IGMs) are far from being fully understood. First insight into the mathematical problem can be obtained by considering the non-stratified regime, which has been extensively examined after Greenspan [2]. For homogeneous fluids viewed in a frame rotating at angular velocity $\Omega$ (with respect to an inertial frame), inertial modes in a bounded domain $V$ are governed by the Poincaré problem (called after Cartan [24] who revisited Poincaré's pioneering paper [25]). The latter is given by

$$ \omega^2 \nabla^2 \Phi = (2\Omega \cdot \nabla)^2 \Phi, \quad \nabla \Phi \cdot n|_{\partial V} = n \cdot (u \wedge 2\Omega)|_{\partial V}, \quad (1.1a-b) $$

where $u$ is the incompressible fluid velocity, $\Phi$ is the (reduced) hydrodynamic pressure and $n$ is the unit vector normal to the fluid boundary $\partial V$. Time dependence was assumed to be $\exp(i\omega t)$, where $|\omega| < 2|\Omega|$ is the real-valued angular frequency [2]. Equation (1.1a) is hyperbolic,
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Sketch of the mathematical model for geophysical vortices (not to scale). (a) Meridional cross-section of a vortex centred at colatitude $θ$, which is embedded within an ocean or an atmosphere (blue region). (b) Front view of a flattened ellipsoidal vortex of semi-axes $a \gg b \gg c$, which is stratified under the uniform planet's gravity $g$ (i.e. the background density $\rho_0$ decreases with increasing $z$).</figcaption></figure>

but boundary condition (1.1b) is neither of Dirichlet-type or Neumann-type (it is a mixed-typed condition, sometimes called oblique condition). Hence, problem (1.1) is an ill-posed Cauchy problem [26]. Nonetheless, it admits smooth polynomial eigenvectors when the geometry is a triaxial ellipsoid [27,28]. Moreover, these eigenvectors form a complete set in the Hilbert subspace of complex-valued divergenceless fields tangent to the boundary with the $L^2$ norm in ellipsoids [29]. Given the formal analogy between rotation and stratification [30], internal (gravity) modes also obey an ill-posed problem in non-rotating stratified fluids [31]. In the presence of stratification and rotation, the mathematical problem becomes even more complicated because IGMs obey a mixed hyperbolic-elliptic operator [32]. As for pure inertial [33,34] and internal (gravity) waves [35], IGWs can also converge after multiple reflections on solid boundaries to wave attractors [36]. Attractors are interesting geometrical structures [37], which are capable of focusing the wave energy at small length scales.

This work builds upon and extends previous studies of IGMs in confined geometries. Allen [38] and Friedlander & Siegmann [39] paid attention to IGMs in spherical (and cylindrical) containers when gravity is constant and parallel to the rotation axis. Misaligned cases were later considered in [32] with arbitrary gravity fields. Here, we aim to study IGMs using a simple model retaining the key ingredients to describe pancake-like geophysical vortices. Briefly, we consider a fluid ellipsoid subject to a constant (ambient) gravity field and uniformly rotating along an axis that is tilted with respect to gravity (to account for the full planetary rotation). This setup allows us to extend prior results about pure inertial modes in ellipsoids [27–29], while preserving polynomial eigenvectors with stratification. Finally, we characterise the wave spectrum using the mathematical theory recently presented in Colin de Verdière & Vidal [40] for inertial modes in ellipsoids. The manuscript is divided as follows. We first formulate the wave problem for general fluid volumes in §2. Then, we assume that the domain is an ellipsoid in §3 to model geophysical vortices. The wave spectrum for an ellipsoid is mathematically described and compared to numerical solutions in §4. The results are discussed in light of geophysical applications in §5, and we end the paper in §6. Additional (technical) details are provided in the Electronic Supplementary Material (ESM).

## 2. Formulation of the general model

### (a) Primitive fluid dynamics equations

As illustrated in figure 2(a), we model a vortex by a fluid domain of volume $V$ and (fictive) boundary $\partial V$, embedded within an ambient fluid (e.g. a planetary atmosphere or ocean) that is rotating at the angular velocity $\Omega = \Omega_s 1_\Omega$ with respect to an inertial frame (where $1_\Omega$ is a unit vector). We work in the frame rotating at angular velocity $\Omega$, and denote the position vector by $\boldsymbol{r} = (x, y, z)^\top$ using the Cartesian coordinates. The vortex is internally stratified in density under
the action of gravity $g$. The vortex's stratification is quantified by the Brunt-Väisälä (BV) frequency $N$. Homogeneous fluids are such that $N^2 = 0$, whereas stably stratified fluids are characterised by $N^2 > 0$. As a starting point, we employ the Boussinesq approximation to account for density variations [41]. The vortex is differentially rotating with respect to the ambient fluid outside $V$, which can be modelled by a background flow $U_0$. Laboratory experiments [10,11] and numerical calculations [9] show that, for incompressible fluids, $U_0$ is often close to a uniform-vorticity flow. The latter can be sought as $U_0 = \omega \times r + \nabla \Psi$, where $\omega$ is the vortex's rotation vector (either cyclonic or anticyclonic) and $\Psi$ is a scalar ensuring that the flow does not cross $\partial V$ (e.g. [42] in an ellipsoid). The amplitude of $U_0$ is characterised by the Rossby number $Ro = |\omega|/\Omega_s$, and we are interested in modelling stratified vortices in the regime $|Ro| \ll (N/\Omega_s)^2$. Thus, we neglect the effects of $U_0$ onto the background density and pressure. We further simplify the model by discarding baroclinic (and centrifugal effects), the very weak vortex's self-gravity, and the (weak) lateral variations of the planet's gravity at the size of a vortex. Hence, the fluid is only subject to the ambient gravity $g = -g(z)\mathbf{1}_z$ within $V$, where $\mathbf{1}_z$ is the unit vector along the vertical $z$-axis. We denote by $\rho_+ + \rho_-(z)$ the reference density within $V$ (with the mean density $\rho_*$), which is given by the hydrostatic equilibrium $\nabla P_0 = (\rho_+ + \rho_-)g$ where $P_0$ is the hydrostatic pressure.

Now, we seek small-amplitude perturbations around the reference state [$U_0, \rho_0, P_0$]. Without diffusion, the velocity $v(r, t)$ obeys the linearised Euler equation given in the rotating frame by

$$ \partial_t v + (U_0 \cdot \nabla) v + (v \cdot \nabla) U_0 + 2\Omega \times v = -\nabla\pi + (\rho/\rho_*) g \quad (2.1) $$

together with the incompressible condition $\nabla \cdot v = 0$ for an incompressible fluid, and where $\pi = p/\rho_*$ is the reduced pressure. The density $\rho(r, t)$ is then governed by

$$ \partial_t \rho/\rho_* + U_0 \cdot \nabla \rho/\rho_* = (v \cdot \mathbf{1}_z) N^2(z)/g(z), \quad (2.2) $$

where $N^2(z) = g \cdot \nabla(\rho_0/\rho_*)$ is the BV frequency within the Boussinesq approximation. The background flow has two main effects in the linear theory. First, it is responsible for a shift in frequency of the wave motions (e.g. see in [43] for a discussion in unbounded fluids). Second, it can be responsible for the onset of hydrodynamic instabilities. In particular, when $U_0$ is an elliptical flow [10,11], the so-called elliptical instability could develop [19] and be responsible for the transition towards a wave turbulence regime [44]. For mathematical simplicity, we assume $U_0 = 0$ in equations (2.1)-(2.2) to focus on the free wave motions given by

$$ \partial_t v + 2\Omega \times v = -\nabla\pi - N(z)\xi \mathbf{1}_z, \quad \nabla \cdot v = 0, \quad \partial_t \xi = (v \cdot \mathbf{1}_z) N(z), \quad (2.3a-c) $$

where we have introduced $\xi = (\rho/\rho_*) g(z)/N(z)$. The background flow will be reintroduced in future work (to investigate the outcome of the elliptical instability).

Finally, equations (2.3a-c) are supplemented with boundary conditions (BCs). We impose that the fields are regular at $r = 0$. Then, realistic conditions would be to couple the dynamics inside the vortex with the exterior flow on $\partial V$. This would amount to considering continuity conditions across $\partial V$, and decaying solutions at infinity (e.g for Rankine [45] or Lamb-Oseen [46] vortices). For mathematical simplicity, we here neglect the coupling with the exterior fluid and assume that $\partial V$ is a rigid boundary. Consequently, the velocity obeys the no-penetration BC $v \cdot n|_{\partial V} = 0$ on $\partial V$, where $n$ is the unit (outward) vector normal to the boundary. Without diffusion, no other physical BCs have to be explicitly enforced in the model. Indeed, the density and pressure BCs automatically follow from the velocity BC. In practice, the pressure satisfies a mixed BC that is given by the continuity of the normal component of equation (2.3a) on $\partial V$ (see the ESM).

## (b) Mixed hyperbolic-elliptic problem

We seek modal solutions of equations (2.3a-c) given by

$$ [v, \pi, \xi](r, t) = [u, \Phi, \zeta](r) \exp(i\omega t), \quad (2.4) $$

where $[u, \Phi, \zeta]$ are complex-valued amplitudes depending on space, and with $\omega \in \mathbb{C}$. The wave properties are deeply tied to the nature of the mathematical problem. Instead of solving the
primitive equations, it is worth considering a master equation that governs the evolution of either **u** or **Φ**. If we multiply Euler equation (2.3a) by i**ω** and replace i**ωζ** by buoyancy equation (2.3c), we obtain for **ω** ≠ 0

$$ -\omega^2 \mathbf{u} + 2\mathbf{\Omega} \times (i\omega\mathbf{u}) + N^2(z)(\mathbf{u} \cdot \mathbf{1}_z)\mathbf{1}_z = -i\omega\nabla\Phi, \quad \nabla \cdot \mathbf{u} = 0, \qquad (2.5a,b) $$

Equations (2.5a,b) and the no-penetration BC $\mathbf{u} \cdot \mathbf{n}|_{\partial V} = 0$ show that the velocity is an appropriate variable to explore the wave properties. Note that a similar equation can be obtained for the fluid particle displacement vector [47]. On the contrary, non-oscillatory modes with $\omega = 0$ require a specific consideration. For unstratified fluids with $N = 0$, they are called geostrophic modes because they obey the geostrophic balance $2\mathbf{\Omega} \times \mathbf{u} = -\nabla\Phi$. When $N(z) > 0$, the steady modes are two-dimensional (2D) such that $\mathbf{u} \cdot \mathbf{1}_z = 0$, and the density perturbation obeys the diagnostic equation $\rho_*(2\mathbf{\Omega} \cdot \nabla\mathbf{u}) = \nabla\zeta \times \mathbf{1}_z$ (resulting from the thermal wind balance).

Because of historical reasons (dating back to Poincaré [25]), the problem is generally recast into a scalar equation for the pressure $\Phi$ [32,39]. As shown in the ESM, the pressure is governed in V by a second-order equation whose higher-order terms are given by

$$ [\omega^2 - N^2(z)]\nabla^2\Phi + N^2(z)(\mathbf{1}_z \cdot \nabla)^2\Phi - (2\mathbf{\Omega} \cdot \nabla)^2\Phi + \dots = 0 \qquad (2.6) $$

when $|\omega| \neq \omega_{\pm}(z)$ with

$$ 2\omega_{\pm}^2(z) = [N^2(z) + 4|\mathbf{\Omega}|^2] \pm \sqrt{[N^2(z) + 4|\mathbf{\Omega}|^2]^2 - 16N^2(z)(\mathbf{\Omega} \cdot \mathbf{1}_z)^2}. \qquad (2.7) $$

The lower-order terms in equation (2.6), which vanish when $N$ is constant, are here omitted for concision but are given in the ESM. When $|\omega| = \omega_{\pm}(z)$, the pressure instead obeys a first-order equation (see the ESM). Finally, the pressure obeys a mixed BC on $\partial V$ (i.e. neither a pure Dirichlet nor Neumann BC), which is difficult to take into account for numerical computations. However, the nature of the mathematical problem can be simply determined from the pressure equation by using microlocal analysis [48]. This is a branch of mathematics, which studies the solutions of PDE using (notably) asymptotic geometric techniques. Microlocal tools have already proven useful in physics to study wave problems [49,50]. In particular, the nature of equation (2.6) is governed by a scalar quantity that is called the principal symbol. The latter is obtained from the equation by freezing the variable coefficients and keeping only the highest-degree terms. In practice, it amounts to transforming the derivatives as $\nabla \leftrightarrow i\mathbf{k}$, where $\mathbf{k} = (k_1, k_2, k_3)^{\top}$ is a real-valued wave vector at the position $(x, y, z)$. The principal symbol is then a homogeneous polynomial in $\mathbf{k}$, which encapsulates many properties of the spectral problem [37,51]. The problem is elliptic when the principal symbol is invertible, and hyperbolic when the principal symbol vanishes. The vanishing condition on the principal symbol gives

$$ |\mathbf{k}|^2 \omega^2 - [N^2(z) |\mathbf{1}_z \times \mathbf{k}|^2 + (2\mathbf{\Omega} \cdot \mathbf{k})^2] = 0, \qquad (2.8) $$

which is the dispersion relation of IGWs in unbounded fluids when $N$ is constant [3]. The wave-like equation is thus hyperbolic when equation (2.8) is satisfied for some non-zero wave vectors $\mathbf{k} \in \mathbb{R}^3$. Otherwise, it is elliptic (except when $\omega = 0$). Next, the bounds for the hyperbolic domain can be obtained from equation (2.8). The latter defines a quadric cone given by $(\mathbf{g}^*)^{ij} k_i k_j = 0$ (using Einstein notation) when $\omega \neq 0$, where $\mathbf{g}^*$ is the associated metric tensor. The nature of this quadric cone depends on the eigenvalues of $\mathbf{g}^*$, which are given by $\mu_0 = \omega^2 - N^2(z)$ and $\mu_{\pm} = \omega^2 - \omega_{\pm}^2(z)$. When the eigenvalues are non-zero and of the same sign, the surface is elliptic. On the contrary, the surface is hyperbolic when the eigenvalues are non-zero but of different signs. The above microlocal description agrees with the analysis presented in Section 3 of Friedlander & Siegmann [32]. The wave-like equation is thus hyperbolic in V when $\omega_{-}(z) < |\omega| < \omega_{+}(z)$, whereas it is elliptic in V when $0 < |\omega| < \omega_{-}(z)$ or $|\omega| > \omega_{+}(z)$.
(c) Spectral problem for the velocity

Further wave properties can be obtained by considering the velocity equation. We introduce the Hilbert space $\mathcal{V}$ spanned by complex-valued vector fields that are square-integrable in $V$. The complex-valued inner product in this Hilbert space is given by

$$ \langle u, v \rangle := \int_V u^\dagger \cdot v \, dV, \qquad (2.9) $$

where $\dagger$ is the complex conjugate. The associated norm is $||u|| = \langle u, u \rangle^{1/2}$. Then, we denote by $\mathcal{V}^0 \subset \mathcal{V}$ the closed subspace that is orthogonal, with respect to inner product (2.9), to the space spanned by vector fields made of gradients of smooth functions. A smooth element in $\mathcal{V}^0$ is divergenceless and tangent to the boundary [29], that is

$$ \mathcal{V}^0 := \{u \in \mathcal{V} \mid \nabla \cdot u = 0 \text{ in } V, u \cdot n = 0 \text{ on } \partial V\}. \qquad (2.10) $$

Finally, we introduce the orthogonal projector $\mathbb{L}$ from $\mathcal{V}$ onto $\mathcal{V}^0$ (called Leray projector after [52]). For any vector $e \in \mathcal{V}$ written as $e = \nabla \Psi + \nabla \times A$ (using a Helmholtz decomposition) for some potentials $[\Psi, A]$, the projected vector is defined by $\mathbb{L}(e) := e - \nabla \Psi$. This projector is used in the mathematical analysis of incompressible flows [53], but also in some numerical studies [54].

Equipped with the above definitions, we can seek solutions of equation (2.5) with $u \in \mathcal{V}^0$ in ellipsoids. We apply the orthogonal projector to the equation to remove the pressure term, and find that $(\omega, u)$ obey a quadratic eigenvalue problem (QEP). The latter is given by

$$ \mathcal{Q}_{\omega}(u) = 0, \quad \mathcal{Q}_{\omega} := -\omega^2 \mathcal{I} + \omega i \mathcal{C} + \mathcal{K}, \qquad (2.11a,b) $$

where $\mathcal{I}$ is the identity operator, and with the two operators

$$ i\mathcal{C}(u) := i\mathbb{L}[2\Omega \times u], \quad \mathcal{K}(u) := \mathbb{L}[N^2(z)(u \cdot 1_z)1_z]. \qquad (2.12a,b) $$

The operator $i\mathcal{C}$, called the Poincaré operator [40], is a bounded self-adjoint operator in $\mathcal{V}^0$ [29]. The buoyancy operator $\mathcal{K}$ has properties given by Proposition 2.1.

**Proposition 2.1.** *The operator $\mathcal{K}$ is bounded, self-adjoint and positive in $\mathcal{V}^0$.*

*Proof.* We must show that $\mathcal{K}$ is symmetric, bounded, and positive in $\mathcal{V}^0$. We verify that $\mathcal{K}$ is symmetric since $\langle a, \mathcal{K}(b) \rangle = \langle \mathcal{K}(a), b \rangle$ for any $(a, b) \in \mathcal{V}^0 \times \mathcal{V}^0$. Moreover, we have

$$ ||\mathcal{K}(u)||^2 = \langle \mathcal{K}(u), N^2(u \cdot 1_z)1_z \rangle - \langle \mathcal{K}(u), \nabla\Phi \rangle = \langle \mathcal{K}(u), N^2(u \cdot 1_z)1_z \rangle $$

since $\langle \mathcal{K}(u), \nabla\Phi \rangle = 0$ using the definition of $\mathbb{L}$. Using Cauchy-Schwarz inequality, we then obtain $||\mathcal{K}(u)||^2 \le ||\mathcal{K}(u)|| ||N^2(u \cdot 1_z)1_z|| \le \max(N^2)||\mathcal{K}(u)|| ||u||$, showing that $\mathcal{K}$ is bounded. Finally, $\mathcal{K}$ is positive because $\langle u, \mathcal{K}(u) \rangle \ge 0$ if $N^2 \ge 0$. ■

We want to characterise the spectrum $\sigma$ of the quadratic spectral form $\mathcal{Q}_{\omega}$ in ellipsoids, which is the set of complex numbers $\omega$ for which $\mathcal{Q}_{\omega}$ is not continuously invertible. Actually, the spectral problem $\mathcal{Q}_{\omega}(u) = 0$ can be converted into a standard eigenvalue problem (SEP) for a bounded self-adjoint operator in $\mathcal{V}^0 \times \mathcal{V}^0$.

*Proof.* This is a general result for spectral families of the form $\mathcal{Q}_{\omega} := -\omega^2 + \omega A + B$ acting on a Hilbert space $\mathcal{H}$, where $[A, B]$ are bounded self-adjoint operators and with $B \ge 0$ [55,56]. When $\omega \ne 0$, the QEP $\mathcal{Q}_{\omega}(u) = 0$ can be recast as a SEP for the operators $\mathcal{L}$ or $\mathcal{T}$ given by

$$ \mathcal{L} := \begin{pmatrix} 0 & \mathcal{I} \\ \mathcal{B} & \mathcal{A} \end{pmatrix}, \quad \mathcal{T} := \begin{pmatrix} \mathcal{A} & \mathcal{B}^{1/2} \\ \mathcal{B}^{1/2} & 0 \end{pmatrix}, \qquad (2.13a,b) $$

which are both acting on the Hilbert space $\mathcal{H} \times \mathcal{H}$ equipped with the inner product $\langle\langle z, \zeta \rangle\rangle = \langle z_1, \zeta_1 \rangle + \langle z_2, \zeta_2 \rangle$ with $z = (z_1, z_2)^\top$ and $\zeta = (\zeta_1, \zeta_2)^\top$. Here, $\mathcal{T}$ is a bounded self-adjoint operator (where $B^{1/2}$ is the square root of $B$, which is also self-adjoint [57]), whereas $\mathcal{L}$ has no particular symmetries. The spectra of $\mathcal{Q}_{\omega}$, $\mathcal{L}$ and $\mathcal{T}$ are identical outside $\omega = 0$. ■
The spectrum $\sigma(Q_\omega)$ is thus the disjoint union of the point spectrum given by

$$ \sigma_P(Q_\omega) := \{\omega \in \mathbb{C} \mid \exists \mathbf{u} \neq \mathbf{0} \ Q_\omega(\mathbf{u}) = \mathbf{0}\}, \qquad (2.14) $$

which is the subset spanned by the eigenvalues of QEP (2.11), and the remaining continuous spectrum (i.e. the set of $\omega$ for which $T - \omega I$ is injective and has a dense range, but is not surjective). This mathematical distinction is of direct physical interest to IGMs. Indeed, the point spectrum is associated with (regular) square-integrable eigenvectors in diffusionless fluids. On the contrary, a non-empty continuous spectrum is often characterised by (almost) periodic attracting trajectories obtained from ray tracing techniques [37]. These geometrical structures, called attractors [37], are associated with nearly singular (i.e. non-square-integrable) velocity fields when diffusion is vanishingly small [58]. Given the symmetries of the different operators, some general properties of the spectrum are given in Proposition 2.2.

**Proposition 2.2.** The spectrum of QEP (2.11) is real-valued and even. Moreover, the spectrum is bounded by $|\omega| \le \omega_{\max}$ with $\omega_{\max} \le \Omega_s + \sqrt{\Omega_s^2 + \max(N^2)}$. Finally, two eigen-pairs $(\omega_1, \mathbf{u}_1)$ and $(\omega_2, \mathbf{u}_2)$ with $\omega_1 \neq \omega_2$ are not orthogonal with respect to inner product (2.9), but satisfy the two integral properties

$$ \omega_1\omega_2\langle\mathbf{u}_2, \mathbf{u}_1\rangle = -\langle\mathbf{u}_2, \mathcal{K}(\mathbf{u}_1)\rangle, \quad (\omega_2 + \omega_1)\langle\mathbf{u}_2, \mathbf{u}_1\rangle = \langle\mathbf{u}_2, i\mathcal{C}(\mathbf{u}_1)\rangle. \qquad (2.15a,b) $$

*Proof.* $\omega \in \mathbb{R}$ result from the fact that the operator $T$ defined in equation (2.13b) is self-adjoint when $\mathcal{K} \ge 0$. Then, the complex-conjugate of QEP (2.11) gives

$$ -\omega^2 \mathbf{u}^{\dagger} - \omega i \mathcal{C}(\mathbf{u}^{\dagger}) + \mathcal{K}(\mathbf{u}^{\dagger}) = \mathbf{0}, $$

which shows that $(-\omega, \mathbf{u}^\dagger)$ is an eigen-pair of the QEP. The upper bound of $\omega_{\max}$ can be obtained using a min-max principle [55]. Finally, the orthogonality conditions can be found using vector manipulations of the QEP [59]. ■

## 3. Ellipsoidal model for a constant BV frequency

We now simplify the general model to consider an idealised configuration of geophysical vortices. To get physical insight into the problem, it is useful to consider a shape that is amenable to mathematical analysis. For simplicity, we assume that the boundary $\partial V$ is a triaxial ellipsoid of semi-axes $[a, b, c]$, where the c-axis is aligned with gravity (figure 2b). The use of ellipsoids has a long-standing history in the study of geophysical vortices [14], and is supported by recent fluid dynamic experiments [10,11]. In the rotating frame, the boundary is expressed using the Cartesian coordinates as $F(x, y, z) = 1$ with $F(x, y, z) := (x/a)^2 + (y/b)^2 + (z/c)^2$, such that the unit normal is $n := \nabla F/||\nabla F||$. Moreover, geophysical vortices with $N \gg 2\Omega_s$ are usually strongly flattened with $a, b \gg c$ (e.g. Jupiter's GRS is expected to have a depth of only a few hundreds of km [11]). In such vortices, the ambient gravity is expected to weakly varies with depth. Thus, we also assume that $g = -g_1z$ is uniform and consider a background density $\rho_* + \rho_0(z)$ that varies linearly with $z$ as

$$ \rho_0(z)/\rho_* = -(N^2/g) z, \qquad (3.1) $$

where $N$ is a constant BV frequency. With these assumptions, the eigensolutions can be characterised using finite-dimensional arguments in ellipsoids.

### (a) Polynomial description

Since an ellipsoid is a quadratic surface and expression (3.1) is a polynomial, we can seek eigenvectors using polynomial functions in the Cartesian coordinates. We denote by $P_n: \mathbb{R}^3 \to \mathbb{C}$ the space spanned by the scalar Cartesian monomials $x^i y^j z^k$ of degree $i + j + k$ smaller or equal to $n$ with $n \ge 1$. Its dimension is $N(P_n) = (n+1)(n+2)(n+3)/6$. Then, we denote by $V_n$ the space of vectors in ellipsoids with coefficients in $P_n$ and define the subspace $V_n^0 := V_n \cap V^0$
of dimension $N(\mathcal{V}_n^0) = n(n+1)(2n+7)/6$, whose elements are divergenceless vector fields that are tangent to the boundary $\partial V$. Finally, we introduce the subspaces $\mathcal{W}_n^0 := (\mathcal{V}_{n-1}^0)^\perp \cap \mathcal{V}_n^0$ of dimension $N(\mathcal{W}_n^0) = n(n+2)$, which are spanned by vectors in $\mathcal{V}_n^0$ of degree $n$ that are orthogonal to polynomial vectors of degree $\le n-1$. These polynomial spaces are key to characterise the spectrum because of Proposition 3.1.

**Proposition 3.1.** The spaces $\mathcal{V}_n^0$ with $n \ge 1$ are invariant by the action of $\mathcal{C}$ and $\mathcal{K}$ in triaxial ellipsoids, that is $\mathcal{C}|_{\mathcal{V}_n^0} \subseteq \mathcal{V}_n^0$ and $\mathcal{K}|_{\mathcal{V}_n^0} \subseteq \mathcal{V}_n^0$.

*Proof.* The invariance of $\mathcal{V}_n^0$ under the action of $\mathcal{C}$ has been demonstrated elsewhere [29,40]. Here, we can use similar arguments to show that it is invariant by the action of $\mathcal{K}$. Briefly, we must prove that the Leray projector $\mathbb{L}$ is well-defined on $\mathcal{V}_n^0$. This is ensured by Theorem 2.4 in [60], which shows that there is a unique scalar $\Psi \in \mathcal{P}_{n+1}$ in the definition of $\mathbb{L}$.

We deduce from Proposition 3.1 that $\mathcal{V}_n^0$ is invariant by the action of $\mathcal{Q}_\omega$. Moreover, since $\oplus_{n \in \mathbb{N}^*} \mathcal{V}_n^0$ is dense in $\mathcal{V}^0$ [40], it shows that $\sigma(\mathcal{Q}_\omega)$ is pure point in ellipsoids (i.e. $\sigma = \sigma_P$). The continuous spectrum is thus empty, and equation (2.5) admits square-integrable eigenvectors in $\mathcal{V}_n^0$. From a physical viewpoint, an empty continuous spectrum means that there are no singular velocity fields associated with wave attractors. This is a difference with other geometries, in which wave attractors can exist even when $N$ is constant (e.g. in a trapezoid [36]). However, this does not preclude the existence of (almost) periodic wave trajectories using ray tracing techniques. The classical example is the pure inertial wave problem within a sphere. Periodic trajectories have been found in this geometry using ray theory [61] but, since the continuous spectrum is also empty in a full sphere [29,62], these structures are not associated with singular global modes for vanishingly small diffusion. Here, because the spectrum is pure point, we can solve the eigenvalue problem by restricting $u$ to $\mathcal{V}_n^0$ in ellipsoids.

## (b) Galerkin algorithm

An accurate numerical description can be developed by seeking $u \in \mathcal{V}_n^0$ as $u = \sum_{j \ge 1} \alpha_j e_j$, where $\{\alpha_j\}_{j \ge 1}$ are complex-valued coefficients, and $\{e_j\}_{j \ge 1}$ are real-valued basis vector elements of $\mathcal{V}_n^0$ [63] normalised such that $||e_j|| = 1$. Then, we substitute the above expansion into the QEP $\mathcal{Q}_\omega(u) = 0$, and project the entire equation on every $e_i$ using inner-product (2.9) to cancel out the spatial dependence (Galerkin method [63]). This gives the finite-dimensional QEP

$$(-\omega^2 M_n + \omega i C_n + K_n) \alpha = 0 \qquad (3.2)$$

where $\alpha = (\alpha_1, \alpha_2, \dots)^\top$ is the eigenvector, and where $[M_n, C_n, K_n]$ are three real-valued matrices of size $N(\mathcal{V}_n^0) \times N(\mathcal{V}_n^0)$ and of elements $M_{n,ij} = \langle e_i, e_j \rangle$, $C_{n,ij} = \langle e_i, 2\Omega \times e_j \rangle$ and $K_{n,ij} = \langle e_i, N^2(e_j \cdot 1_z)1_z \rangle$. The matrix $M_n$ is positive definite but non-diagonal (because the chosen basis elements are not mutually orthogonal by construction in triaxial ellipsoids [63]). The matrix $C_n$ is antisymmetric, and $K_n$ is a real-valued positive matrix. Further properties of $[C_n, K_n]$ are given below in Proposition 3.2.

**Proposition 3.2.** $\dim(\ker C_n) = \lceil n/2 \rceil$ and $\dim(\ker K_n) = n(n+1)(n+2)/6$.

*Proof.* When $N=0$, $\ker(C_n)$ is associated with geostrophic modes of degree $n$. Its dimension was proved in [40]. When $N \ne 0$, $\ker(K_n)$ is associated with 2D motions such that $u \cdot 1_z = 0$. The latter can be sought as $u = \nabla \times (\Psi 1_z)$ where $\Psi = p[1-F(x, y, z)]$ is a stream function and with $p \in \mathcal{P}_{n-1}$. The dimension of $\ker(K_n)$ is thus $\dim(\mathcal{P}_{n-1}) = n(n+1)(n+2)/6$.

For numerical convenience, we convert the finite-dimensional QEP into an equivalent finite-dimensional SEP when $\omega \ne 0$. The finite-dimensional counterpart of the self-adjoint operator
<figure><img src="image_6.png" /><figcaption>Figure 3. Sparse upper triangular block matrices $A_n$ (left) and $B_n$ (right) for $n = 5$ and $1_\Omega = 1_z$.</figcaption></figure>

defined in equation (2.13b) is not well suited for numerical computations (because of the square root). Instead, we consider the SEP $L_n z = \omega z$ with $z = (\alpha, \beta)^\top$ and the matrix $L_n$ given by

$$L_n := \begin{pmatrix} 0 & I \\ B_n & iA_n \end{pmatrix}, \quad A_n := M_n^{-1}C_n, \quad B_n := M_n^{-1}K_n, \qquad (3.3a-c)$$

where $I$ is the identity matrix. As illustrated in figure 3, the matrices $A_n$ and $B_n$ are sparse and have an upper triangular block structure (which directly results from Proposition 3.1). It has a practical consequence in the unstratified case $N = 0$. The eigenvalues of the pure inertial modes, given by $A_n \alpha = \omega \alpha$, can indeed be computed separately for every space $W_n^0$ by seeking the eigenvalues of every diagonal block of $A_n$. Yet, since the matrix $L_n$ has not an upper triangular block structure, we cannot compute separately the eigenvalues for every degree $n$ when $N \neq 0$.

We have modified the bespoke numerical code SHINE, initially developed for compressible fluids [64,65], to implement the above algorithm for Boussinesq fluids. The matrices $[M_n, C_n, K_n]$ are first built analytically. Since the basis vector elements are not mutually orthogonal in an ellipsoid [63], the matrices become numerically ill-conditioned when the polynomial degree is large. Hence, we perform the computations using extended-precision algorithms. We can illustrate the Galerkin algorithm by considering the six polynomial solutions $u \in \mathcal{V}_1^0$. The latter are sought in the form of uniform-vorticity flows as [42]

$$u = \alpha \times r + \nabla \psi, \quad \psi = \alpha_x \frac{c^2 - b^2}{b^2 + c^2} yz + \alpha_y \frac{a^2 - c^2}{a^2 + c^2} xz + \alpha_z \frac{b^2 - a^2}{a^2 + b^2} xy, \qquad (3.4a-b)$$

and with $\alpha = (\alpha_x, \alpha_y, \alpha_z)^\top$. Next, the block matrices of $L_1$ are given by

$$B_1 = \begin{pmatrix} \frac{N^2 c^2}{b^2 + c^2} & 0 & 0 \\ 0 & \frac{N^2 c^2}{a^2 + c^2} & 0 \\ 0 & 0 & 0 \end{pmatrix}, \quad A_1 = \begin{pmatrix} 0 & -\frac{2\Omega_z a^2}{a^2 + c^2} & \frac{2\Omega_y a^2}{a^2 + b^2} \\ \frac{2\Omega_z b^2}{b^2 + c^2} & 0 & -\frac{2\Omega_x b^2}{a^2 + b^2} \\ -\frac{2\Omega_y c^2}{b^2 + c^2} & \frac{2\Omega_x c^2}{a^2 + c^2} & 0 \end{pmatrix}, \qquad (3.5a,b)$$

and with $\Omega = (\Omega_x, \Omega_y, \Omega_z)^\top$. Finally, the eigenvalue-eigenvector solutions of $L_1$ can be solved analytically to extend the linear solutions obtained by Vantieghem [27] when $N = 0$ and $1_\Omega = 1_z$. However, we must rely on numerical computations to obtain $L_n$ for higher degrees $n$.
(c) Enumeration of the eigenvalues

We can explicitly count the number of non-zero and zero eigenvalues of the QEP when $u \in \mathcal{V}_n^0$. Since the spectra of the primitive equations, the one of $\mathcal{Q}_\omega$ and the one of $\mathcal{L}$ are identical when $\omega \neq 0$, their finite-dimensional counterparts have the same number of non-zero eigenvalues.

**Proposition 3.3.** The number of non-zero eigenvalues of QEP (3.2) is $n(n+1)(n+5)/3$.

*Proof.* The number of non-zero eigenvalues can be estimated from the primitive equations

$$i\omega u + 2\Omega \times u = -\nabla\pi - \zeta N \mathbf{1}_z, \quad \nabla \cdot u = 0, \quad i\omega\zeta = (u \cdot \mathbf{1}_z) N.$$

Solutions of the latter equations belong to the space spanned by $u \in \mathcal{V}_n^0$ and $\zeta \in \mathcal{P}_n$, whose dimension is $N(\mathcal{V}_n^0) + N(\mathcal{P}_n) = (n+1)(n^2+4n+2)/2$. Next, we count the number of zero eigenvalues. From the density equation, we deduce that $\omega = 0$ is associated with 2D flows such that $u \cdot \mathbf{1}_z = 0$. Such 2D flows $u \in \mathcal{V}_n^0$ can be sought as $u = \nabla \times (\Psi \mathbf{1}_z)$ where $\Psi = p[1 - F(x, y, z)]$ is a stream function and with $p \in \mathcal{P}_{n-1}$. Assuming that $\Omega = (\Omega_x, \Omega_y, \Omega_z)^\top$, we obtain $2\Omega_z \partial_x \Psi = -\partial_x \pi$ and $2\Omega_z \partial_y \Psi = -\partial_y \pi$ from the horizontal components of the momentum equation. Then, we deduce that $\pi = -2\Omega_z \Psi + h(z)$, where $h(z)$ is an arbitrary polynomial function of $z$ and of degree $n+1$ (with $n \ge 1$). Finally, the density perturbation is obtained from the vertical component of the momentum equation, which gives $\zeta N = 2(\Omega_y \partial_y \Psi + \Omega_x \partial_x \Psi) - \partial_z \pi$. Thus, $\omega = 0$ is associated with an eigenspace of dimension $\dim(\mathcal{P}_{n-1}) + (n+1) = (n+1)(n^2+2n+6)/6$. Therefore, the number of non-zero eigenvalues is given by $n(n+1)(n+5)/3$. ■

Equipped with Proposition 3.3, we can enumerate the zero eigenvalues of the matrix $L_n$, which has $2N(\mathcal{V}_n^0)$ eigenvalues. Hence, the algebraic multiplicity of the zero eigenvalue $\omega = 0$ is given by $n(n+1)(n+2)/3$. However, $\omega = 0$ is not only associated with standard eigenvectors. The geometric multiplicity $m_L$ of $\omega = 0$ is the dimension of $\ker L_n$. The null space of the equation $L_n(\alpha, \beta)^\top = 0$ is spanned by the vector elements such that $\beta = 0$ and $K_n\alpha = 0$. Therefore, we have $\ker L_n = \ker K_n \oplus \{0\}$ and $\omega = 0$ has the geometric multiplicity $m_L = n(n+1)(n+2)/6$ according to Proposition 3.2. However, since $L_n$ is not Hermitian, the eigenvalue $\omega = 0$ has an algebraic multiplicity $r_L > m_L$. The missing $r_L - m_L = n(n+1)(n+2)/6$ polynomial solutions are associated with generalized eigenvectors, which belong to $\ker(L_n^2)$.

**Proposition 3.4.** The eigenvalue $\omega = 0$ of the matrix $L_n$ has the geometric multiplicity $m_L = n(n+1)(n+2)/6$, and the algebraic multiplicity $r_L = 2m_L$ (because of generalized eigenvectors).

The vector space associated with $\omega = 0$ is thus partly spanned by generalized eigenvectors, which here satisfy $L_n z \neq 0$ and $L_n^2 z = 0$. Back in the temporal domain, they correspond to solutions that are not steady but linear in $t$. However, these admissible solutions are mathematical artefacts due to the formulation of the quadratic problem in terms of the matrix $L_n$. Indeed, QEP (3.2) can be converted into an equivalent self-adjoint SEP [59], which is the finite-dimensional counterpart of the self-adjoint operator defined in equation (2.13b). This shows that, from a physical viewpoint, $\omega = 0$ is only associated with standard eigenvectors. Therefore, we will not further elaborate on the (artificial) generalized eigenvectors of $L_n$ below. Following our prior investigation of the non-stratified problem [40], we can combine polynomial computations for fixed values of $n$ and spectral theory for $n \to \infty$ to study the eigensolutions in ellipsoids.

4. Point spectrum in ellipsoids

We have shown that $\sigma(\mathcal{Q}_\omega)$ is pure point in ellipsoids when $N$ is constant. For this configuration, we can further characterise the eigensolutions. In particular, the spectrum can be divided into two subsets. If an eigenvalue is isolated in the spectrum and of finite multiplicity (when $n \to \infty$), it belongs to the discrete spectrum $\sigma_D$. Otherwise, it belongs to the essential spectrum $\sigma_E := \sigma \setminus \sigma_D$.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Sketch of the bounded pure point spectrum in triaxial ellipsoids (not to scale) when $N$ is constant. The spectrum is symmetric with respect to 0 and bounded by $|\omega| \le \omega_+$ (Proposition 4.1). The velocity equation is either hyperbolic or elliptic in volume. The spectrum is spanned by IGMs and low-frequency (surface) modes.</figcaption></figure>

(i.e. the complement to the discrete spectrum). Since $\oplus_{n \in \mathbb{N}^*} \mathcal{V}_n^0$ is dense in $\mathcal{V}^0$, we conclude from Proposition 3.4 that the eigenvalue $\omega = 0$ has an infinite multiplicity when $n \to \infty$ (i.e. $\omega = 0$ belongs to $\sigma_E$). Now, it remains to characterise the other eigensolutions.

When $N$ is constant, the frequencies $\omega_\pm$ defined in equation (2.7) are constant in the volume. Moreover, it is possible to find a better estimate of $\omega_{\max}$ (compared to the one given in Proposition 2.2). In the aligned case, it was shown in [39] that $\omega_{\max} = \omega_+$ with $\omega_+ = \max(N, 2\Omega_s)$. Thus, $\sigma_3$ is empty in ellipsoids when $\Omega$ and $1_z$ are collinear. We also show below in Proposition 4.1 that $\omega_{\max} = \omega_+$, even when $\Omega$ and $1_z$ are not collinear.

**Proposition 4.1.** When $N$ is constant, the pure point spectrum is bounded by $\omega_+$ in an ellipsoid.

*Proof.* We introduce the projector operator $\mathbb{P}$ acting on a pair $z = (u, \zeta)^\top$ by $\mathbb{P}(z) := (\mathbb{L}u, \zeta)^\top$, where $\mathbb{L}$ is the Leray projector introduced in §2. Then, we can reformulate the spectral problem as the spectral properties of the linear operator $\mathcal{P}$ defined by $\mathcal{P}(z) := \mathbb{P}(\mathcal{E}\mathbb{P}(z))$ with

$$ \mathcal{E} = i \begin{pmatrix} 0 & -2\Omega_3 & 2\Omega_2 & 0 \\ 2\Omega_3 & 0 & -2\Omega_1 & 0 \\ -2\Omega_2 & 2\Omega_1 & 0 & N \\ 0 & 0 & -N & 0 \end{pmatrix}. $$

$\mathcal{P}$ is a self-adjoint operator for square-integrable fields. We have $||\mathbb{P}|| = 1$, and the eigenvalues of the Hermitian matrix $\mathcal{E}$ are $\pm\omega_\pm$ such that $||\mathcal{E}|| \le \omega_+$. Thus, the spectrum is bounded by $\omega_+$. ■

Therefore, we can divide the allowable frequency range into disjoint intervals with $I_1 = [-\omega_+, -\omega_-] \cup [-\omega_-, \omega_+]$ and $I_2 = -\omega_-, \omega_- \setminus \{0\}$. The velocity equation is hyperbolic in $V$ when $\omega \in I_1$, and elliptic in $V$ when $\omega \in I_2$. We also define the associated spectra $\sigma_i = \sigma_P \cap I_i$ for $i = \{1, 2\}$ and, when $u$ is restricted to $\mathcal{V}_n^0$, the subsets are denoted by $[\sigma_{1,n}, \sigma_{2,n}]$ such that $\sigma_{i,n-1} \subseteq \sigma_{i,n}$. In unbounded fluids, a continuum of IGWs exists with angular frequencies $\omega \in I_1$ given by dispersion relation (2.8). We will show below that the spectrum $\sigma_1 \subset \sigma_E$ is dense in $I_1$ within ellipsoids. On the contrary, $\sigma_2$ is known to be empty in unbounded fluids (which agrees with prior experimental observations [66]). This is consistent with the fact that the velocity satisfies an elliptic-type equation when $\omega \in I_2$ (i.e. plane-wave solutions are evanescent in these intervals). Yet, a countable number of eigenfrequencies can exist in some bounded geometries when $|\omega| < \omega_-$ [38,39]. The eigensolutions with $\omega \in \sigma_1$ are referred as class-I solutions, and those with $\omega \in \sigma_2$ are called class-II solutions in [38,39]. We describe the two spectra below.

## (a) Low-frequency (surface) modes

### (i) Illustrative examples

Contrary to the free-space case, it is known that $\sigma_2$ is not empty when $\Omega \propto 1_z$ in some bounded geometries (e.g. in cylinders or spheres [21,39]). In triaxial ellipsoids, we can also show that $\sigma_2 \neq \emptyset$.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Spectrum $\sigma_2$ in the aligned case $1_\Omega = 1_z$. (a) Eigenvalues $\omega_{3,5} \ge 0$ in $\sigma_{1,1}$ for a sphere with $b/a = c/a = 1$ (solid purple curves) and an ellipsoid with $b/a = 0.7$ and $c/a = 0.1$ (dashed olive curves). Red region shows $\sigma_2$ in which the problem is hyperbolic. (b) Number of eigenvalues in $\sigma_{1,n}$ for a sphere ($b/a = c/a = 1$) and a flattened ellipsoid ($b/a = 0.7$ and $c/a = 0.1$). Colour bar shows the value of $N/\Omega_s$.</figcaption></figure>

when $\Omega \cdot 1_z \ne 0$ (otherwise, $\omega_- = 0$ and $\sigma_1$ fills the entire interval $0 < |\omega| < \omega_+$). To this end, we first consider the eigenvalues for $u \in \mathcal{V}_1^0$, which are the roots of the characteristic equation $\det(L_1 - \omega I) = 0$. The six roots $\omega_j$ of the characteristic polynomial are the two degenerate solutions $\omega_{1,2} = 0$ and the four solutions

$$ \omega_{3,4} = \pm \sqrt{\frac{\beta_4}{2\Delta} + \frac{\sqrt{(\beta_4/\Delta)^2 - 4\beta_2/\Delta}}{2}}, \quad \omega_{5,6} = \pm \sqrt{\frac{\beta_4}{2\Delta} - \frac{\sqrt{(\beta_4/\Delta)^2 - 4\beta_2/\Delta}}{2}}, \quad (4.1a-b) $$

where we have introduced the positive coefficients

$$ \begin{aligned} \beta_4 &= N^2 c^2 (a^2 + b^2)(a^2 + b^2 + 2c^2) + 4\Omega_x^2 b^2 c^2 (b^2 + c^2) + 4\Omega_y^2 a^2 c^2 (a^2 + c^2) + 4\Omega_z^2 a^2 b^2 (a^2 + b^2), \\ \beta_2 &= c^4 N^2 [4(\Omega_x^2 b^2 + \Omega_y^2 a^2) + N^2 (a^2 + b^2)], \\ \Delta &= a^4 b^2 + a^4 c^2 + a^2 b^4 + 2a^2 b^2 c^2 + a^2 c^4 + b^4 c^2 + b^2 c^4. \end{aligned} $$

The only non-zero solutions in the non-stratified case $N = 0$ are $\omega_{3,4} = \pm \sqrt{\beta_4/\Delta}$, which reduce to the eigenvalues of the pure inertial modes in $\mathcal{V}_1^0$ obtained by Vantieghem [27] when $\Omega_x = \Omega_y = 0$. Solutions (4.1) are illustrated in figure 5(a) for a sphere and a flattened (triaxial) ellipsoid. Several points are worth commenting on in the figure. First, $\omega_3$ and $\omega_5$ smoothly vary as a function of $N$. We observe that $\omega_3 \in \sigma_{1,1}$ for any value of $N$ but, on the contrary, $\omega_5$ can either belong to $\sigma_{1,1}$ or $\sigma_{2,1}$ when $N$ is varied. For instance, the transition occurs here at $N/\Omega_s > 4$ in a sphere and at a much larger value $N/\Omega_s \gtrsim 25$ in the pancake-like ellipsoid.

Next, we show in figure 5(b) the evolution of $\dim(\sigma_{2,n})$ as a function of $n$ for some illustrative cases. We see that we have always $\sigma_{2,n} \ne \emptyset$ for a sufficiently large degree $n$, and that $\dim(\sigma_{2,n})$ increases when $n$ is increased. This suggests that $\sigma_2$ has an infinite (but countable) number of eigenvalues when $n \to \infty$. However, $\dim(\sigma_{2,n})$ is found to strongly depend on the geometry and the value of $N/\Omega_s$. More specifically, $\dim(\sigma_{2,n})$ is reduced when $N$ is increased, but the variation is more pronounced in the sphere than in pancake-like ellipsoids with $c/a \ll 1$. Finally, the number of eigenvalues seems to be bounded by $n(n+1)$, which corresponds to the number of spherical harmonics $Y_l^m$ of degree $2 \le l \le n+1$ and order $0 < |m| < l$. This shows that $\sigma_2$ is related to solid (spherical) harmonics.

Further insight into the link between $\sigma_2$ and the theory of spherical harmonics can be gained by considering the aligned case $1_\Omega = 1_z$ with $N = 2\Omega_s$. For this particular configuration, we have $\omega_- = \omega_+$ from equation (2.7). Then, pressure equation (2.6) and its associated BC for eigenvalues
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Spatial structure of the low-frequency (surface) modes for $u \in \mathcal{V}_{10}^0$ in the aligned case $1_{\Omega} = 1_z$ as a function of $N/\Omega_s$, in a flattened ellipsoid with $b/a = 1$ and $c/a = 0.1$. Meridional view of the absolute value of the normalised density perturbation $|\zeta|$ (in logarithmic scale) in the $(Oxz)$ plane (with rescaled axes to have a 1:1 ratio).</figcaption></figure>

$\omega \in \sigma_2$ reduce to

$$ \nabla^2 \Phi = 0, \quad n \cdot \nabla \Phi = -\frac{1}{i\omega} (2\Omega \times n) \cdot \nabla \Phi \quad \text{on } \partial V. \qquad (4.2a,b) $$

Laplace equation (4.2a) can then be solved using harmonic functions in spherical or ellipsoidal geometries, and the eigenvalues $\omega$ are the roots of a transcendental equation obtained from BC (4.2b). In a sphere, the pressure associated with $u \in \mathcal{V}_n^0$ is given in spherical coordinates $(r, \theta, \phi)$ by the solid harmonics $\Phi(r) = r^l Y_l^m(\theta, \phi)$ with $\omega = \pm 2\Omega_s m/l$, where $Y_l^m$ is the spherical harmonic of degree $l \le n + 1$ and azimuthal order $m < l$. The modes with $m > 0$ have a negative frequency $\omega < 0$ when $\Omega_s > 0$ (respectively positive $\omega > 0$ if $m < 0$). Given our phase convention $\propto \exp[i(\omega t + m\phi)]$, these modes have a positive azimuthal phase velocity $c_\phi = -\omega/m > 0$ when $\Omega_s > 0$ (i.e. prograde direction). The explicit solutions show that the spectrum $\sigma_{2,n}$ is here isomorph to the set of rational numbers $\mathbb{Q}$, such that $\dim(\sigma_{2,n}) = n(n + 1)$ in this case. Moreover, it can explicitly be shown that there is a one-to-one correspondence between the solutions in the sphere and those in the ellipsoid for this specific configuration [39]. Since $\mathbb{Q}$ is dense in $\mathbb{R}$, $\sigma_2$ belongs to the essential spectrum $\sigma_E$ when $1_{\Omega} = 1_z$ and $N = 2\Omega_s$.

## (ii) General properties

The above examples strongly suggest that $\sigma_2 \subset \sigma_E$. This can be shown by further looking at the pressure problem. Assuming that $|\omega| \neq \omega_\pm$, the pressure is governed by a boundary-value problem given by equation (2.6) in $V$ and a mixed-type BC on $\partial V$. This boundary-value scalar problem is elliptic if the equation is elliptic in $V$ (which occurs here when $0 < |\omega| < \omega_-$), but also if some invertibility conditions are satisfied by the BC on $\partial V$ [48] (see the ESM). If both requirements were met in $V$ and everywhere on $\partial V$ in a given interval $[\omega_1, \omega_2]$, then the spectrum would be discrete in this interval (see the ESM). Otherwise, the spectrum would be essential in such an interval. Here, we can prove that $\sigma_2$ belongs to the essential spectrum in an ellipsoid when $\Omega$ and $1_z$ are aligned (see the ESM). We also conjecture that $\sigma_2 \in \sigma_E$ in the misaligned case. Because the spectrum is pure point in ellipsoids, an infinite number of modes would thus exist with $0 < |\omega| < \omega_-$ when $n \to \infty$.

Finally, we explore the properties of the corresponding eigenvectors. We remind the reader that these modes only exist when $\omega_- \neq 0$. Otherwise, pure inertial modes (when $N = 0$), pure internal (gravity) modes (when $\Omega = 0$) or pure IGMs (when $\Omega \cdot 1_z = 0$) fill the entire interval $0 < |\omega| < \omega_+$. Two modal families are shown in figure 6. When $N/\Omega_s \le 2$, the energy of the lowest-frequency modes with $|\omega| \to 0$ is nearly aligned with the rotation axis. This structure is reminiscent of geostrophic modes without stratification, and fully agrees with asymptotic analysis
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Spatial structure of IGMs for $u \in V_{10}^0$ in the aligned case $1_\Omega = 1_z$ as a function of $N/\Omega_s$, in a flattened ellipsoid with $b/a = 1$ and $c/a = 0.1$. Meridional view of the absolute value of the normalised density perturbation $|\zeta|$ (in logarithmic scale) in the $(Oxz)$ plane (with rescaled axes to have a 1:1 ratio).</figcaption></figure>

when $N/\Omega_s \ll 1$. As investigated by Allen [38] when $\Omega \propto 1_z$, these modes appear because a weak stratification imposes some spatial restrictions on the steady geostrophic modes without stratification. The modes are nearly $z$-independent throughout the fluid (at the leading order), and have a phase travelling in the prograde direction. Note that this origin seems similar to the appearance of low-frequency Rossby modes in a volume $V$ without closed geostrophic contours [2]. On the contrary, the higher-frequency modes with $|\omega| \lesssim \omega_-$ have an energy that is maximum near the boundary and decays more or less significantly away from the boundary. These low-frequency (surface) modes seem similar to the (trapped) Kelvin waves in shallow-water models [67]. Indeed, as for the Kelvin waves, the wave energy of these low-frequency (surface) modes is mainly trapped near the boundary and decays away from $\partial V$ (see figure 6). Therefore, the presence of a boundary introduces a new subset of low-frequency (surface) modes in a rotating stratified fluid (which differ from the higher-frequency IGMs).

## (b) High-frequency IGMs

We investigate IGMs with angular frequencies belonging to $\sigma_1$. A few illustrative modes are shown in figure 7. IGMs are mainly inertial modes modified by gravity when $N \ll 2\Omega_s$, whereas they are gravity modes more or less affected by rotation when $N \gg 2\Omega_s$. In between these two limits, rotation and stratification can have competing effects. Moreover, some of their characteristics are directly controlled by the dispersion relation of IGWs. For example, the modes exhibit a pancake-like structure when $\omega/\Omega_s \to 2$ and $1_\Omega = 1_z$. From dispersion relation (2.8), the angle $\vartheta$ between $1_z$ and $k$ is such that $\vartheta \to 0$ when $\omega/\Omega_s \to 2$. Since the group velocity is orthogonal to $k$ for IGWs [3], the energy propagates in the horizontal direction (in agreement with the observed pancake-like structure). On the contrary, we have $\vartheta \to \pi/2$ for high-frequency IGWs with $\omega \to N$, such that the energy varies in the vertical direction. IGMs with $\omega \in \sigma_1$ are thus the direct counterparts in ellipsoids of IGWs in unbounded fluids.

Such as IGWs without boundaries, IGMs are dense in this frequency interval when $n \to \infty$. Yet, the eigenfrequencies are likely not equiprobable within this interval in ellipsoids. The frequency density is worth obtaining for physical applications, since it will strongly constrain the likelihood of obtaining (resonant) IGMs at a given frequency in geophysical vortices. The asymptotic density of the eigenvalues in $\sigma_1$ can be obtained using mathematical analysis when $n \to \infty$, and then compared to the density obtained from the polynomial computations at fixed values of $n$. We define the probability measure $\pi_n$ of the eigenvalues $\{\omega_j\}$ of $Q_\omega|W_n^0$ given by
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Density of IGMs in $\sigma_1$. Comparison between asymptotic density $d\pi_\infty/d\omega$ given in Proposition 4.2 and $\sigma_{2,35}$ obtained from polynomial computations. Left: Flattened ellipsoid $b/a = 1$ and $c/a = 0.1$ in the aligned case $\mathbf{1}_\Omega = \mathbf{1}_z$ with $N/\Omega_s = 5$. Right: Sphere $b/a = c/a = 1$ in the perpendicular case $\mathbf{1}_\Omega \cdot \mathbf{1}_z = 0$ with $N/\Omega_s = 10$.</figcaption></figure>

$\pi_n := (1/d_n) \sum_{j=1}^{d_n} \delta(\omega_j)$ with $d_n = 2N(W_n^0)$, where $\delta$ is the Dirac function. As explained above, the eigenvalues of SEP (3.3) cannot be easily separated in different subsets for every degree $n$. Therefore, we rather consider the joint repartition of the eigenvalues of $Q_\omega|V_n^0$. Naturally, the two repartitions have the same asymptotic distribution when $n \to \infty$. Asymptotic properties of the essential spectrum of PDEs can generally be investigated using microlocal analysis. Such an approach usually neglects the boundary effects but, here, it is possible to account for them while keeping a microlocal description. The asymptotic behaviour is given in Proposition 4.2.

**Proposition 4.2.** The measure $\pi_n$ converges weakly to the asymptotic measure $\pi_\infty$ when $n \to \infty$. The asymptotic measure is characterised by a cumulative distribution function given by

$$ \int_{-\infty}^{\lambda} d\pi_{\infty} = \frac{1}{8\pi} \mathrm{Area}(\mathcal{S}_{\lambda} \cap S^2) $$

where $S^2$ is the surface of the unit sphere, $\mathcal{S}_\lambda$ is the surface in $\mathbb{R}^3$ defined by

$$ \mathcal{S}_{\lambda} := \{\tilde{k} \in \mathbb{R}^3 \mid \omega(\tilde{k}) \le \lambda\}, \quad \tilde{k} = (\sqrt{A_1}k_x, \sqrt{A_2}k_y, \sqrt{A_3}k_z)^{\top}, $$

with $A_1 = 1/a^2$, $A_2 = 1/b^2$ and $A_3 = 1/c^2$, and where $\omega(\tilde{k})$ is the plane-wave dispersion relation given by equation (2.8) for the rescaled wave vector $\tilde{k}$. Hence, the asymptotic measure has a probability density function $f_\infty := d\pi_\infty/d\omega$ that is even and non-vanishing when $\omega_-^2 \le \omega^2 \le \omega_+^2$.

**Proof.** This is the most technical part of the paper. The proof closely follows the strategy presented in Colin de Verdière & Vidal [40] for the non-stratified problem, so we only outline the key steps below. The proof can be extended to matrix operators that only involve pseudo-differential operators of degree 0. This condition is not satisfied by the operator $\mathcal{T}$ defined in equation (2.13b). Thus, the starting point is to consider the matrix operator $\mathcal{L}$ defined in equation (2.13a). The fact that this operator is not self-adjoint is not a problem in the analysis (i.e. in the trace formula). Then, we apply the tools of microlocal analysis in the presence of boundary effects. They allow us to get the asymptotics of $\sum_{j=1}^{d_n} f(\omega_j)$ when $n \to \infty$ for any function $f \in C_0(\mathbb{R} \setminus 0)$ (i.e. the set of real-valued continuous functions vanishing near zero). The analysis heavily relies on the principal symbol of the operator $\mathcal{L}$ when $\omega \ne 0$. This gives the asymptotic measure in the interval $I_1$ (where the velocity equation is hyperbolic). Note that we avoid the value $\omega = 0$ because ker $\mathcal{L}$ is not clearly related to the steady solutions of the original physical problem (due to the existence of generalized eigenvectors). ■
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Asymptotic density $d\pi_\infty/d\tilde{\omega}$ of IGMs with normalised eigenvalues $\tilde{\omega} := \omega/\omega_+$, given by Proposition 4.2. Left: Sphere $b/a = c/a = 1$ in the aligned case $1_\Omega = 1_z$. Right: Axisymmetric ellipsoids with $b/a = 1$ and varying $c/a$ ratios. Top panel computed for $N/\Omega_s = 0.5$, and bottom one for $N/\Omega_s = 5$.</figcaption></figure>

Proposition 4.2 shows that $\sigma_1$ is dense in $I_1$ when $n \to \infty$. Next, we compare in figure 8 the asymptotic density with the density obtained from polynomial computations at a fixed degree $n=35$. The asymptotic density is normalised such that $\int_{-\infty}^{\infty} d\pi_\infty = 1$ but, since Proposition 4.2 heavily relies on the dispersion relation of IGWs, it is assumed that $d\pi_\infty/d\omega = 0$ when $\omega \in I_2$. To have a correct comparison, we have thus removed the spectra $\sigma_{2,n}$ before normalising the density. An excellent agreement is then obtained with the asymptotic density (even for the moderate polynomial degree $n=35$ considered here). We have also represented a particular configuration for which $1_\Omega \cdot 1_z = 0$, such that $\omega_- = 0$ from expression (2.7). In this case, IGMs are dense in the entire interval $]0, \omega_+[$. Note that a similar agreement is found for other configurations (not shown). Hence, Proposition 4.2 can be used to obtain the density of $\sigma_1$.

The asymptotic behaviour is further illustrated in figure 9 for different parameters and geometries. The density of pure inertial modes is known to be uniform within the frequency interval $]-2\Omega_s, 2\Omega_s[$ for a sphere [40]. Yet, the density of pure internal gravity modes without rotation is not uniform in the sphere (it is larger for the high-frequency modes than for the low-frequency modes). Between these two limits, the density is very close to the non-rotating case when $N/\Omega_s \gg 1$. On the contrary, the density of the high-frequency IGMs with $|\omega|/\omega_{max} \to 1$ approaches the density of the non-stratified modes when $N/\Omega_s \le 2$. We also observe that the geometry does modify the spectral density. When $N/\Omega_s < 2$, the density is maximal near $|\omega| \to \omega_-$ in elongated (i.e. prolate) ellipsoids with $c/a \gg 1$, whereas it is maximal near $|\omega| \to \omega_+$ in flattened (i.e. oblate) ellipsoids with $c/a \ll 1$. Instead, low-frequency IGMs with $|\omega| \to \omega_-$ are favoured in flattened ellipsoids when $N/\Omega_s > 2$.

Finally, the situation $\omega_+ = \omega_- \text{ occurring when } 1_\Omega = 1_z \text{ and } N/\Omega_s = 2$ deserves a specific consideration. The spectrum $\sigma_1$ is then reduced to the degenerate eigenvalue $\omega/\Omega_s = \pm 2$, whose multiplicity can be determined. Indeed, we know that $\dim(\sigma_{2,n}) = n(n+1)$ from the explicit solutions of equations (4.2a,b). Because the total number of non-zero eigenvalues is fixed for a degree $n$ (Proposition 3.3), we have $\dim(\sigma_{1,n}) = n(n+1)(n+2)/3$ when $\omega_+ = \omega_-.$

## 5. Physical discussion

The idealised model with a constant $N$ has allowed us to identify the key properties of the bounded oscillations in rotating stratified ellipsoids. Next, physical insight into geophysical vortices could be obtained by plugging realistic parameters into the model. Typical geophysical estimates for some large-scale vortices are given in table 1. These vortices are strongly stratified
<table><caption>**Table 1.** Typical estimates for some Jovian vortices and Mediterranean eddies (obtained from [10,11]). Rotation is usually described using the *f*-plane approximation with the Coriolis parameter $f = 2\Omega \cdot \mathbf{1}_z = 2\Omega_z \cos(\theta)$, where $\theta$ the mean colatitude of the vortex (figure 2).</caption><thead><tr><th></th><th><i>f</i> [rad.s<sup>-1</sup>]</th><th><i>N</i>/<i>f</i></th><th><i>a</i> [km]</th><th><i>b</i> [km]</th><th><i>c</i> [km]</th></tr></thead><tbody><tr><td>Jupiter's GRS</td><td>~ 1.4 × 10<sup>-4</sup></td><td><i>O</i>(10<sup>2</sup>)</td><td>~ 9000</td><td>~ 5000</td><td>~ 100</td></tr><tr><td>Jupiter's Oval BA</td><td>~ 1.9 × 10<sup>-4</sup></td><td><i>O</i>(10<sup>2</sup>)</td><td>~ 3000</td><td>~ 3000</td><td>~ 50</td></tr><tr><td>Mediterranean eddies</td><td>~ 8.5 × 10<sup>-5</sup></td><td>10–50</td><td>10–10<sup>2</sup></td><td>10–10<sup>2</sup></td><td>0.5–1</td></tr></tbody></table>

with $N/\Omega_s \gg 1$ and, consequently, they usually have a very small aspect ratio $c/a \le 10^{-2}$. However, the value of $N$ is still subject to strong uncertainties. Indeed, the remote determination of the ellipsoidal geometry is only sensitive to the difference between the BV frequency within the vortex and that of the ambient fluid [9–11]. The detection of IGWs (or IGMs) might be another way to estimate the internal stratification of such vortices. Various mechanisms could indeed excite such wave motions within a vortex (e.g. as interactions with the surrounding fluid, or bulk turbulence [68]). The bounds given by formula (2.7), which depend on the orientation of $\Omega$ and on the ratio $N/\Omega_s$, can be simplified in the relevant regime $N/\Omega_s \gg 1$. It gives

$$ \omega_{+}^{2} \approx N^{2} + 4 \left[ |\Omega|^{2} - (\Omega \cdot \mathbf{1}_{z})^{2} \right], \quad \omega_{-}^{2} \approx 4 |\Omega \cdot \mathbf{1}_{z}|^{2}. \qquad (5.1a,b) $$

We see that the upper bound approaches $\omega_+^2 \to N^2$ for realistic values $N/\Omega_s \gg 1$. Detecting IGMs within a vortex may thus yield an estimate of the vortex’s stratification.

A striking property of the global modes is that they admit, under our assumptions, polynomial solutions. This behaviour directly results from Proposition 3.1, and IGMs in other configurations can be very different. To illustrate this point, we consider a sphere of radius $R$ that is stratified under the radial gravity $g = -(g/R)r$. The background hydrostatic density is then $\nabla\rho_0 = -(\beta/R)r$ (where $\beta > 0$ is a constant). The BV frequency is thus given by $N^2 = N_0^2(r/R)^2$ using the Boussinesq approximation, where $N_0 = \sqrt{\beta g/\rho_*}$ is the constant value of $N$ at $r = R$. Then, the IGMs obey a mixed hyperbolic-elliptic equation that changes of type when [32,58]

$$ s^2/(\omega^2 - 4\Omega_s^2) + z^2/\omega^2 = 1/N_0^2, \qquad (5.2) $$

where $s = \sqrt{x^2 + y^2}$ is the cylindrical radius measured from the rotation axis $\Omega = \Omega_s \mathbf{1}_z$. For a given angular frequency $\omega$, equation (5.2) defines a turning surface that separates regions where the equation is hyperbolic and elliptic in $V$. Such turning surfaces strongly impact the spatial structure of IGMs, which are mainly trapped in the hyperbolic region [58]. Moreover, IGMs with $|\omega| \to 2\Omega_s < N$ become trapped in the equatorial waveguide when $N/\Omega_s \gg 1$ [69]. Finally, such turning surfaces are also likely responsible for a non-empty continuous spectrum without diffusion [58]. On the contrary, the eigenmodes in a stratified ellipsoid with a constant $N$ are always smooth and can penetrate deep inside the volume.

We have also shown the existence of low-frequency (surface) modes when $0 < |\omega| < \omega_-$*. Their properties are very similar to those of the (trapped) Kelvin waves in rotating shallow-water models [67], which owe their existence to topological features [70,71]. Our study brings here a complementary viewpoint. Indeed, we show that such low-frequency modes can exist in a bounded geometry with rigid walls if the pressure problem fails to be elliptic (due to the BC, see the ESM). Since this situation is expected to be rather generic, such low-frequency (surface) modes could also exist in more complicated models of stratified vortices.

Finally, we turn our attention to the lifetime of geophysical vortices. Without bulk motions, the slow temporal evolution of a stratified vortex is governed by diffusion processes that operate on the slow time scale [16,72] $\tau_\nu = (a^2/\nu)(2\Omega_s/N)^2$, where $\nu$ is the kinematic viscosity of the fluid within the vortex. Yet, it is known that bulk dissipation can be significantly increased in the presence of IGMs shaped by wave attractors [31,58]. Hence, we could wonder how much energy could be dissipated if some IGMs were forced inside an elliptical vortex. Without stratification, a
pure inertial mode (ω, u) satisfies in an ellipsoid the intriguing integral property [27,28]

$$ \int_V u^\dagger \cdot \nabla^2 u \, dV = \int_V u^\dagger \cdot \mathbb{L}(\nabla^2 u) \, dV = 0. \qquad (5.3) $$

This could be interpreted as a vanishing viscous dissipation of the inertial modes in an ellipsoid if there were no boundary layers. However, when viscous BCs are properly taken into account on the boundary (either with no-slip or stress-free BCs), the inertial modes have a non-zero viscous dissipation [2,73]. From a mathematical viewpoint, formula (5.3) is valid in an ellipsoid because (i) $\mathbb{L}(\nabla^2 u) \in \mathcal{V}_{n-2}^0$ when $u \in \mathcal{V}_n^0$ and (ii) two inertial eigenvectors ($u_i, u_j$) are orthogonal with respect to inner product (2.9) such that $\langle u_i, u_j \rangle = \delta_{ij}$ (when properly normalised). Here, the IGMs do not satisfy property (5.3) because they are not mutually orthogonal with respect to inner product (2.9) as explained in Proposition 2.2. This shows that the IGMs are (at least) subject to viscous dissipation in the bulk on the slow time scale $1/E$ in a stratified vortex. In addition, they will also be affected by thermal diffusion (if the stratification has a thermal origin, otherwise by molecular diffusion). Since the IGMs are here smooth polynomial solutions, it is possible to account for the leading-order dissipative effects using boundary-layer theory (see the ESM). This also shows that the IGMs would be dissipated on a typical time scale that is of the same order as $\tau_\nu$ in the absence of motions. The presence of IGMs may thus affect the slow temporal evolution of stratified vortices. Moreover, small-scale turbulence is often present inside stratified vortices, and IGMs could be involved in nonlinear interactions to generate (weakly) turbulent motions [19,68,74]. Therefore, IGMs may play a role in the transition towards turbulence and the long-term evolution of stratified vortices.

## 6. Concluding remarks

We have investigated the normal modes of a diffusionless fluid enclosed in a uniformly rotating ellipsoidal vortex, which is stratified in density under a background gravity. Indeed, they could be important for the dynamics of pancake-like stratified vortices (which are ubiquitous in geophysical flows). We have employed the Boussinesq approximation for simplicity, and shown that the problem can be recast as a mixed hyperbolic-elliptic equation for the fluid velocity. Next, we have simplified the model by considering a constant BV frequency. This idealised model retains the key ingredients to account for the dynamics of geophysical vortices. Under these assumptions, we have proved that the problem admits smooth polynomial eigenvectors in triaxial ellipsoids. Given the polynomial form of the eigenvectors, we have characterised the spectrum by using a bespoke numerical algorithm and asymptotic spectral theory following our prior study without stratification [40]. The normal modes consist of high-frequency IGMs and low-frequency (surface) modes. Using microlocal analysis, we have shown that the boundary plays a key role in sustaining such low-frequency (surface) modes. Finally, we have discussed our results in the light of geophysical applications, arguing that accounting for IGMs in the models could be useful to (better) understand the dynamics of stratified vortices.

Several questions have remained unanswered by our work (even in the linear theory). A more in-depth investigation of the properties of the low-frequency (surface) modes will be presented in a mathematical work. From a physical viewpoint, the model could also be extended to account for additional ingredients (while keeping polynomial solutions). For instance, a uniform-vorticity background flow $U_0$ will be taken into account to study the effects of the vortex's differential rotation (as measured by the Rossby number $Ro$). Large-scale vortices (e.g. Jupiter's vortices or Mediterranean eddies) are expected to be in the low-Rossby regime $|Ro| \lesssim \mathcal{O}(10^{-1})$, whereas smaller-scale vortices are in the opposite regime $|Ro| \gtrsim \mathcal{O}(1)$ (e.g. vortices at sub-mesoscales $\lesssim 1 - 10$ km in the Earth's oceans [6]). We have also entirely discarded magnetic effects. This is a reasonable assumption since, for realistic parameters, magnetic fields are expected to only weakly modify the high-frequency IGMs at the size of a vortex. However, magnetic fields are known to affect the low-frequency spectrum in rotating, stratified and electrically conducting fluids [21,75]. Moreover, there is a class of magnetic fields leaving invariant the three-dimensional
polynomial vector spaces in an ellipsoid [76,77]. Therefore, our ellipsoidal model could be used as a toy model to investigate the properties of these Magneto-Archimedean-Coriolis (MAC) modes, by combining microlocal analysis and polynomial computations. Other physical ingredients (e.g. a variable BV frequency) would break the exact polynomial nature of the eigenvectors. Hence, the spectrum should be investigated using another numerical method.

Finally, future studies could build upon our work to understand the onset of turbulence within geophysical vortices. Turbulent mixing and dissipation may indeed (partly) explain why some stratified vortices are remarkably stable over long time scales, while some others have much shorter lifetimes. Several routes exist for the excitation of turbulence in stratified fluids, but direct applications of our work include the study of elliptical instabilities. So far, the flows driven by elliptical instabilities have only received scant attention in stratified environments [78–81]. The polynomial Galerkin algorithm presented in this work can be extended in the linear theory to determine the onset of the instabilities upon an elliptical flow. Therefore, elliptical instabilities will be considered in future work to explore, as a function of the Rossby number, the transition towards turbulence in stratified vortices.

**Data Accessibility.** The paper has an Electronic Supplementary Material (ESM). The source code is available at https://bitbucket.org/vidalje/.

**Authors' Contributions.** This work is an idea of JV to pursue the interdisciplinary collaboration with YCdV after their joint work [40]. JV designed the study and performed the Galerkin computations, whereas YCdV conducted the microlocal analysis. Both authors discussed the results presented in the paper and drafted the manuscript before giving final approval for submission.

**Competing Interests.** The authors declare that they have no competing interests.

**Declaration of AI use.** The authors did not use AI-assisted technologies in creating this article.

**Funding.** JV received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 847433, THEIA project). The numerical computations were performed using both laboratory resources and using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. A CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant's open access conditions.

**Acknowledgements.** JV would like to dedicate this work to the memory of Norman R. Lebovitz (1935-2022) and Stjin Vantieghem (1983-2023), who both revivified the topic and made pioneering contributions. JV acknowledges Bernard Valette and David Cébron for fruitful discussions over the years. YCdV acknowledges Gerd Grubb for her clarification of the definition of a boundary-value elliptic problem. The authors acknowledge the three anonymous referees, whose valuable comments helped them to improve the clarity and quality of the manuscript.

## References

1. Mowbray DE, Rarity BSH. 1967 A theoretical and experimental investigation of the phase configuration of internal waves of small amplitude in a density stratified liquid. *J. Fluid Mech.* **28**, 1–16. (doi:10.1017/S0022112067001867).

2. Greenspan HP. 1968 *The theory of rotating fluids*. Cambridge, UK: Cambridge University Press.

3. LeBlond PH, Mysak LA. 1981 *Waves in the ocean*. Amsterdam: Elsevier.

4. Dauxois T, Joubaud S, Odier P, Venaille A. 2018 Instabilities of internal gravity wave beams. *Annu. Rev. Fluid Mech.* **50**, 131–156. (doi:10.1146/annurev-fluid-122316-044539).

5. Le Bars M. 2023 Wave turbulence in geophysical flows. *J. Fluid Mech.* **962**, F1. (doi:10.1017/jfm.2023.219).

6. McWilliams JC. 2016 Submesoscale currents in the ocean. *Proc. R. Soc. A* **472**, 20160117. (doi:10.1098/rspa.2016.0117).

7. Hopfinger EJ, Van Heijst GJF. 1993 Vortices in rotating fluids. *Annu. Rev. Fluid Mech.* **25**, 241–289. (doi:10.1146/annurev.fl.25.010193.001325).

8. Billant P, Chomaz JM. 2001 Self-similarity of strongly stratified inviscid flows. *Phys. Fluids* **13**, 1645–1651. (doi:10.1063/1.1369125).
9. Hassanzadeh P, Marcus PS, Le Gal P. 2012 The universal aspect ratio of vortices in rotating stratified flows: theory and simulation. *J. Fluid Mech.* **706**, 46–57. (doi:10.1017/jfm.2012.180).

10. Aubert O, Le Bars M, Le Gal P, Marcus PS. 2012 The universal aspect ratio of vortices in rotating stratified flows: experiments and observations. *J. Fluid Mech.* **706**, 34–45. (doi:10.1017/jfm.2012.176).

11. Lemasquerier D, Facchini G, Favier B, Le Bars M. 2020 Remote determination of the shape of Jupiter's vortices from laboratory experiments. *Nat. Phys.* **16**, 695–700. (doi:10.1038/s41567-020-0833-9).

12. Richardson PL, Bower AS, Zenk W. 2000 A census of Meddies tracked by floats. *Prog. Oceanogr.* **45**, 209–250. (doi:10.1016/S0079-6611(99)00053-1).

13. Marcus PS. 1993 Jupiter's Great Red Spot and other vortices. *Annu. Rev. Astron. Astrophys.* **31**, 523–569. (doi:10.1146/annurev.aa.31.090193.002515).

14. McKiver WJ. 2015 The ellipsoidal vortex: A novel approach to geophysical turbulence. *Adv. Math. Phys.* **2015**. (doi:10.1155/2015/613683).

15. Tsang YK, Dritschel DG. 2015 Ellipsoidal vortices in rotating stratified fluids: beyond the quasi-geostrophic approximation. *J. Fluid Mech.* **762**, 196–231. (doi:10.1017/jfm.2014.630).

16. Facchini G, Le Bars M. 2016 On the lifetime of a pancake anticyclone in a rotating stratified flow. *J. Fluid Mech.* **804**, 688–711. (doi:10.1017/jfm.2016.549).

17. Sipp D, Lauga E, Jacquin L. 1999 Vortices in rotating systems: centrifugal, elliptic and hyperbolic type instabilities. *Phys. Fluids* **11**, 3716–3728. (doi:10.1063/1.870180).

18. Yim E, Billant P, Ménesguen C. 2016 Stability of an isolated pancake vortex in continuously stratified-rotating fluids. *J. Fluid Mech.* **801**, 508–553. (doi:10.1017/jfm.2016.402).

19. Kerswell RR. 2002 Elliptical instability. *Annu. Rev. Fluid Mech.* **34**, 83–113. (doi:10.1146/annurev.fluid.34.081701.171829).

20. Miyazaki T, Fukumoto Y. 1992 Three-dimensional instability of strained vortices in a stably stratified fluid. *Phys. Fluids* **4**, 2515–2522. (doi:10.1063/1.858438).

21. Kerswell RR. 1993 Elliptical instabilities of stratified, hydromagnetic waves. *Geophys. Astrophys. Fluid Dyn.* **71**, 105–143. (doi:10.1080/03091929308203599).

22. Le Bars M, Le Dizès S. 2006 Thermo-elliptical instability in a rotating cylindrical shell. *J. Fluid Mech.* **563**, 189–198. (doi:10.1017/S0022112006001674).

23. Vidal J, Cébron D, Ud-Doula A, Alecian E. 2019 Fossil field decay due to nonlinear tides in massive binaries. *Astron. Astrophys.* **629**, A142. (doi:10.1051/0004-6361/201935658).

24. Cartan ME. 1922 Sur les petites oscillations d'une masse fluide. *Bull. Sci. Math.* **46**, 317–369.

25. Poincaré H. 1885 Sur l'équilibre d'une masse fluide animée d'un mouvement de rotation. *Acta Math.* **7**, 259–380.

26. Rieutord M, Georgeot B, Valdettaro L. 2000 Wave attractors in rotating fluids: a paradigm for ill-posed Cauchy problems. *Phys. Rev. lett.* **85**, 4277. (doi:10.1103/PhysRevLett.85.4277).

27. Vantieghem S. 2014 Inertial modes in a rotating triaxial ellipsoid. *Proc. R. Soc. A* **470**, 20140093. (doi:10.1098/rspa.2014.0093).

28. Ivers DJ. 2017 Enumeration, orthogonality and completeness of the incompressible Coriolis modes in a tri-axial ellipsoid. *Geophys. Astrophys. Fluid Dyn.* **111**, 333–354. (doi:10.1080/03091929.2017.1330412).

29. Backus G, Rieutord M. 2017 Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid. *Phys. Rev. E* **95**, 053116. (doi:10.1103/PhysRevE.95.053116).

30. Veronis G. 1970 The analogy between rotating and stratified fluids. *Annu. Rev. Fluid Mech.* **2**, 37–66. (doi:10.1146/annurev.fl.02.010170.000345).

31. Rieutord M, Noui K. 1999 On the analogy between gravity modes and inertial modes in spherical geometry. *Eur. Phys. J. B* **9**, 731–738. (doi:10.1007/s100510050818).

32. Friedlander S, Siegmann WL. 1982 Internal waves in a rotating stratified fluid in an arbitrary gravitational field. *Geophys. Astrophys. Fluid Dyn.* **19**, 267–291. (doi:10.1080/03091928208208959).

33. Brunet M, Dauxois T, Cortet PP. 2019 Linear and nonlinear regimes of an inertial wave attractor. *Phys. Rev. Fluids* **4**, 034801. (doi:10.1103/PhysRevFluids.4.034801).

34. Favier B, Le Dizès S. 2024 Inertial wave super-attractor in a truncated elliptic cone. *J. Fluid Mech.* **980**, A6. (doi:10.1017/jfm.2024.5).

35. Maas LRM, Benielli D, Sommeria J, Lam FPA. 1997 Observation of an internal wave attractor in a confined, stably stratified fluid. *Nature* **388**, 557–561. (doi:10.1038/41509).

36. Pacary C, Dauxois T, Ermanyuk E, Metz P, Moulin M, Joubaud S. 2023 Observation of inertia-gravity wave attractors in an axisymmetric enclosed basin. *Phys. Rev. Fluids* **8**, 104802. (doi:10.1103/PhysRevFluids.8.104802).
37. Colin De Verdière Y, Saint-Raymond L. 2020 Attractors for two-dimensional waves with homogeneous Hamiltonians of degree 0. *Commun. Pure Appl. Math.* **73**, 421–462. (doi:10.1002/cpa.21845).

38. Allen JS. 1971 Some aspects of the initial-value problem for the inviscid motion of a contained, rotating, weakly-stratified fluid. *J. Fluid Mech.* **46**, 1–23. (doi:10.1017/S0022112071000363).

39. Friedlander S, Siegmann WL. 1982 Internal waves in a contained rotating stratified fluid. *J. Fluid Mech.* **114**, 123–156. (doi:10.1017/S002211208200007X).

40. Colin de Verdière Y, Vidal J. 2023 The spectrum of the Poincaré operator in an ellipsoid. *arxiv:2305.01369*. (doi:10.48550/arXiv:2305.01369).

41. Spiegel EA, Veronis G. 1960 On the Boussinesq approximation for a compressible fluid. *Astrophys. J.* **131**, 442–447. (doi:10.1086/146849).

42. Noir J, Cébron D. 2013 Precession-driven flows in non-axisymmetric ellipsoids. *J. Fluid Mec.* **737**, 412–439. (doi:10.1017/jfm.2013.524).

43. Lam H, Delache A, Godeferd FS. 2023 Supply mechanisms of the geostrophic mode in rotating turbulence: interactions with self, waves and eddies. *J. Fluid Mech.* **971**, A10. (doi:10.1017/jfm.2023.644).

44. Le Reun T, Favier B, Barker AJ, Le Bars M. 2017 Inertial wave turbulence driven by elliptical instability. *Phys. Rev. Lett.* **119**, 034502. (doi:10.1103/PhysRevLett.119.034502).

45. Park J, Billant P. 2013 Instabilities and waves on a columnar vortex in a strongly stratified and rotating fluid. *Phys. Fluids* **25**. (doi:10.1063/1.4816512).

46. Le Dizès S. 2008 Inviscid waves on a Lamb-Oseen vortex in a rotating stratified fluid: consequences for the elliptic instability. *J. Fluid Mech.* **597**, 283–303. (doi:10.1017/S0022112007009780).

47. Friedlander S. 1985 Internal oscillations in the Earth's fluid core. *Geophys. J. Int.* **80**, 345–361. (doi:10.1111/j.1365-246X.1985.tb05099.x).

48. Grubb G. 2008 *Distributions and operators*. Graduate Texts in Mathematics (USA): Springer.

49. Faure F. 2023 Manifestation of the topological index formula in quantum waves and geophysical waves. *Ann. Henri Lebesgue* **6**, 449–492. (doi:10.5802/ahl.169).

50. Venaille A, Onuki Y, Perez N, Leclerc A. 2023 From ray tracing to waves of topological origin in continuous media. *SciPost Phys.* **14**, 062. (doi:10.1137/14.4.062).

51. Colin de Verdière Y. 2020 Spectral theory of pseudodifferential operators of degree 0 and an application to forced linear waves. *Anal. PDE* **13**, 1521–1537. (doi:10.2140/apde.2020.13.1521).

52. Leray J. 1934 Sur le mouvement d'un liquide visqueux emplissant l'espace. *Acta Math.* **63**, 193–248.

53. Foias C, Manley O, Rosa R, Temam R. 2001 *Navier-Stokes equations and turbulence*. Cambridge, UK: Cambridge University Press.

54. Teed RJ, Dormy E. 2023 Solenoidal force balances in numerical dynamos. *J. Fluid Mech.* **964**, A26. (doi:10.1017/jfm.2023.332).

55. Barston EM. 1967 Eigenvalue Problem for Lagrangian Systems. II. *J. Math. Phys.* **8**, 1886–1892. (doi:10.1063/1.1705433).

56. Valette B. 1989 Étude d'une classe de problèmes spectraux. *C. R. Acad. Sci. Paris* **309**, 785–788.

57. Wouk A. 1966 A note on square roots of positive operators. *SIAM Rev.* **8**, 100–102. (doi:10.1137/1008008).

58. Dintrans B, Rieutord M, Valdettaro L. 1999 Gravito-inertial waves in a rotating stratified sphere or spherical shell. *J. Fluid Mech.* **398**, 271–297. (doi:10.1017/S0022112099006308).

59. Barston EM. 1967 Eigenvalue problem for Lagrangian systems. *J. Math. Phys.* **8**, 523–532. (doi:10.1063/1.1705227).

60. Axler S, Shin PJ. 2018 The Neumann problem on ellipsoids. *J. Appl. Math. Comput.* **57**, 261–278. (doi:10.1007/s12190-017-1105-4).

61. Rabitti A, Maas LRM. 2014 Inertial wave rays in rotating spherical fluid domains. *J. Fluid Mech.* **758**, 621–654. (doi:10.1017/jfm.2014.551).

62. Ivers DJ, Jackson A, Winch D. 2015 Enumeration, orthogonality and completeness of the incompressible Coriolis modes in a sphere. *J. Fluid Mech.* **766**, 468–498. (doi:10.1017/jfm.2015.27).

63. Lebovitz NR. 1989 The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases. *Geophys. Astrophys. Fluid Dyn.* **46**, 221–243. (doi:10.1080/03091928908208913).
64. Vidal J, Su S, Cébron D. 2020 Compressible fluid modes in rigid ellipsoids: towards modal acoustic velocimetry. *J. Fluid Mech.* **885**, A39. (doi:10.1017/jfm.2019.1004).

65. Vidal J, Cébron D. 2020 Acoustic and inertial modes in planetary-like rotating ellipsoids. *Proc. R. Soc. A* **476**, 20200131. (doi:10.1098/rspa.2020.0131).

66. Peacock T, Weidman P. 2005 The effect of rotation on conical wave beams in a stratified fluid. *Exp. Fluids* **39**, 32–37. (doi:10.1007/s00348-005-0955-y).

67. Thomson W. 1880 On gravitational oscillations of rotating water. *Proc. R. Soc. Edinburgh* **10**, 92–100. (doi:10.1017/S0370164600043467).

68. Lanchon N, Mora DO, Monsalve E, Cortet PP. 2023 Internal wave turbulence in a stratified fluid with and without eigenmodes of the experimental domain. *Phys. Rev. Fluids* **8**, 054802. (doi:10.1103/PhysRevFluids.8.054802).

69. Stewartson K, Walton IC. 1976 On waves in a thin shell of stratified rotating fluid. *Proc. R. Soc. Lond. A* **349**, 141–156. (doi:10.1098/rspa.1976.0064).

70. Tauber C, Delplace P, Venaille A. 2019 A bulk-interface correspondence for equatorial waves. *J. Fluid Mech.* **868**, R2. (doi:10.1017/jfm.2019.233).

71. Venaille A, Delplace P. 2021 Wave topology brought to the coast. *Phys. Rev. Res.* **3**, 043002. (doi:10.1103/PhysRevResearch.3.043002).

72. Le Bars M. 2021 Numerical study of the McIntyre instability around Gaussian floating vortices in thermal wind balance. *Phys. Rev. Fluids* **6**, 093801. (doi:10.1103/PhysRevFluids.6.093801).

73. Vidal J, Cébron D. 2023 Precession-driven flows in stress-free ellipsoids. *J. Fluid Mech.* **954**, A5. (doi:10.1017/jfm.2022.976).

74. Boury S, Maurer P, Joubaud S, Peacock T, Odier P. 2023 Triadic resonant instability in confined and unconfined axisymmetric geometries. *J. Fluid Mech.* **957**, A20. (doi:10.1017/jfm.2023.58).

75. Friedlander S. 1987 Hydromagnetic waves in the Earth’s fluid core. *Geophys. Astrophys. Fluid Dyn.* **39**, 315–333. (doi:10.1080/03091928708208816).

76. Malkus WVR. 1967 Hydromagnetic planetary waves. *J. Fluid Mech.* **28**, 793–802. (doi:10.1017/S0022112067002447).

77. Gerick F, Jault D, Noir J, Vidal J. 2020 Pressure torque of torsional Alfvén modes acting on an ellipsoidal mantle. *Geophys. J. Int.* **222**, 338–351. (doi:10.1093/gji/ggaa166).

78. Cébron D, Maubert P, Le Bars M. 2010 Tidal instability in a rotating and differentially heated ellipsoidal shell. *Geophys. J. Int.* **182**, 1311–1318. (doi:10.1111/j.1365-246X.2010.04712.x).

79. Le Reun T, Favier B, Le Bars M. 2018 Parametric instability and wave turbulence driven by tidal excitation of internal waves. *J. Fluid Mech.* **840**, 498–529. (doi:10.1017/jfm.2018.18).

80. Vidal J, Cébron D, Schaeffer N, Hollerbach R. 2018 Magnetic fields driven by tidal mixing in radiative stars. *Mon. Not. R. Astron. Soc.* **475**, 4579–4594. (doi:10.1093/mnras/sty080).

81. Onuki Y, Joubaud S, Dauxois T. 2023 Breaking of internal waves parametrically excited by ageostrophic anticyclonic instability. *J. Phys. Oceanogr.* **53**, 1591–1613. (doi:10.1175/JPO-D-22-0152.1).
Electronic Supplementary Material of "Inertia-gravity waves
in geophysical vortices"

Jérémie Vidal* & Yves Colin de Verdière†

February 14, 2024

Abstract

We present the derivation of the boundary-value pressure problem governing the nor-
mal modes of an incompressible fluid, which is rotating and stratified in density under the
Bousinesq approximation. Next, we discuss for which conditions this is an elliptic boundary-
value problem (using microlocal analysis). Applications to ellipsoidal vortices with a constant
Brunt-Väisälä frequency are then considered. Finally, we introduce a boundary-layer analysis
to determine the (leading-order) diffusive decay rates of the normal modes in an ellipsoid.

1 Boundary-value pressure problem

We consider a diffusionless and incompressible fluid, which is enclosed in a smooth bounded
domain $V$ (whose boundary is denoted by $\partial V$). The fluid is rotating at the angular velocity
$\Omega = (\Omega_x, \Omega_y, \Omega_z)^\top$ with respect to an inertial frame, and is stratified in density under the imposed
gravity field $g = -g(z)\mathbf{1}_z$ within the Boussinesq approximation. At rest, the pressure $P_0$ and the
density satisfy the hydrostatic balance $\nabla P_0 = [\rho_* + \rho_0(z)]g$, where $\rho_*$ is the typical (mean) den-
sity of the fluid within the Boussinesq approximation. Then, we seek small-amplitude oscillatory
perturbations (where $\omega \in \mathbb{R}$ is the angular frequency) for the velocity $\mathbf{u} \exp(i\omega t)$, the pressure
$\Phi \exp(i\omega t)$, and the density $\zeta \exp(i\omega t)$ upon the hydrostatic basic state. These perturbations are
governed in the rotating frame by the incompressible Euler equations

$$
i\omega u + 2(\Omega \times u) = -\nabla\Phi - (\zeta/\rho^*)g(z)\mathbf{1}_z, \quad \nabla \cdot u = 0, \qquad (1a,b)
$$

and the buoyancy equation

$$
i\omega\zeta = \frac{\rho_* N^2(z)}{g(z)} (\mathbf{u} \cdot \mathbf{1}_z), \quad N^2(z) := -\frac{g(z)}{\rho_*} \frac{d\rho_0}{dz}, \qquad (2a,b)
$$

where $N(z)$ is the Brunt-Väisälä (BV) frequency. The latter equations can be reduced, as long as
$\omega \neq 0$, to a quadratic equation in $\omega$ for the velocity given by

$$
-\omega^2 \mathbf{u} + i\omega (2\Omega \times \mathbf{u}) + N^2(z) (\mathbf{u} \cdot \mathbf{1}_z) \mathbf{1}_z = -i\omega \nabla \Phi. \quad (3)
$$

<footnote>*Université Grenoble Alpes, CNRS, ISTerre, 38000 Grenoble, France</footnote>

<footnote>†Université Grenoble Alpes, CNRS, Institut Fourier, 38000 Grenoble, France</footnote>
BOX 1.1: Poincaré problem

When $N = 0$, the velocity is given by (see formula 2.7.2 in Greenspan [1])

$$ (4|\Omega|^2 - \omega^2)u = -i\omega\nabla\Phi - \frac{4}{i\omega}(\Omega \cdot \nabla\Phi)\Omega + 2\Omega \times \nabla\Phi. \quad (B1.1.1) $$

The divergence of the latter equation yields the so-called Poincaré equation [1]

$$ \omega^2 \nabla^2 \Phi = (2\Omega \cdot \nabla)^2 \Phi. \quad (B1.1.2) $$

The associated BC, obtained from equation (15) below, is given by

$$ \omega^2 \nabla \Phi \cdot n = i\omega(2\Omega \times n) \cdot \nabla \Phi + 4(\Omega \cdot \nabla \Phi)(\Omega \cdot n) \quad \text{on } \partial V, \quad (B1.1.3) $$

which is identical to equation (2.7.4) in Greenspan [1].

Finally, equation (3) is supplemented with the no-penetration boundary condition (BC) given by $u \cdot n|_{\partial V} = 0$ (i.e. $u$ is tangent to $\partial V$), where $n$ is the unit vector normal to the boundary.

Our goal is to find an alternative formulation of the problem in terms of the pressure $\Phi$. An equation for the pressure can be obtained by taking the divergence of equation (3), which gives

$$ -i\omega\nabla^2\Phi = i\omega\nabla \cdot (2\Omega \times u) + (u \cdot 1_z)\partial_z N^2 + N^2(z)\partial_z(u \cdot 1_z). \quad (4) $$

The associated boundary condition (BC) is then obtained by taking the scalar product of equation (3) with $n$, which gives

$$ -i\omega(\nabla\Phi \cdot n) = i\omega(2\Omega \times u) \cdot n + N^2(z)(u \cdot 1_z)(1_z \cdot n) \quad \text{on } \partial V. \quad (5) $$

However, equations (4)-(5) still depend on $u$. Therefore, an expression of $u$ in terms of $\Phi$ only is required before obtaining the equation for the pressure (e.g. the Poincaré problem, see Box 1.1). We will give below the explicit formula(s) accounting for misaligned rotation and gravity with a varying BV frequency $N(z)$.

## 1.1 Second-order equation in the volume

To find the equation satisfied by the pressure $\Phi$ in $V$, we first find the equation relating $u$ and $\Phi$ (Proposition 1.1). Then, we enforce the incompressibility condition to obtain the pressure equation (Proposition 1.2). Finally, we give in Proposition 1.3 the pressure equation when $N$ is constant. It will be of prime importance for the microlocal analysis presented in Section 2.

**Proposition 1.1.** The velocity $u$ is related to the pressure $\Phi$ by the equation $D_{\omega}(z)u = f$ when $D_{\omega}(z) = (\omega^2 - \omega_+^2)(\omega^2 - \omega_-^2) \neq 0$ with

$$ \omega_{\pm}^2 = \frac{1}{2}(N^2(z) + 4|\Omega|^2) \pm \sqrt{\frac{1}{4}(N^2(z) + 4|\Omega|^2)^2 - 4N^2(z)(\Omega \cdot 1_z)^2}, \quad (6) $$

and with

$$ f = i\omega^3 \nabla\Phi - 2\omega^2(\Omega \times \nabla\Phi) - 4i\omega(\Omega \cdot \nabla\Phi)\Omega + N^2(z) [(2\Omega \cdot 1_z)(1_z \times \nabla\Phi) + i\omega 1_z \times (1_z \times \nabla\Phi)]. \quad (7) $$
Proof. We recast Euler equation (1a) and buoyancy equation (2) as a linear system given by

$$ M \begin{pmatrix} u_x \\ u_y \\ u_z \\ \zeta/\rho_* \end{pmatrix} = - \begin{pmatrix} \partial_x \Phi \\ \partial_y \Phi \\ \partial_z \Phi \\ 0 \end{pmatrix}, \quad M := \begin{pmatrix} i\omega & -2\Omega_z & 2\Omega_y & 0 \\ 2\Omega_z & i\omega & -2\Omega_x & 0 \\ -2\Omega_y & 2\Omega_x & i\omega & g(z) \\ 0 & 0 & -N^2(z)/g(z) & i\omega \end{pmatrix}, \quad (8a,b) $$

and with $\mathbf{u} = (u_x, u_y, u_z)^\top$. The system admits non-zero solutions when $\det(M) := \mathcal{D}_\omega(z) \neq 0$. In this case, we can invert the linear system and the result follows. Note that it reduces to equation (2.7.2) in Greenspan [1] when $N = 0$. ■

**Proposition 1.2.** When $\mathcal{D}_\omega(z) \neq 0$, the pressure obeys the second-order equation given by

$$ \nabla \cdot (\mathbf{f}/\mathcal{D}_\omega) = 0 \iff \nabla \cdot \mathbf{f} - (1/\mathcal{D}_\omega) f_z \partial_z D_\omega = 0, \quad (9) $$

$$ \text{with } \partial_z \mathcal{D}_\omega = -(\omega^2 - \omega_-^2)\partial_z \omega_+^2 - (\omega^2 - \omega_+^2)\partial_z \omega_-^2. $$

When $N(z)$ is variable and $\Omega = \Omega_z \mathbf{1}_z$, equation (9) reduces to

$$ \nabla^2 \Phi + \frac{\partial}{\partial z} \left( \frac{N^2 - 4\Omega_z^2}{\omega^2 - N^2} \frac{\partial \Phi}{\partial z} \right) = 0, \quad (10) $$

which is identical to equation (5.1) in Friedlander & Siegmann [2]. For oscillatory motions with $\omega \neq 0$, equation (11) is thus equivalent to equation (3) when $|\omega| \neq \omega_\pm(z)$. However, the pressure must obey another equation $\mathcal{D}_\omega(z) = 0$ (see below).

**Proposition 1.3.** When $N$ is constant, equation (9) reduces to $\mathcal{P}_\omega(\Phi) = 0$ with

$$ \mathcal{P}_\omega := -\nabla^2 - \frac{N^2}{\omega^2 - N^2} (\mathbf{1}_z \cdot \nabla)^2 + \frac{4}{\omega^2 - N^2} (\mathbf{\Omega} \cdot \nabla)^2. \quad (11) $$

The operator $\mathcal{P}_\omega$ is elliptic when $0 < |\omega| < |\omega_\pm|$, and its principal symbol $\mathbf{p}$ is given by

$$ \mathbf{p} := Ak_x^2 + Bk_y^2 + Ck_z^2 + 2Dk_yk_z + 2Ek_xk_z + 2Fk_xky \quad (12) $$

with

$$ \begin{aligned} A &= \frac{\omega^2 - N^2 - 4\Omega_x^2}{\omega^2 - N^2}, & B &= \frac{\omega^2 - N^2 - 4\Omega_y^2}{\omega^2 - N^2}, & C &= \frac{\omega^2 - 4\Omega_z^2}{\omega^2 - N^2}, \\ D &= \frac{-4\Omega_y\Omega_z}{\omega^2 - N^2}, & E &= \frac{-4\Omega_x\Omega_z}{\omega^2 - N^2}, & F &= \frac{-4\Omega_x\Omega_y}{\omega^2 - N^2}. \end{aligned} $$

Finally, it is worth explicitly writing the Cartesian components of $\mathbf{f} = (f_x, f_y, f_z)^\top$ when $N$ is constant (as we will reuse them below). They are given by

$$ \frac{f_x}{\omega^2 - N^2} = i\omega[A\partial_x\Phi + F\partial_y\Phi + E\partial_z\Phi] - \frac{2\omega^2\Omega_y}{\omega^2 - N^2}\partial_z\Phi \quad + 2\Omega_z\partial_y\Phi, \quad (13a) $$

$$ \frac{f_y}{\omega^2 - N^2} = i\omega[F\partial_x\Phi + B\partial_y\Phi + D\partial_z\Phi] + \frac{2\omega^2\Omega_x}{\omega^2 - N^2}\partial_z\Phi \quad - 2\Omega_z\partial_x\Phi, \quad (13b) $$

$$ \frac{f_z}{\omega^2 - N^2} = i\omega[E\partial_x\Phi + D\partial_y\Phi + C\partial_z\Phi] + \frac{2\omega^2}{\omega^2 - N^2}(\Omega_y\partial_x\Phi - \Omega_x\partial_y\Phi). \quad (13c) $$

where the coefficients $(A, B, C, D, E, F)$ appear in the principal symbol of $\mathcal{P}_\omega$ (Proposition 1.3).
When $\Omega = \Omega_z \mathbf{1}_z$ and $N$ is constant, pressure equation (10) simplifies into

$$ -\nabla^2 \Phi - \frac{N^2 - 4\Omega_z^2}{\omega^2 - N^2} (\mathbf{1}_z \cdot \nabla)^2 \Phi = 0, \qquad (B1.2.1) $$

which is identical to equation (2.21) in Friedlander & Siegmann [3]. If required, the velocity is related to the pressure as

$$ [4\Omega_z^2 - \omega^2]\mathbf{u} = -i\omega\nabla\Phi + \frac{i\omega}{\omega^2 - N^2}[4\Omega_z^2 - N^2](\mathbf{1}_z \cdot \nabla\Phi)\mathbf{1}_z + 2\Omega_z(\mathbf{1}_z \times \nabla\Phi). \quad (B1.2.2) $$

To find the pressure BC, we rearrange equation (15) as

$$ i\omega[\omega^2 - N^2]\nabla\Phi \cdot \mathbf{n} = 2\Omega_z[\omega^2 - N^2](\mathbf{1}_z \times \nabla\Phi) \cdot \mathbf{n} + i\omega[4\Omega_z^2 - N^2](\mathbf{1}_z \cdot \mathbf{n})\frac{\partial\Phi}{\partial z}, $$

which simplifies into

$$ \nabla\Phi \cdot \mathbf{n} = \frac{2\Omega_z}{i\omega}(\mathbf{n} \times \mathbf{1}_z) \cdot \nabla\Phi - \left[\frac{N^2 - 4\Omega_z^2}{\omega^2 - N^2}\right](\mathbf{1}_z \cdot \mathbf{n})\frac{\partial\Phi}{\partial z} \quad \text{on } \partial V. \quad (B1.2.3) $$

BC (B1.2.3) is identical to formula (2.22) in Friedlander & Siegmann [3].

## 1.2 Boundary condition

**Proposition 1.4.** The pressure BC is given by $f \cdot n|_{\partial V} = 0$, where the components of $f$ are given in equations (13a-c).

For the sake of comparison with prior studies, the pressure BC can be rewritten as

$$ i\omega^3 \nabla\Phi \cdot \mathbf{n} = 2\omega^2(\Omega \times \nabla\Phi) \cdot \mathbf{n} + 4i\omega(\Omega \cdot \nabla\Phi)(\Omega \cdot \mathbf{n}) \\ - N^2(z) \left[ (2\Omega \cdot \mathbf{1}_z)(\mathbf{1}_z \times \nabla\Phi) + i\omega \mathbf{1}_z \times (\mathbf{1}_z \times \nabla\Phi) \right] \cdot \mathbf{n}. \quad (14) $$

Since we have $\mathbf{1}_z \times (\mathbf{1}_z \times \nabla\Phi) = (\mathbf{1}_z \cdot \nabla\Phi)\mathbf{1}_z - \nabla\Phi$, we can rewrite the pressure BC as

$$ i\omega[\omega^2 - N^2(z)]\nabla\Phi \cdot \mathbf{n} = 2[\omega^2(\Omega \times \nabla\Phi) - N^2(z)(\Omega \cdot \mathbf{1}_z)(\mathbf{1}_z \times \nabla\Phi)] \cdot \mathbf{n} \\ + i\omega[4(\Omega \cdot \nabla\Phi)\Omega - N^2(z)(\mathbf{1}_z \cdot \nabla\Phi)\mathbf{1}_z] \cdot \mathbf{n}. \quad (15) $$

BC (15) can be strongly simplified in some cases (see Box 1.2).

## 1.3 First-order equation when $D_\omega(z) = 0$

The pressure is no longer given by equation (9) when $D_\omega(z) = 0$. The latter defines the equation of turning surfaces, on which the mathematical problem change of type. Onto such surfaces, the pressure obeys another partial differential equation. Assuming that the field unknowns remain non-singular across such surfaces, this equation is found by using vector manipulations of the primitive equations (Proposition 1.5).
**Proposition 1.5.** When $D_\omega(z) = 0$ (i.e. on turning surfaces), the pressure is given by the first-order equation

$$ (2\Omega \cdot \mathbf{1}_z)(2\Omega \cdot \nabla\Phi) + 2i\omega \Omega \cdot (\mathbf{1}_z \times \nabla\Phi) - \omega^2(\mathbf{1}_z \cdot \nabla\Phi) = 0, \quad (16) $$

which is a particular case of equation (3.5) in Friedlander & Siegmann [2].

*Proof.* We apply $2\Omega \times u$ to equation (1a), which gives

$$ \begin{aligned} i\omega(2\Omega \times u) &= -2\Omega \times (2\Omega \times u) - 2\Omega \times \nabla\Phi - g(z)(\zeta/\rho^*)(2\Omega \times \mathbf{1}_z), \\ &= -4[(\Omega \cdot u)\Omega - |\Omega|^2 u] - 2\Omega \times \nabla\Phi - g(z)(\zeta/\rho^*)(2\Omega \times \mathbf{1}_z). \end{aligned} $$

The latter equation can be rearranged to give

$$ (4|\Omega|^2 - \omega^2) u = -i\omega\nabla\Phi + 4(\Omega \cdot u)\Omega + 2\Omega \times \nabla\Phi + \frac{N^2(z)}{i\omega}(\mathbf{1}_z \cdot u)[2\Omega \times \mathbf{1}_z - i\omega\mathbf{1}_z], \quad (17) $$

where we have used equation (1a) to simplify $2\Omega \times u$ and used equation (2a) to simplify the buoyancy term. Thus, it only remains to express $\Omega \cdot u$ and $\mathbf{1}_z \cdot u$ as a function of $\Phi$ only. We take the scalar product of equation (1a) with $i\omega\Omega$, which gives

$$ \omega^2(\Omega \cdot u) = i\omega(\Omega \cdot \nabla\Phi) + N^2(z)(\Omega \cdot \mathbf{1}_z)(\mathbf{1}_z \cdot u). \quad (18) $$

Then, we take the scalar product of equation (1a) with $i\omega\mathbf{1}_z$ to obtain

$$ [N^2(z) - \omega^2](\mathbf{1}_z \cdot u) = i\omega 2\Omega \cdot (\mathbf{1}_z \times u) - i\omega(\mathbf{1}_z \cdot \nabla\Phi). \quad (19) $$

Since we have also from equation (1a)

$$ i\omega(\mathbf{1}_z \times u) = -\mathbf{1}_z \times (2\Omega \times u) - \mathbf{1}_z \times \nabla\Phi = -2[(\mathbf{1}_z \cdot u)\Omega - (\mathbf{1}_z \cdot \Omega)u] - \mathbf{1}_z \times \nabla\Phi, \quad (20) $$

we combine equations (18) to (20) to obtain

$$ -D_{\omega}(z)(\mathbf{1}_z \cdot u) = i\omega(2\Omega \cdot \mathbf{1}_z)(2\Omega \cdot \nabla\Phi) - 2\omega^2 \Omega \cdot (\mathbf{1}_z \times \nabla\Phi) - i\omega^3(\mathbf{1}_z \cdot \nabla\Phi). \quad (21) $$

From equation (21), we see that the pressure must satisfy equation (16) when $D_\omega(z) = 0$. ■

# 2 Microlocal analysis of boundary-value problems

We investigate, using microlocal analysis, under which conditions a quite general boundary-value scalar problem is elliptic in §2.1. Then, we revisit the ellipsoidal model with a constant $N$ in §2.2 to prove that the low-frequency spectrum $\sigma_2$ is essential.

## 2.1 Elliptic boundary-value problem

**Remark.** In this subsection, the Cartesian coordinates are denoted by $(x_1, x_2, x_3)$ to make use of index notations.
We consider a smooth bounded domain $V$ in $\mathbb{R}^d$, and a symmetric elliptic differential operator $P_\omega$ of order 2 with smooth coefficients and of principal symbol $p$. We also assume that the Dirichlet extension of $P_\omega$ is $> 0$. The scalar $\Phi$ is governed by the boundary-value problem given by

$$ P_{\omega}(\Phi) = 0 \text{ in } V, \quad \mathcal{B}(\Phi) := n^* \cdot \nabla \Phi + \mathcal{V}(\Phi) = 0 \text{ on } \partial V, \qquad (22a,b) $$

where $n^*$ is the unit vector normal computed with respect to the metric $g$ induced on $\partial V$ by the principal symbol of $P_\omega$ (see Box 2.1), and $\mathcal{V}$ is a differential operator of degree 1 tangent<fnref n="1" /> to $\partial V$. To solve system (22), we introduce the Dirichlet-to-Neumann operator $\mathcal{N}$ in Proposition 2.1.

**Proposition 2.1.** If $\Psi : \partial V \to \mathbb{C}$ is a field defined on the boundary, then the Dirichlet-to-Neumann operator is defined by $\mathcal{N}(\Psi) := \partial_{n^*} \Phi$ where $\partial_{n^*} = n^* \cdot \nabla$ is the normal derivative with respect to the vector $n^*$, and where $\Phi$ is the solution of the Dirichlet problem in $V$ given by $P_\omega(\Phi) = 0$ in $V$ and $\Phi|_{\partial V} = \Psi$ on $\partial V$.

Note that $\mathcal{N}$ is a self-adjoint operator acting on functions defined in the space $L^2(\partial V, d\sigma)$ for the area $d\sigma$ defined by the metric $g$ restricted on $\partial V$ [5]. To solve problem (22), we first solve the boundary problem $\mathcal{N}(\Psi) + \mathcal{V}(\Psi) = 0$ on $\partial V$, and then consider the Dirichlet problem in $V$ given the boundary field $\Psi$. In practice, this amounts to solving the problem on the boundary $\partial V$ and, then, propagating the boundary solution in the interior $V$ by using the Dirichlet problem. Actually, the nature of the boundary-value problem crucially depends on the nature of the boundary equation $\mathcal{N}(\Psi) + \mathcal{V}(\Psi) = 0$ on $\partial V$. This can be analysed using microlocal analysis. Indeed, $\mathcal{N}$ is a nice pseudo-differential operator of principal symbol is given in Proposition 2.2. The ellipticity of boundary-value problem (22) is then given by Proposition 2.3. Finally, the nature is given by Proposition 2.4 below.

**Proposition 2.2.** The principal symbol of the Dirichlet-to-Neumann operator $\mathcal{N}$, defined for an elliptic operator with a positive principal symbol $p > 0$, is given by $\sqrt{g^*}|_{\partial V}$, where $g^*|_{\partial V}$ is the dual of the restriction on $\partial V$ of the metric $g$ induced by $p$ in $V$.

*Proof.* It is given in the Appendices of Taylor's book [6]. Note that the principal symbol is defined as $-\sqrt{-g}|_{\partial V}$ in [6]. The first minus sign comes from the fact that an inward vector $n^*$ was considered, and the second from the fact that the other convention $p \le 0$ was used. ■

**Proposition 2.3.** Let us denote by $v$ the principal symbol of $\mathcal{V}$. Then, boundary-value problem (22) is elliptic if and only if the principal symbol $\sqrt{g^*}|_{\partial V} + v$ is elliptic on $\partial V$, that is if it does not vanish for non-zero covectors $k \in \mathbb{R}^2$ on the boundary $\partial V$.

*Proof.* See Chapter 9 in [7]. This is also known as the Shapiro-Lopatinskii condition [7]. ■

**Proposition 2.4.** If boundary-value problem (22) is elliptic in the interval $[\omega_1, \omega_2]$, then its spectrum is discrete (or empty) in this interval. Otherwise, the spectrum is essential in $[\omega_1, \omega_2]$.

*Proof.* If the boundary-value pressure problem is elliptic, the operator $P_\omega$ admits a parametrix $\mathcal{L}$ (i.e. an approximate inverse) such that $P_\omega \circ \mathcal{L} = \text{Id} + \mathcal{R}$, where $\mathcal{R}$ is a pseudo-differential operator of degree $-1$. Hence, $\mathcal{R}$ is a compact operator. Moreover, $P_\omega$ and $\mathcal{R}$ are also analytic functions. Then, the result follows from the Fredholm analytic theorem (see Theorem VI.14 in [9]). Otherwise, the values $\omega \in [\omega_1, \omega_2]$ belong to the essential spectrum [10]. ■

<footnote n="1">¹tangent is here understood in the language of differential geometry [4]. For instance, if $\partial V$ is described by $q_1$ = constant for a set of curvilinear coordinates $(q_1, q_2, q_3)$, then $\mathcal{V}$ only involves derivatives $(\partial_{q_2}, \partial_{q_3})$ on the tangent plane at a given point.</footnote>
BOX 2.1: Normal vector induced by the metric on $\partial V$

The normal vector $\hat{n}$ defined from the metric on $\partial V$ is computed as follows. The principal symbol $\mathfrak{p}$ defines a metric given by $\mathfrak{g}^* := (\mathfrak{g}^*)^{ij}k_ik_j$ for the covectors $k \in \mathbb{R}^3$ in the cotangent space, whose inverse $\mathfrak{g} := (\mathfrak{g}^*)^{-1}$ induces a metric in $V$ defined by $\mathfrak{g} := \mathfrak{g}_{ij}dx_i dx_j$. Then, the normal vector $\mathfrak{n}^*$ computed with respect to $\mathfrak{g}^*$ is given by

$$ \mathfrak{n}^* \propto \mathfrak{g}^* \mathfrak{n}, \quad \mathfrak{g}_{ij} n_i^* n_j^* = 1, $$

where $\mathfrak{n} = (n_1, n_2, n_3)^\top$ is the unit vector normal to $\partial V$ in Cartesian coordinates with $n_1^2 + n_2^2 + n_3^2 = 1$. The constant of proportionality is such that $\mathfrak{n}^*$ has a unit norm with respect to the metric $\mathfrak{g}$. The vector $\mathfrak{n}^*$ is also called the **conormal vector** (with respect to $\mathfrak{g}^*$) in the study of partial differential equations [8]. Note that $\mathfrak{g}^*$ is a contravariant metric tensor associated with the operator $(\mathfrak{g}^*)^{ij}\partial_{x_i}\partial_{x_j}$. Thus, the conormal $\mathfrak{n}^*$ is the contravariant description (in the space with the metric defined by $\mathfrak{g}$) of the covariant normal vector $\mathfrak{n}$ to $\partial V$ (in the space with the Euclidean metric).

In practice, the difficult step is to evaluate the principal symbol of the Dirichlet-to-Neumann operator (Proposition 2.2). In practice, this can be done as follows. The principal symbol $\mathfrak{p}$ induces a metric in $V$ given by $\mathfrak{g} := \mathfrak{g}_{ij}dx_i dx_j$. Let us assume that $x_i = x_i(q_j)$ with the curvilinear coordinates $(q_1, q_2, q_3)$ such that the boundary $\partial V$ is given by $q_1 = \text{constant}$ (for instance). Then, we can rewrite the two-dimensional metric restricted on $\partial V$ as $\mathfrak{g}|_{\partial V} = \tilde{\mathfrak{g}}_{ij}dq_i dq_j$. Then, the dual of the metric induced by $\mathfrak{g}$ on $\partial V$ is given by $\mathfrak{g}^*|_{\partial V} = (\tilde{\mathfrak{g}}^*)^{ij}k_i k_j$ with $\tilde{\mathfrak{g}}^* = \tilde{\mathfrak{g}}^{-1}$.

## 2.2 Back to the pressure equation

We aim to apply the above microlocal analysis to the boundary-value pressure problem when $N$ is constant, which is elliptic in $V$ when $0 < |\omega| < \omega_-$ . The principal symbol $\mathfrak{p}$, defined in equation (12), induces the metric $\mathfrak{g}$ in $V$ and its dual $\mathfrak{g}^*$ characterised by the metric tensors

$$ \mathfrak{g}^* = \begin{pmatrix} A & F & E \\ F & B & D \\ E & D & C \end{pmatrix}, \quad \mathfrak{g} = \frac{1}{\det(\mathfrak{g}^*)} \begin{pmatrix} BC - D^2 & -CF + DE & -BE + DF \\ -CF + DE & AC - E^2 & -AD + EF \\ -BE + DF & -AD + EF & AB - F^2 \end{pmatrix}, \quad (23a,b) $$

with $\det(\mathfrak{g}^*) = ABC - AD^2 - BE^2 - CF^2 + 2DEF$. For a given normal vector $\mathfrak{n} = (n_1, n_2, n_3)^\top$ in Cartesian coordinates, the conormal vector $\mathfrak{n}^*$ is thus given by

$$ \mathfrak{n}^* = \frac{1}{\alpha} \left[ n_1 \begin{pmatrix} A \\ F \\ E \end{pmatrix} + n_2 \begin{pmatrix} F \\ B \\ D \end{pmatrix} + n_3 \begin{pmatrix} E \\ D \\ C \end{pmatrix} \right], \quad (24) $$

with $\alpha = [An_1^2 + Bn_2^2 + Cn_3^2 + 2(Dn_2n_3 + En_1n_3 + Fn_1n_2)]^{1/2}$ such that $\mathfrak{g}_{ij}n_i^* n_j^* = 1$. Thanks to expression (24), we clearly see that the conormal naturally appears in equations (13a-c). Thus, we can apply Proposition 2.3 to determine the nature of the boundary-value pressure problem.

Let us consider here the aligned case with $\Omega = (0, 0, \Omega_z)^\top$, for which the equations are then greatly simplified (see Box 1.2). In this case, the two metrics reduce to

$$ \mathfrak{g}^* = k_x^2 + k_y^2 + Ck_z^2, \quad \mathfrak{g} = dx^2 + dy^2 + C^{-1}dz^2, \quad (25a,b) $$
and the conormal vector is given by $\mathbf{n}^* = (1/\alpha) (n_1, n_2, n_3 C)^\top$ with $\alpha = \sqrt{1 + n_3^2(C - 1)}$, where the Euclidian normal vector is $\mathbf{n} = (n_1, n_2, n_3)^\top$ with $n_1^2 + n_2^2 + n_3^2 = 1$. To demonstrate that the pressure problem is not elliptic for a given frequency $0 < |\omega| < \omega_-$, it is sufficient to find one location on $\partial V$ where Proposition 2.3 is not valid. We look at the behaviour at the equator, where the tangent plane is vertical with $\mathbf{n} = \mathbf{1}_x$ (for instance) such that the pressure BC is simply $f_x = 0$. The restriction of the metrics on $\partial V$ gives

$$ \mathfrak{g}|_{\partial V} = dy^2 + C^{-1}dz^2, \quad \mathfrak{g}^*|_{\partial V} = k_y^2 + Ck_z^2, \qquad (26a,b) $$

and the conormal vector is simply given by $\mathbf{n}^* = \mathbf{1}_x$. The principal symbol of the pressure BC is then given by

$$ \sqrt{k_y^2 + Ck_z^2} = -\frac{2\Omega_z}{\omega}k_y, \qquad (27) $$

which reduces to

$$ \underbrace{\left( \frac{4\Omega_z^2 - \omega^2}{\omega^2} \right)}_{\ge 0} k_y^2 = \underbrace{C}_{\ge 0} k_z^2 \qquad (28) $$

where the two prefactors are positive when $0 < |\omega| < \omega_-$ = $\min(N, 2\Omega_z)$ in the aligned case. Therefore, we can find a real-valued covector $\mathbf{k} = (0, k_y, k_z)^\top$ such that the principal symbol of the BC vanishes. The boundary-value pressure problem is thus not elliptic in this case (Proposition 2.2), and the low-frequency spectrum is essential (Proposition 2.4). Similarly, we conjecture that the spectrum is always essential when $0 < |\omega| < \omega_-$ when $\Omega$ and $\mathbf{1}_z$ are misaligned. The corresponding calculations (which are more lengthy) will be included in a forthcoming paper [5], in which we will also further explore the properties of the low-frequency (surface) modes.

## 3 Boundary-layer analysis

We introduce diffusion in the problem to investigate the diffusive decay rates of IGMs for a constant BV frequency in an ellipsoid. We use the same model as in the main text except that, here, the fluid has a non-zero kinematic viscosity $\nu$ and diffusivity $\kappa$. We employ boundary-layer theory (BLT) to approach the low-diffusive regime (relevant for geophysical vortices), and introduce dimensionless variables below. We take $\Omega_s^{-1}$ as the time scale, $L$ as the length scale (e.g. either $a$ or $c$), and $\rho_*LN^2/g$ as the density scale. For the sake of concision, the dimensionless variables are written below using the same symbols as the dimensional variables in the main text. We seek modal solutions as $[\mathbf{v}, \pi, \rho](\mathbf{r}, t) = [\mathbf{u}, \Phi, \zeta](\mathbf{r}) \exp(\lambda t)$, where $\lambda \in \mathbb{C}$ is the diffusive eigenvalue. The dimensionless linearised Navier-Stokes equations are then given in the rotating frame by

$$ \lambda \mathbf{u} + 2(\mathbf{1}_{\Omega} \times \mathbf{u}) = -\nabla\Phi - \tilde{N}^2\zeta\mathbf{1}_z + E\nabla^2\mathbf{u}, \quad \lambda\zeta = (\mathbf{u} \cdot \mathbf{1}_z) + E_{\kappa}\nabla^2\zeta, \qquad (29a,b) $$

with the (dimensionless) Ekman number $E = \nu/(\Omega_s L^2)$ and $E_\kappa = E/Pr$ where $Pr = \nu/\kappa$ is the Prandtl number, and with the normalised BV frequency $\tilde{N} = N/\Omega_s$. Multiplying equation (29a) by $\mathbf{u}^\dagger$ and the complex-conjugate of equation (29b) by $\tilde{N}^2\zeta$, we obtain after volume integration

$$ \lambda ||\mathbf{u}||^2 + \lambda^\dagger \tilde{N}^2 ||\zeta||^2 = E \langle \mathbf{u}, \nabla^2 \mathbf{u} \rangle + \tilde{N}^2 E_\kappa \langle \nabla^2 \zeta, \zeta \rangle, \qquad (30) $$

where the scalar product between two scalar quantities is defined as $\langle \zeta_1, \zeta_2 \rangle = \int_V \zeta_1^\dagger \cdot \zeta_2 dV$. Note that formula (30) is also valid in a sphere with a radial gravity [11]. Finally, the problem is supplemented with BC. We assume that the velocity field satisfies the stress-free BCs

$$ \mathbf{u} \cdot \mathbf{n}|_{\partial V} = 0, \quad [\epsilon(\mathbf{u}) \cdot \mathbf{n}] \times \mathbf{n}|_{\partial V} = \mathbf{0} \qquad (31a,b) $$
where $\epsilon(\mathbf{u}) = (1/2) [\nabla\mathbf{u} + (\nabla\mathbf{u})^\top]$ is the strain-rate tensor. SF-BCs (31a,b) allow the tangential velocity to slip on $\partial V$. For the density perturbation, we enforce a Neumann BC $\nabla\zeta \cdot \mathbf{n}|_{\partial V} = 0$. A density jump is thus possible between the vortex and the ambient fluid (as considered in idealised models [12, 13, 14]), but without exchange of mass between the two fluids.

We can simplify equation (31) using BLT when $E \ll 1$. Formally, classical BLT for rotating fluids is not expected to be strongly modified by stratification when $\tilde{N}$ is not too large (e.g. for rigid boundaries [15]). Moreover, it is known that the boundary-layer flow is $E^{1/2}$ smaller than the bulk flow for SF-BCs [16, 17]. This considerably simplifies the BLT, because explicit expressions for the boundary-layer solutions are no longer required to estimate the decay rate of the modes. We seek the variables at the leading order in $E$ for the SF-BCs as

$$ \lambda \approx i\omega_0 + E\lambda_1, \quad [\mathbf{u}, \zeta] \approx [\mathbf{u}_0, \zeta_0] + E^{1/2} [\mathbf{u}_1, \zeta_1], \qquad (32a,b) $$

where $[\omega_0, \mathbf{u}_0]$ is the eigenvalue-eigenvector pair of a diffusionless IGM, $\lambda_1 \in \mathbb{C}$ is the first-order correction of the eigenvalue, and $[\mathbf{u}_1, \zeta_1]$ are the first-order corrections within the boundary layer such that $\mathbf{u}_0 + E^{1/2}\mathbf{u}_1$ satisfies SF-BCs (31a,b) and $\zeta_0 + E^{1/2}\zeta_1$ satisfies the Neumann BC for the density. We substitute the above asymptotic expansions into equation (30) and we obtain

$$ \lambda_1 ||\mathbf{u}_0||^2 + \lambda_1^\dagger \tilde{N}^2 ||\zeta_0||^2 \approx - \left( 2 \int_V \epsilon(\mathbf{u}_0) : \epsilon(\mathbf{u}_0^\dagger) dV + \frac{\tilde{N}^2}{Pr} \int_V (\nabla\zeta_0)^2 dV \right) \qquad (33) $$

at the order $E$, where we have used the SF-BCs and the Neumann BC to simplify the volume integrals on the right-hand side. Since the latter is real-valued and negative, there is no frequency correction at the leading order in $E$ due to SF-BCs such that $\Im_m(\lambda_1) = 0$ and $\Re_e(\lambda_1) = \tau_1 \le 0$. We recover from equation (33) the viscous decay rate of pure inertial modes with SF-BCs [17].

Finally, we can crudely estimate whether forced IGMs would quickly decay or not in large-scale stratified vortices by plugging geophysical estimates into equation (33). Typical values for the kinematic viscosity are $\nu \sim 10^{-6} \text{ m}^2\cdot\text{s}^{-1}$ for water and $\nu \sim 4 \times 10^{-7} \text{ m}^2\cdot\text{s}^{-1}$ for gas giants (according to ab-initio simulations [18]). Using Table 1 in the main text, typical values for the Ekman number are thus $E \sim 10^{-8} - 10^{-12}$ for Mediterranean eddies and $E \sim 10^{-12} - 10^{-17}$ for Jovian vortices. Estimating the Prandtl number depends on whether the stratification is due to thermal effects or compositional ones (for which $Pr \gg 1$ in both water and gas giants). Heat diffusion in Mediterranean eddies is characterised by $Pr \approx 0.7$, and typical values $Pr \sim 10^{-2} - 1$ are expected for thermal diffusion in gas giants [18]. Hence, we have $\tilde{N}E_\kappa \sim 10^{-4} - 10^{-8}$ for Mediterranean eddies and $\tilde{N}E_\kappa \sim 10^{-6} - 10^{-11}$ for Jovian vortices. Consequently, the largest-scale IGMs are not expected to be strongly damped by diffusion in stratified vortices.

## References

[1] Greenspan HP. 1968 *The theory of rotating fluids*. Cambridge (UK): Cambridge University Press.

[2] Friedlander S, Siegmann WL. 1982a Internal waves in a rotating stratified fluid in an arbitrary gravitational field. *Geophys. Astrophys. Fluid Dyn.* **19**, 267–291. (doi:10.1080/03091928208208959).

[3] Friedlander S, Siegmann WL. 1982b Internal waves in a contained rotating stratified fluid. *J. Fluid Mech.* **114**, 123–156. (doi:10.1017/S002211208200007X).
[4] Schutz BF. 1980 *Geometrical methods of Mathematical Physics*. Cambridge (UK): Cambridge University Press.

[5] Colin de Verdière Y, Vidal J. 2024 On gravito-inertial surface waves. Preprint.

[6] Taylor M. 2013 *Partial differential equations II: Qualitative studies of linear equations*. New York (USA): Springer.

[7] Grubb G. 2008 *Distributions and operators*. Graduate Texts in Mathematics. Springer, New York (USA): Springer.

[8] Miranda C. 1970 *Partial differential equations of elliptic type*. Springer-Verlag.

[9] Reed M, Simon B. 1972 *Methods of modern Mathematical Physics. I. Functional analysis*. San Diego (USA): Academic Press.

[10] Colin de Verdière Y. 2020 Spectral theory of pseudodifferential operators of degree 0 and an application to forced linear waves. *Anal. PDE* **13**, 1521–1537.

[11] Dintrans B, Rieutord M, Valdettaro L. 1999 Gravito-inertial waves in a rotating stratified sphere or spherical shell. *J. Fluid Mech.* **398**, 271–297. (doi:10.1017/S0022112099006308).

[12] Hassanzadeh P, Marcus PS, Le Gal P. 2012 The universal aspect ratio of vortices in rotating stratified flows: theory and simulation. *J. Fluid Mech.* **706**, 46–57. (doi:10.1017/jfm.2012.180).

[13] Aubert O, Le Bars M, Le Gal P, Marcus PS. 2012 The universal aspect ratio of vortices in rotating stratified flows: experiments and observations. *J. Fluid Mech.* **706**, 34–45. (doi:10.1017/jfm.2012.176).

[14] Lemasquerier D, Facchini G, Favier B, Le Bars M. 2020 Remote determination of the shape of Jupiter's vortices from laboratory experiments. *Nat. Phys.* **16**, 695–700. (doi:10.1038/s41567-020-0833-9).

[15] Friedlander S. 1989 Asymptotic behaviour of decay rates of internal waves in a rotating stratified spherical shell. *Geophys. J. Int.* **96**, 245–252. (doi:10.1111/j.1365-246X.1989.tb04448.x).

[16] Rieutord M, Zahn JP. 1997 Ekman pumping and tidal dissipation in close binaries: a refutation of Tassoul's mechanism. *Astrophys. J.* **474**, 760. (doi:10.1086/303494).

[17] Vidal J, Cébron D. 2023 Precession-driven flows in stress-free ellipsoids. *J. Fluid Mech.* **954**, A5. (doi:10.1017/jfm.2022.976).

[18] French M, Becker A, Lorenzen W, Nettelmann N, Bethkenhagen M, Wicht J, Redmer R. 2012 Ab initio simulations for material properties along the Jupiter adiabat. *Astrophys. J. Suppl.* *Ser.* **202**, 5. (doi:10.1088/0067-0049/202/1/5).