2024

geophysics, fluid mechanics, mathematical physics

Jérémie Vidal

# Inertia-gravity waves in geophysical vortices

Jérémie Vidal1 & Yves Colin de Verdière2 1 Université Grenoble Alpes, CNRS,
ISTerre, 38000 Grenoble, France
2 Université Grenoble Alpes, CNRS, Institut Fourier, 38000 Grenoble, France
ID JV, 0000-0002-3654-6633<EMAIL_ADDRESS>
###### Abstract

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).

###### keywords:

waves, stratification, rotating flows, triaxial ellipsoids

††journal: rspa

##  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].

| \begin{overpic}[height=117.07924pt]{Fig1a.pdf} \put(-10.0,0.0){(a)}
\end{overpic}
---
| \begin{overpic}[height=117.07924pt]{Fig1b.pdf} \put(-20.0,0.0){(b)}
\end{overpic}
---
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 $\boldsymbol{\Omega}$.

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, 10, 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, 21, 22, 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
$\boldsymbol{\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\boldsymbol{\Omega}\cdot\nabla)^{2}\Phi,\quad\left.\nabla\Phi\boldsymbol{\cdot}\boldsymbol{n}\right|_{\partial
V}=\left.\boldsymbol{n}\boldsymbol{\cdot}(\boldsymbol{u}\wedge
2\boldsymbol{\Omega})\right|_{\partial V},$ (1a–b)

where $\boldsymbol{u}$ is the incompressible fluid velocity, $\Phi$ is the
(reduced) hydrodynamic pressure and $\boldsymbol{n}$ is the unit vector normal
to the fluid boundary $\partial V$. Time dependence was assumed to be
$\exp(\mathrm{i}\omega t)$, where $|\omega|<2|\boldsymbol{\Omega}|$ is the
real-valued angular frequency [2]. Equation (1a) is hyperbolic, but boundary
condition (1b) is neither of Dirichlet-type or Neumann-type (it is a mixed-
typed condition, sometimes called oblique condition). Hence, problem (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 $\mathrm{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, 28, 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

###  2.1 Primitive fluid dynamics equations

| \begin{overpic}[height=108.405pt]{./Fig2a.pdf} \put(-20.0,0.0){(a)}
\end{overpic}
---
| \begin{overpic}[height=108.405pt]{./Fig2b.jpg} \put(-15.0,0.0){(b)}
\end{overpic}
---
Figure 2: Sketch of the mathematical model for geophysical vortices (not to
scale). (a) Meridional cross-section of a vortex centred at colatitude
$\theta$, 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 $\boldsymbol{g}$ (i.e.
the background density $\rho_{0}$ decreases with increasing $z$).

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 $\boldsymbol{\Omega}=\Omega_{s}\boldsymbol{1}_{\Omega}$ with respect
to an inertial frame (where $\boldsymbol{1}_{\Omega}$ is a unit vector). We
work in the frame rotating at angular velocity $\boldsymbol{\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 $\boldsymbol{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 $\boldsymbol{U}_{0}$. Laboratory experiments [10, 11] and
numerical calculations [9] show that, for incompressible fluids,
$\boldsymbol{U}_{0}$ is often close to a uniform-vorticity flow. The latter
can be sought as
$\boldsymbol{U}_{0}=\boldsymbol{\omega}\times\boldsymbol{r}+\nabla\Psi$, where
$\boldsymbol{\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
$\boldsymbol{U}_{0}$ is characterised by the Rossby number
$Ro=|\boldsymbol{\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 $\boldsymbol{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
$\boldsymbol{g}=-g(z)\boldsymbol{1}_{z}$ within $V$, where
$\boldsymbol{1}_{z}$ is the unit vector along the vertical $z-$axis. We denote
by $\rho_{*}+\rho_{0}(z)$ the reference density within $V$ (with the mean
density $\rho_{*}$), which is given by the hydrostatic equilibrium $\nabla
P_{0}=(\rho_{*}+\rho_{0})\boldsymbol{g}$ where $P_{0}$ is the hydrostatic
pressure.

Now, we seek small-amplitude perturbations around the reference state
$[\boldsymbol{U}_{0},\rho_{0},P_{0}]$. Without diffusion, the velocity
$\boldsymbol{v}(\boldsymbol{r},t)$ obeys the linearised Euler equation given
in the rotating frame by

$\partial_{t}\boldsymbol{v}+(\boldsymbol{U}_{0}\boldsymbol{\cdot}\nabla)\boldsymbol{v}+(\boldsymbol{v}\boldsymbol{\cdot}\nabla)\boldsymbol{U}_{0}+2\boldsymbol{\Omega}\times\boldsymbol{v}=-\nabla\pi+({\rho}/{\rho_{*}})\,\boldsymbol{g}$
(2)

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

$\partial_{t}\rho/\rho_{*}+\boldsymbol{U}_{0}\cdot\nabla\rho/\rho_{*}=(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{1}_{z})\,N^{2}(z)/g(z),$
(3)

where $N^{2}(z)=\boldsymbol{g}\boldsymbol{\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 $\boldsymbol{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 $\boldsymbol{U}_{0}=\boldsymbol{0}$ in
equations (2)-(3) to focus on the free wave motions given by

$\partial_{t}\boldsymbol{v}+2\boldsymbol{\Omega}\times\boldsymbol{v}=-\nabla\pi-N(z)\xi\,\boldsymbol{1}_{z},\quad\nabla\boldsymbol{\cdot}\boldsymbol{v}=0,\quad\partial_{t}\xi=(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{1}_{z})\,N(z),$
(4a–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 (4a-c) are supplemented with boundary conditions (BCs). We
impose that the fields are regular at $\boldsymbol{r}=\boldsymbol{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
$\left.\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{n}\right|_{\partial V}=0$
on $\partial V$, where $\boldsymbol{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 (4a) on
$\partial V$ (see the ESM).

###  2.2 Mixed hyperbolic-elliptic problem

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

$\left[\boldsymbol{v},\pi,\xi\right](\boldsymbol{r},t)=\left[\boldsymbol{u},\Phi,\zeta\right](\boldsymbol{r})\exp(\mathrm{i}\omega
t),$ (5)

where $\left[\boldsymbol{u},\Phi,\zeta\right]$ 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 $\boldsymbol{u}$ or $\Phi$. If we multiply Euler
equation (4a) by $\mathrm{i}\omega$ and replace $\mathrm{i}\omega\zeta$ by
buoyancy equation (4c), we obtain for $\omega\neq 0$

$-\omega^{2}\boldsymbol{u}+2\boldsymbol{\Omega}\times(\mathrm{i}\omega\boldsymbol{u})+N^{2}(z)(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{1}_{z})\boldsymbol{1}_{z}=-\mathrm{i}\omega\nabla\Phi,\quad\nabla\boldsymbol{\cdot}\boldsymbol{u}=0,$
(6a,b)

Equations (6a,b) and the no-penetration BC
$\left.\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n}\right|_{\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\boldsymbol{\Omega}\times\boldsymbol{u}=-\nabla\Phi$. When $N(z)>0$,
the steady modes are two-dimensional (2D) such that
$\boldsymbol{u}\boldsymbol{\cdot}~{}\boldsymbol{1}_{z}=0$, and the density
perturbation obeys the diagnostic equation
$\rho_{*}\left(2\boldsymbol{\Omega}\boldsymbol{\cdot}\nabla\boldsymbol{u}\right)=~{}\nabla\zeta\times\boldsymbol{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$ [39, 32]. 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)\left(\boldsymbol{1}_{z}\boldsymbol{\cdot}\nabla\right)^{2}\Phi-\left(2\boldsymbol{\Omega}\boldsymbol{\cdot}\nabla\right)^{2}\Phi+\dots=0$
(7)

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

$2\omega_{\pm}^{2}(z)=[N^{2}(z)+4|\boldsymbol{\Omega}|^{2}]\pm\sqrt{\left[N^{2}(z)+4|\boldsymbol{\Omega}|^{2}\right]^{2}-16N^{2}(z)(\boldsymbol{\Omega}\boldsymbol{\cdot}\boldsymbol{1}_{z})^{2}}.$
(8)

The lower-order terms in equation (7), 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 (7) 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~{}\mathrm{i}\boldsymbol{k}$, where
$\boldsymbol{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 $\boldsymbol{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

$|\boldsymbol{k}|^{2}\omega^{2}-\left[N^{2}(z)|\boldsymbol{1}_{z}\times\boldsymbol{k}|^{2}+(2\boldsymbol{\Omega}\boldsymbol{\cdot}\boldsymbol{k})^{2}\right]=0,$
(9)

which is the dispersion relation of IGWs in unbounded fluids when $N$ is
constant [3]. The wave-like equation is thus hyperbolic when equation (9) is
satisfied for some non-zero wave vectors $\boldsymbol{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 (9). The latter defines a
quadric cone given by $(\boldsymbol{\mathfrak{g}}^{\star})^{ij}k_{i}k_{j}=0$
(using Einstein notation) when $\omega\neq 0$, where
$\boldsymbol{\mathfrak{g}}^{\star}$ is the associated metric tensor. The
nature of this quadric cone depends on the eigenvalues of
$\boldsymbol{\mathfrak{g}}^{\star}$, 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)$.

###  2.3 Spectral problem for the velocity

Further wave properties can be obtained by considering the velocity equation.
We introduce the Hilbert space $\boldsymbol{\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\boldsymbol{u},\boldsymbol{v}\rangle:=\int_{V}\boldsymbol{u}^{\dagger}\boldsymbol{\cdot}\boldsymbol{v}\,\mathrm{d}V,$
(10)

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

$\boldsymbol{\mathcal{V}}^{0}:=\\{\boldsymbol{u}\in\boldsymbol{\mathcal{V}}\
|\ \nabla\boldsymbol{\cdot}\boldsymbol{u}=0\ \text{in V},\
\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{n}=0\ \text{on $\partial V$}\\}.$
(11)

Finally, we introduce the orthogonal projector $\mathbb{L}$ from
$\boldsymbol{\mathcal{V}}$ onto $\boldsymbol{\mathcal{V}}^{0}$ (called Leray
projector after [52]). For any vector
$\boldsymbol{e}\in\boldsymbol{\mathcal{V}}$ written as
$\boldsymbol{e}=\nabla\Psi+\nabla\times\boldsymbol{A}$ (using a Helmholtz
decomposition) for some potentials $[\Psi,\boldsymbol{A}]$, the projected
vector is defined by $\mathbb{L}(\boldsymbol{e}):=\boldsymbol{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 (6)
with $\boldsymbol{u}\in\boldsymbol{\mathcal{V}}^{0}$ in ellipsoids. We apply
the orthogonal projector to the equation to remove the pressure term, and find
that $(\omega,\boldsymbol{u})$ obey a quadratic eigenvalue problem (QEP). The
latter is given by

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

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

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

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

######  Proposition 2.1.

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

###### Proof.

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

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

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

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

###### Proof.

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

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

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

The spectrum $\sigma(\boldsymbol{\mathcal{Q}}_{\omega})$ is thus the disjoint
union of the point spectrum given by

$\sigma_{P}(\boldsymbol{\mathcal{Q}}_{\omega}):=\\{\omega\in\mathbb{C}\ |\
\exists\,\boldsymbol{u}\neq\boldsymbol{0}\
\,\boldsymbol{\mathcal{Q}}_{\omega}(\boldsymbol{u})=\boldsymbol{0}\\},$ (15)

which is the subset spanned by the eigenvalues of QEP (12), and the remaining
continuous spectrum (i.e. the set of $\omega$ for which
$\boldsymbol{\mathcal{T}}-\omega\boldsymbol{\mathcal{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 (12) is real-valued and even. Moreover, the spectrum is
bounded by $|\omega|\leq\omega_{\max}$ with
$\omega_{\max}\leq\Omega_{s}+\sqrt{\Omega_{s}^{2}+\max(N^{2})}$. Finally, two
eigen-pairs $(\omega_{1},\boldsymbol{u}_{1})$ and
$(\omega_{2},\boldsymbol{u}_{2})$ with $\omega_{1}\neq\omega_{2}$ are not
orthogonal with respect to inner product (10), but satisfy the two integral
properties

$\omega_{1}\omega_{2}\langle\boldsymbol{u}_{2},\boldsymbol{u}_{1}\rangle=-\langle\boldsymbol{u}_{2},\boldsymbol{\mathcal{K}}(\boldsymbol{u}_{1})\rangle,\quad(\omega_{2}+\omega_{1})\langle\boldsymbol{u}_{2},\boldsymbol{u}_{1}\rangle=\langle\boldsymbol{u}_{2},\mathrm{i}\boldsymbol{\mathcal{C}}(\boldsymbol{u}_{1})\rangle.$
(16a,b)

###### Proof.

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

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

which shows that $(-\omega,\boldsymbol{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 $\boldsymbol{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
$\boldsymbol{g}=-g\boldsymbol{1}_{z}$ 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,$ (17)

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

###  3.1 Polynomial description

Since an ellipsoid is a quadratic surface and expression (17) is a polynomial,
we can seek eigenvectors using polynomial functions in the Cartesian
coordinates. We denote by $\mathcal{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\geq 1$. Its dimension is
$N(\mathcal{P}_{n})=(n+1)(n+2)(n+3)/6$. Then, we denote by
$\boldsymbol{\mathcal{V}}_{n}$ the space of vectors in ellipsoids with
coefficients in $\mathcal{P}_{n}$ and define the subspace
$\boldsymbol{\mathcal{V}}_{n}^{0}:=\boldsymbol{\mathcal{V}}_{n}\cap\boldsymbol{\mathcal{V}}^{0}$
of dimension $N(\boldsymbol{\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
$\boldsymbol{\mathcal{W}}_{n}^{0}:=~{}(\boldsymbol{\mathcal{V}}_{n-1}^{0})^{\perp}\cap\boldsymbol{\mathcal{V}}_{n}^{0}$
of dimension $N(\boldsymbol{\mathcal{W}}_{n}^{0})=n(n+2)$, which are spanned
by vectors in $\boldsymbol{\mathcal{V}}_{n}^{0}$ of degree $n$ that are
orthogonal to polynomial vectors of degree $\leq n-1$. These polynomial spaces
are key to characterise the spectrum because of Proposition 3.1.

######  Proposition 3.1.

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

###### Proof.

The invariance of $\boldsymbol{\mathcal{V}}^{0}_{n}$ under the action of
$\boldsymbol{\mathcal{C}}$ has been demonstrated elsewhere [29, 40]. Here, we
can use similar arguments to show that it is invariant by the action of
$\boldsymbol{\mathcal{K}}$. Briefly, we must prove that the Leray projector
$\mathbb{L}$ is well-defined on $\boldsymbol{\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 $\boldsymbol{\mathcal{V}}^{0}_{n}$ is
invariant by the action of $\boldsymbol{\mathcal{Q}}_{\omega}$. Moreover,
since $\oplus_{n\in\mathbb{N}^{*}}\boldsymbol{\mathcal{V}}_{n}^{0}$ is dense
in $\boldsymbol{\mathcal{V}}^{0}$ [40], it shows that
$\sigma(\boldsymbol{\mathcal{Q}}_{\omega})$ is pure point in ellipsoids (i.e.
$\sigma=\sigma_{P}$). The continuous spectrum is thus empty, and equation (6)
admits square-integrable eigenvectors in $\boldsymbol{\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[62, 29], 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 $\boldsymbol{u}$ to
$\boldsymbol{\mathcal{V}}_{n}^{0}$ in ellipsoids.

###  3.2 Galerkin algorithm

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

$\left(-\omega^{2}\boldsymbol{M}_{n}+\omega\mathrm{i}\boldsymbol{C}_{n}+\boldsymbol{K}_{n}\right)\boldsymbol{\alpha}=\boldsymbol{0}$
(18)

where $\boldsymbol{\alpha}=(\alpha_{1},\alpha_{2},\dots)^{\top}$ is the
eigenvector, and where
$[\boldsymbol{M}_{n},\boldsymbol{C}_{n},\boldsymbol{K}_{n}]$ are three real-
valued matrices of size $N({\boldsymbol{\mathcal{V}}_{n}^{0}})\times
N({\boldsymbol{\mathcal{V}}_{n}^{0}})$ and of elements
$M_{n,ij}=\langle\boldsymbol{e}_{i},\boldsymbol{e}_{j}\rangle$,
$C_{n,ij}=\langle\boldsymbol{e}_{i},2\boldsymbol{\Omega}\times\boldsymbol{e}_{j}\rangle$
and
$K_{n,ij}=~{}\langle\boldsymbol{e}_{i},N^{2}(\boldsymbol{e}_{j}\boldsymbol{\cdot}\boldsymbol{1}_{z})\boldsymbol{1}_{z}\rangle$.
The matrix $\boldsymbol{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 $\boldsymbol{C}_{n}$ is antisymmetric,
and $\boldsymbol{K}_{n}$ is a real-valued positive matrix. Further properties
of $[\boldsymbol{C}_{n},\boldsymbol{K}_{n}]$ are given below in Proposition
3.2.

######  Proposition 3.2.

$\dim(\ker\boldsymbol{C}_{n})=\lceil n/2\rceil$ and
$\dim(\ker\boldsymbol{K}_{n})=n(n+1)(n+2)/6$.

###### Proof.

When $N=0$, $\ker(\boldsymbol{C}_{n})$ is associated with geostrophic modes of
degree $n$. Its dimension was proved in [40]. When $N\neq 0$,
$\ker(\boldsymbol{K}_{n})$ is associated with 2D motions such that
$\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{1}_{z}=0$. The latter can be
sought as $\boldsymbol{u}=\nabla\times(\Psi\boldsymbol{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(\boldsymbol{K}_{n})$ is thus
$\dim(\mathcal{P}_{n-1})=n(n+1)(n+2)/6$. ∎

|
---|---
Figure 3: Sparse upper triangular block matrices $\boldsymbol{A}_{n}$ (_left_)
and $\boldsymbol{B}_{n}$ (_right_) for $n=5$ and
$\boldsymbol{1}_{\Omega}=\boldsymbol{1}_{z}$.

For numerical convenience, we convert the finite-dimensional QEP into an
equivalent finite-dimensional SEP when $\omega\neq 0$. The finite-dimensional
counterpart of the self-adjoint operator defined in equation (14b) is not well
suited for numerical computations (because of the square root). Instead, we
consider the SEP $\boldsymbol{L}_{n}\boldsymbol{z}=\omega\boldsymbol{z}$ with
$\boldsymbol{z}=(\boldsymbol{\alpha},\boldsymbol{\beta})^{\top}$ and the
matrix $\boldsymbol{L}_{n}$ given by

$\boldsymbol{L}_{n}:=\begin{pmatrix}0&\boldsymbol{I}\\\
\boldsymbol{B}_{n}&\mathrm{i}\boldsymbol{A}_{n}\\\
\end{pmatrix},\quad\boldsymbol{A}_{n}:=\boldsymbol{M}_{n}^{-1}\boldsymbol{C}_{n},\quad\boldsymbol{B}_{n}:=\boldsymbol{M}_{n}^{-1}\boldsymbol{K}_{n},$
(19a–c)

where $\boldsymbol{I}$ is the identity matrix. As illustrated in figure 3, the
matrices $\boldsymbol{A}_{n}$ and $\boldsymbol{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
$\boldsymbol{A}_{n}\boldsymbol{\alpha}=\omega\boldsymbol{\alpha}$, can indeed
be computed separately for every space $\boldsymbol{\mathcal{W}}_{n}^{0}$ by
seeking the eigenvalues of every diagonal block of $\boldsymbol{A}_{n}$. Yet,
since the matrix $\boldsymbol{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
$[\boldsymbol{M}_{n},\boldsymbol{C}_{n},\boldsymbol{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
$\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{1}^{0}$. The latter are sought in
the form of uniform-vorticity flows as [42]

$\boldsymbol{u}=\boldsymbol{\alpha}\times\boldsymbol{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,$
(20a–b)

and with $\boldsymbol{\alpha}=(\alpha_{x},\alpha_{y},\alpha_{z})^{\top}$.
Next, the block matrices of $\boldsymbol{L}_{1}$ are given by

$\boldsymbol{B}_{1}=\begin{pmatrix}\dfrac{N^{2}c^{2}}{b^{2}+c^{2}}&0&0\\\
0&\dfrac{N^{2}c^{2}}{a^{2}+c^{2}}&0\\\ 0&0&0\\\
\end{pmatrix},\quad\boldsymbol{A}_{1}=\begin{pmatrix}0&-\dfrac{2\Omega_{z}a^{2}}{a^{2}+c^{2}}&\dfrac{2\Omega_{y}a^{2}}{a^{2}+b^{2}}\\\
\dfrac{2\Omega_{z}b^{2}}{b^{2}+c^{2}}&0&-\dfrac{2\Omega_{x}b^{2}}{a^{2}+b^{2}}\\\
-\frac{2\Omega_{y}c^{2}}{b^{2}+c^{2}}&\dfrac{2\Omega_{x}c^{2}}{a^{2}+c^{2}}&0\\\
\end{pmatrix},$ (21a,b)

and with $\boldsymbol{\Omega}=(\Omega_{x},\Omega_{y},\Omega_{z})^{\top}$.
Finally, the eigenvalue-eigenvector solutions of $\boldsymbol{L}_{1}$ can be
solved analytically to extend the linear solutions obtained by Vantieghem [27]
when $N=0$ and $\boldsymbol{1}_{\Omega}=\boldsymbol{1}_{z}$. However, we must
rely on numerical computations to obtain $\boldsymbol{L}_{n}$ for higher
degrees $n$.

###  3.3 Enumeration of the eigenvalues

We can explicitly count the number of non-zero and zero eigenvalues of the QEP
when $\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{n}^{0}$. Since the spectra of
the primitive equations, the one of $\boldsymbol{\mathcal{Q}}_{\omega}$ and
the one of $\boldsymbol{\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 (18) is $n(n+1)(n+5)/3$.

###### Proof.

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

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

Solutions of the latter equations belong to the space spanned by
$\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{n}^{0}$ and
$\zeta\in\mathcal{P}_{n}$, whose dimension is
$N(\boldsymbol{\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
$\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{1}_{z}=0$. Such 2D flows
$\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{n}^{0}$ can be sought as
$\boldsymbol{u}=\nabla\times(\Psi\boldsymbol{1}_{z})$ where
$\Psi=p[1-F(x,y,z)]$ is a stream function and with $p\in\mathcal{P}_{n-1}$.
Assuming that $\boldsymbol{\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\geq 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 $\boldsymbol{L}_{n}$, which has $2N(\boldsymbol{\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\boldsymbol{L}_{n}$. The null space of the
equation $\boldsymbol{L}_{n}(\boldsymbol{\alpha},\boldsymbol{\beta})^{\top}=0$
is spanned by the vector elements such that
$\boldsymbol{\beta}=\boldsymbol{0}$ and
$\boldsymbol{K}_{n}\boldsymbol{\alpha}=~{}0$. Therefore, we have
$\ker\boldsymbol{L}_{n}=\ker\boldsymbol{K}_{n}\bigoplus\\{0\\}$ and $\omega=0$
has the geometric multiplicity $m_{L}=n(n+1)(n+2)/6$ according to Proposition
3.2. However, since $\boldsymbol{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(\boldsymbol{L}_{n}^{2})$.

######  Proposition 3.4.

The eigenvalue $\omega=0$ of the matrix $\boldsymbol{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
$\boldsymbol{L}_{n}\boldsymbol{z}\neq\boldsymbol{0}$ and
$\boldsymbol{L}_{n}^{2}\boldsymbol{z}=\boldsymbol{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
$\boldsymbol{L}_{n}$. Indeed, QEP (18) can be converted into an equivalent
self-ajoint SEP [59], which is the finite-dimensional counterpart of the self-
adjoint operator defined in equation (14b). 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 $\boldsymbol{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

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|\leq\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.

We have shown that $\sigma(\boldsymbol{\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\backslash\sigma_{D}$ (i.e. the complement to the discrete
spectrum). Since $\oplus_{n\in\mathbb{N}^{*}}\boldsymbol{\mathcal{V}}_{n}^{0}$
is dense in $\boldsymbol{\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 (8)
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 $\boldsymbol{\Omega}$ and $\boldsymbol{1}_{z}$ are collinear. We also
show below in Proposition 4.1 that $\omega_{\max}=\omega_{+}$, even when
$\boldsymbol{\Omega}$ and $\boldsymbol{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
$\boldsymbol{z}=~{}(\boldsymbol{u},\zeta)^{\top}$ by
$\mathbb{P}(\boldsymbol{z}):=(\mathbb{L}\boldsymbol{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
$\boldsymbol{\mathcal{P}}$ defined by
$\boldsymbol{\mathcal{P}}(\boldsymbol{z}):=\mathbb{P}(\boldsymbol{\mathcal{E}}\mathbb{P}(\boldsymbol{z}))$
with

$\boldsymbol{\mathcal{E}}=\mathrm{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}.$

$\boldsymbol{\mathcal{P}}$ is a self-adjoint operator for square-integrable
fields. We have $||\mathbb{P}||=1$, and the eigenvalues of the Hermitian
matrix $\boldsymbol{\mathcal{E}}$ are $\pm\omega_{\pm}$ such that
$||\boldsymbol{\mathcal{E}}||\leq\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 $\boldsymbol{u}$ is restricted to
$\boldsymbol{\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 (9). 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.

###  4.1 Low-frequency (surface) modes

####  4.1.1 Illustrative examples

Contrary to the free-space case, it is known that $\sigma_{2}$ is not empty
when $\boldsymbol{\Omega}\propto\boldsymbol{1}_{z}$ in some bounded geometries
(e.g. in cylinders or spheres [39, 21]). In triaxial ellipsoids, we can also
show that $\sigma_{2}\neq\emptyset$ when
$\boldsymbol{\Omega}\boldsymbol{\cdot}\boldsymbol{1}_{z}\neq 0$ (otherwise,
$\omega_{-}=0$ and $\sigma_{1}$ fills the entire interval
$0<|\omega|<\omega_{+}$). To this end, we first consider the eigenvalues for
$\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{1}^{0}$, which are the roots of
the characteristic equation $\det(\boldsymbol{L}_{1}-\omega\boldsymbol{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}},$
(22a–b)

where we have introduced the positive coefficients

$\displaystyle\allowdisplaybreaks\small\beta_{4}$
$\displaystyle=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}),$
$\displaystyle\beta_{2}$
$\displaystyle=c^{4}N^{2}[4(\Omega_{x}^{2}b^{2}+\Omega_{y}^{2}a^{2})+N^{2}(a^{2}+b^{2})],$
$\displaystyle\Delta$
$\displaystyle=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}.$

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 $\boldsymbol{\mathcal{V}}_{1}^{0}$ obtained by
Vantieghem [27] when $\Omega_{x}=\Omega_{y}=0$. Solutions (22) 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.

\begin{overpic}[width=212.47617pt]{./Fig5a.png} \put(-5.0,70.0){(a)} \end{overpic} | \begin{overpic}[width=212.47617pt]{./Fig5b.pdf} \put(-5.0,70.0){(b)} \end{overpic}
---|---
Figure 5: Spectrum $\sigma_{2}$ in the aligned case
$\boldsymbol{1}_{\Omega}=\boldsymbol{1}_{z}$. (a) Eigenvalues
$\omega_{3,5}\geq 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}$.

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}\neq\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\leq
l\leq 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
$\boldsymbol{1}_{\Omega}=\boldsymbol{1}_{z}$ with $N=2\Omega_{s}$. For this
particular configuration, we have $\omega_{-}=\omega_{+}$ from equation (8).
Then, pressure equation (7) and its associated BC for eigenvalues
$\omega\in\sigma_{2}$ reduce to

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

Laplace equation (23a) 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 (23b). In a sphere, the
pressure associated with $\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{n}^{0}$
is given in spherical coordinates $(r,\theta,\phi)$ by the solid harmonics
$\Phi(\boldsymbol{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\leq
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[\mathrm{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
$\boldsymbol{1}_{\Omega}=\boldsymbol{1}_{z}$ and $N=2\Omega_{s}$.

####  4.1.2 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 (7) 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 $\boldsymbol{\Omega}$ and
$\boldsymbol{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$.

Figure 6: Spatial structure of the low-frequency (surface) modes for
$\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{10}^{0}$ in the aligned case
$\boldsymbol{1}_{\Omega}=\boldsymbol{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).

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 $\boldsymbol{\Omega}=\boldsymbol{0}$) or pure IGMs (when
$\boldsymbol{\Omega}\boldsymbol{\cdot}\boldsymbol{1}_{z}=0$) fill the entire
interval $0<|\omega|<\omega_{+}$. Two modal families are shown in figure 6.
When $N/\Omega_{s}\leq 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 when $N/\Omega_{s}\ll 1$. As investigated by Allen [38]
when $\boldsymbol{\Omega}\propto\boldsymbol{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).

###  4.2 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 $\boldsymbol{1}_{\Omega}=\boldsymbol{1}_{z}$.
From dispersion relation (9), the angle $\vartheta$ between
$\boldsymbol{1}_{z}$ and $\boldsymbol{k}$ is such that $\vartheta\to 0$ when
$\omega/\Omega_{s}\to 2$. Since the group velocity is orthogonal to
$\boldsymbol{k}$ for IGWs [3], the energy propagates in the horizontal
direction (in agreement with the observed pancake-like structure). On the
contrary, we have $\theta\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.

Figure 7: Spatial structure of IGMs for
$\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{10}^{0}$ in the aligned case
$\boldsymbol{1}_{\Omega}=\boldsymbol{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).

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
$\boldsymbol{\mathcal{Q}}_{\omega}|_{\boldsymbol{\mathcal{W}}_{n}^{0}}$ given
by $\pi_{n}:=(1/d_{n})\sum_{j=1}^{d_{n}}\delta(\omega_{j})$ with
$d_{n}=2N(\boldsymbol{\mathcal{W}}_{n}^{0})$, where $\delta$ is the Dirac
function. As explained above, the eigenvalues of SEP (19) cannot be easily
separated in different subsets for every degree $n$. Therefore, we rather
consider the joint repartition of the eigenvalues of
$\boldsymbol{\mathcal{Q}}_{\omega}|_{\boldsymbol{\mathcal{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.

|
---|---
Figure 8: Density of IGMs in $\sigma_{1}$. Comparison between asymptotic
density $\mathrm{d}\pi_{\infty}/\mathrm{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
$\boldsymbol{1}_{\Omega}=\boldsymbol{1}_{z}$ with $N/\Omega_{s}=5$. _Right_ :
Sphere $b/a=c/a=1$ in the perpendicular case
$\boldsymbol{1}_{\Omega}\boldsymbol{\cdot}\boldsymbol{1}_{z}=0$ with
$N/\Omega_{s}=10$.

######  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}\mathrm{d}\pi_{\infty}=\frac{1}{8\pi}\mathrm{Area}(\mathfrak{S}_{\lambda}\cap
S^{2})$

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

$\mathfrak{S}_{\lambda}:=\\{\widetilde{\boldsymbol{k}}\in\mathbb{R}^{3}\ |\
\omega(\widetilde{\boldsymbol{k}})\leq\lambda\\},\quad\widetilde{\boldsymbol{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(\widetilde{\boldsymbol{k}})$ is the plane-wave dispersion relation
given by equation (9) for the rescaled wave vector
$\widetilde{\boldsymbol{k}}$. Hence, the asymptotic measure has a probability
density function $f_{\infty}:=~{}\mathrm{d}\pi_{\infty}/\mathrm{d}\omega$ that
is even and non-vanishing when
$\omega_{-}^{2}\leq\omega^{2}\leq\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
$\boldsymbol{\mathcal{T}}$ defined in equation (14b). Thus, the starting point
is to consider the matrix operator $\boldsymbol{\mathcal{L}}$ defined in
equation (14a). 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 $\boldsymbol{\mathcal{L}}$ when $\omega\neq
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\boldsymbol{\mathcal{L}}$ is not clearly related to the steady
solutions of the original physical problem (due to the existence of
generalized eigenvectors). ∎

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}\mathrm{d}\pi_{\infty}=1$ but,
since Proposition 4.2 heavily relies on the dispersion relation of IGWs, it is
assumed that $\mathrm{d}\pi_{\infty}/\mathrm{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 $\boldsymbol{1}_{\Omega}\boldsymbol{\cdot}\boldsymbol{1}_{z}=0$,
such that $\omega_{-}=0$ from expression (8). 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}\leq 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$.

|
---|---
Figure 9: Asymptotic density
$\mathrm{d}\pi_{\infty}/\mathrm{d}\widetilde{\omega}$ of IGMs with normalised
eigenvalues $\widetilde{\omega}:=\omega/\omega_{+}$, given by Proposition 4.2.
_Left_ : Sphere $b/a=c/a=1$ in the aligned case
$\boldsymbol{1}_{\Omega}=\boldsymbol{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$.

Finally, the situation $\omega_{+}=\omega_{-}$ occurring when
$\boldsymbol{1}_{\Omega}=\boldsymbol{1}_{z}$ 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 (23a,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 with $N/\Omega_{s}\gg 1$ and, consequently, they usually have a
very small aspect ratio $c/a\leq 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, 10, 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 (8), which depend on the
orientation of $\boldsymbol{\Omega}$ and on the ratio $N/\Omega_{s}$, can be
simplified in the relevant regime $N/\Omega_{s}\gg 1$. It gives

$\omega_{+}^{2}\simeq
N^{2}+4\left[|\boldsymbol{\Omega}|^{2}-(\boldsymbol{\Omega}\boldsymbol{\cdot}\boldsymbol{1}_{z})^{2}\right],\quad\omega_{-}^{2}\simeq
4|\boldsymbol{\Omega}\boldsymbol{\cdot}\boldsymbol{1}_{z}|^{2}.$ (24a,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 $\boldsymbol{g}=-(g/R)\,\boldsymbol{r}$. The
background hydrostatic density is then
$\nabla\rho_{0}=-(\beta/R)\,\boldsymbol{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},$ (25)

where $s=\sqrt{x^{2}+y^{2}}$ is the cylindrical radius measured from the
rotation axis $\boldsymbol{\Omega}=\Omega_{s}\boldsymbol{1}_{z}$. For a given
angular frequency $\omega$, equation (25) 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.

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\boldsymbol{\Omega}\boldsymbol{\cdot}\boldsymbol{1}_{z}=2\Omega_{z}\cos(\theta)$, where $\theta$ the mean colatitude of the vortex (figure 2). |  $f$ [rad.s-1] | $N/f$ |  $a$ [km] |  $b$ [km] |  $c$ [km]
---|---|---|---|---|---
Jupiter’s GRS | $\sim 1.4\times 10^{-4}$ | $\mathcal{O}(10^{2})$ | $\sim 9000$ | $\sim 5000$ | $\sim 100$
Jupiter’s Oval BA | $\sim 1.9\times 10^{-4}$ | $\mathcal{O}(10^{2})$ | $\sim 3000$ | $\sim 3000$ | $\sim 50$
Mediterranean eddies | $\sim 8.5\times 10^{-5}$ | $10-50$ | $10-10^{2}$ | $10-10^{2}$ | $0.5-1$

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 $(\omega,\boldsymbol{u})$ satisfies in an
ellipsoid the intriguing integral property [27, 28]

$\int_{V}\boldsymbol{u}^{\dagger}\boldsymbol{\cdot}\nabla^{2}\boldsymbol{u}\,\mathrm{d}V=\int_{V}\boldsymbol{u}^{\dagger}\boldsymbol{\cdot}\mathbb{L}(\nabla^{2}\boldsymbol{u})\,\mathrm{d}V=0.$
(26)

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 (26) is valid in an ellipsoid
because (i)
$\mathbb{L}(\nabla^{2}\boldsymbol{u})\in\boldsymbol{\mathcal{V}}_{n-2}^{0}$
when $\boldsymbol{u}\in\boldsymbol{\mathcal{V}}_{n}^{0}$ and (ii) two inertial
eigenvectors $(\boldsymbol{u}_{i},\boldsymbol{u}_{j})$ are orthogonal with
respect to inner product (10) such that
$\langle\boldsymbol{u}_{i},\boldsymbol{u}_{j}\rangle=\delta_{ij}$ (when
properly normalised). Here, the IGMs do not satisfy property (26) because they
are not mutually orthogonal with respect to inner product (10) 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 [75]. 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 $\boldsymbol{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 [76, 21]. Moreover, there is a class of
magnetic fields leaving invariant the three-dimensional polynomial vector
spaces in an ellipsoid [77, 78]. 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 [79, 80, 81, 82]. 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.

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

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.

The authors declare that they have no competing interests.

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

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.

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. Springer, New York (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:SciPostPhys.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] Colin de Verdière Y, Vidal J. 2024 On gravito-inertial surface waves. Preprint. (doi:hal-04461197).
  * [76] Friedlander S. 1987 Hydromagnetic waves in the Earth’s fluid core. Geophys. Astrophys. Fluid Dyn. 39, 315–333. (doi:10.1080/03091928708208816).
  * [77] Malkus WVR. 1967 Hydromagnetic planetary waves. J. Fluid Mech. 28, 793–802. (doi:10.1017/S0022112067002447).
  * [78] 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).
  * [79] 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).
  * [80] 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).
  * [81] 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).
  * [82] 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).