Jérémie Vidal
email: jeremie.vidal@univ-grenoble-alpes.fr

Yves Colin De Verdière

Inertia-gravity waves in geophysical vortices

Keywords: geophysics, fluid mechanics, mathematical physics waves, stratification, rotating flows, triaxial ellipsoids

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

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]. 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 pancakelike 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 [START_REF] Kerswell | Elliptical instability[END_REF]. 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 [START_REF] Miyazaki | Three-dimensional instability of strained vortices in a stably stratified fluid[END_REF][START_REF] Kerswell | Elliptical instabilities of stratified, hydromagnetic waves[END_REF][START_REF] Bars | Thermo-elliptical instability in a rotating cylindrical shell[END_REF][START_REF] Vidal | Fossil field decay due to nonlinear tides in massive binaries[END_REF]. 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 Ω (with respect to an inertial frame), inertial modes in a bounded domain V are governed by the Poincaré problem (called after Cartan [START_REF] Cartan | Sur les petites oscillations d'une masse fluide[END_REF] who revisited Poincaré's pioneering paper [START_REF] Poincaré | Sur l'équilibre d'une masse fluide animée d'un mouvement de rotation[END_REF]). The latter is given by

ω 2 ∇ 2 Φ = (2Ω • ∇) 2 Φ, ∇Φ • n| ∂V = n • (u ∧ 2Ω)| ∂V , (1.1a-b)
where u is the incompressible fluid velocity, Φ is the (reduced) hydrodynamic pressure and n is the unit vector normal to the fluid boundary ∂V . Time dependence was assumed to be exp(iωt), where |ω| < 2|Ω| is the real-valued angular frequency [2]. Equation (1.1a) is hyperbolic, but boundary condition (1.1b) is neither of Dirichlet-type or Neumann-type (it is a mixed-typed condition, sometimes called oblique condition). Hence, problem (1.1) is an ill-posed Cauchy problem [START_REF] Rieutord | Wave attractors in rotating fluids: a paradigm for ill-posed Cauchy problems[END_REF]. Nonetheless, it admits smooth polynomial eigenvectors when the geometry is a triaxial ellipsoid [START_REF] Vantieghem | Inertial modes in a rotating triaxial ellipsoid[END_REF][START_REF] Ivers | Enumeration, orthogonality and completeness of the incompressible Coriolis modes in a tri-axial ellipsoid[END_REF]. Moreover, these eigenvectors form a complete set in the Hilbert subspace of complex-valued divergenceless fields tangent to the boundary with the L 2 norm in ellipsoids [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF]. Given the formal analogy between rotation and stratification [START_REF] Veronis | The analogy between rotating and stratified fluids[END_REF], internal (gravity) modes also obey an ill-posed problem in non-rotating stratified fluids [START_REF] Rieutord | On the analogy between gravity modes and inertial modes in spherical geometry[END_REF]. In the presence of stratification and rotation, the mathematical problem becomes even more complicated because IGMs obey a mixed hyperbolic-elliptic operator [START_REF] Friedlander | Internal waves in a rotating stratified fluid in an arbitrary gravitational field[END_REF]. As for pure inertial [START_REF] Brunet | Linear and nonlinear regimes of an inertial wave attractor[END_REF][START_REF] Favier | Inertial wave super-attractor in a truncated elliptic cone[END_REF] and internal (gravity) waves [START_REF] Maas | Observation of an internal wave attractor in a confined, stably stratified fluid[END_REF], IGWs can also converge after multiple reflections on solid boundaries to wave attractors [START_REF] Pacary | Observation of inertia-gravity wave attractors in an axisymmetric enclosed basin[END_REF].

Attractors are interesting geometrical structures [START_REF] De Verdière | Attractors for two-dimensional waves with homogeneous Hamiltonians of degree 0[END_REF], 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 [START_REF] Allen | Some aspects of the initial-value problem for the inviscid motion of a contained, rotating, weakly-stratified fluid[END_REF] and Friedlander & Siegmann [START_REF] Friedlander | Internal waves in a contained rotating stratified fluid[END_REF] 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 [START_REF] Friedlander | Internal waves in a rotating stratified fluid in an arbitrary gravitational field[END_REF] 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 [START_REF] Vantieghem | Inertial modes in a rotating triaxial ellipsoid[END_REF][START_REF] Ivers | Enumeration, orthogonality and completeness of the incompressible Coriolis modes in a tri-axial ellipsoid[END_REF][START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF], while preserving polynomial eigenvectors with stratification. Finally, we characterise the wave spectrum using the mathematical theory recently presented in Colin de Verdière & Vidal [START_REF] De Verdière | The spectrum of the Poincaré operator in an ellipsoid[END_REF] 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).

Formulation of the general model (a) Primitive fluid dynamics equations

As illustrated in figure 2(a), we model a vortex by a fluid domain of volume V and (fictive) boundary ∂V , embedded within an ambient fluid (e.g. a planetary atmosphere or ocean) that is rotating at the angular velocity Ω = Ωs1 Ω with respect to an inertial frame (where 1 Ω is a unit vector). We work in the frame rotating at angular velocity Ω, and denote the position vector by r = (x, y, z) ⊤ using the Cartesian coordinates. The vortex is internally stratified in density under the action of gravity g. The vortex's stratification is quantified by the Brunt-Väisälä (BV) frequency N . Homogeneous fluids are such that N 2 = 0, whereas stably stratified fluids are characterised by N 2 > 0. As a starting point, we employ the Boussinesq approximation to account for density variations [START_REF] Spiegel | On the Boussinesq approximation for a compressible fluid[END_REF]. The vortex is differentially rotating with respect to the ambient fluid outside V , which can be modelled by a background flow U 0 . Laboratory experiments [10,11] and numerical calculations [9] show that, for incompressible fluids, U 0 is often close to a uniform-vorticity flow. The latter can be sought as U 0 = ω × r + ∇Ψ , where ω is the vortex's rotation vector (either cyclonic or anticyclonic) and Ψ is a scalar ensuring that the flow does not cross ∂V (e.g. [START_REF] Noir | Precession-driven flows in non-axisymmetric ellipsoids[END_REF] in an ellipsoid). The amplitude of U 0 is characterised by the Rossby number Ro = |ω|/Ωs, and we are interested in modelling stratified vortices in the regime |Ro| ≪ (N/Ωs) 2 . Thus, we neglect the effects of U 0 onto the background density and pressure. We further simplify the model by discarding baroclinic (and centrifugal effects), the very weak vortex's self-gravity, and the (weak) lateral variations of the planet's gravity at the size of a vortex. Hence, the fluid is only subject to the ambient gravity g = -g(z)1z within V , where 1z is the unit vector along the vertical z-axis. We denote by ρ * + ρ 0 (z) the reference density within V (with the mean density ρ * ), which is given by the hydrostatic equilibrium ∇P 0 = (ρ * + ρ 0 )g where P 0 is the hydrostatic pressure. Now, we seek small-amplitude perturbations around the reference state [U 0 , ρ 0 , P 0 ]. Without diffusion, the velocity v(r, t) obeys the linearised Euler equation given in the rotating frame by

∂ t v + (U 0 • ∇)v + (v • ∇)U 0 + 2Ω × v = -∇π + (ρ/ρ * ) g (2.1)
together with the incompressible condition ∇ • v = 0 for an incompressible fluid, and where π = p/ρ * is the reduced pressure. The density ρ(r, t) is then governed by

∂ t ρ/ρ * + U 0 • ∇ρ/ρ * = (v • 1z) N 2 (z)/g(z), (2.2)
where N 2 (z) = g • ∇(ρ 0 /ρ * ) 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 [START_REF] Lam | Supply mechanisms of the geostrophic mode in rotating turbulence: interactions with self, waves and eddies[END_REF] for a discussion in unbounded fluids). Second, it can be responsible for the onset of hydrodynamic instabilities. In particular, when U 0 is an elliptical flow [10,11], the so-called elliptical instability could develop [START_REF] Kerswell | Elliptical instability[END_REF] and be responsible for the transition towards a wave turbulence regime [START_REF] Reun | Inertial wave turbulence driven by elliptical instability[END_REF]. For mathematical simplicity, we assume U 0 = 0 in equations (2.1)-(2.2) to focus on the free wave motions given by

∂ t v + 2Ω × v = -∇π -N (z)ξ 1z, ∇ • v = 0, ∂ t ξ = (v • 1z) N (z), (2.3a-c)
where we have introduced ξ = (ρ/ρ * ) g(z)/N (z). The background flow will be reintroduced in future work (to investigate the outcome of the elliptical instability). Finally, equations (2.3a-c) are supplemented with boundary conditions (BCs). We impose that the fields are regular at r = 0. Then, realistic conditions would be to couple the dynamics inside the vortex with the exterior flow on ∂V . This would amount to considering continuity conditions across ∂V , and decaying solutions at infinity (e.g for Rankine [START_REF] Park | Instabilities and waves on a columnar vortex in a strongly stratified and rotating fluid[END_REF] or Lamb-Oseen [START_REF] Dizès | Inviscid waves on a Lamb-Oseen vortex in a rotating stratified fluid: consequences for the elliptic instability[END_REF] vortices). For mathematical simplicity, we here neglect the coupling with the exterior fluid and assume that ∂V is a rigid boundary. Consequently, the velocity obeys the no-penetration BC v • n| ∂V = 0 on ∂V , where n is the unit (outward) vector normal to the boundary. Without diffusion, no other physical BCs have to be explicitly enforced in the model. Indeed, the density and pressure BCs automatically follow from the velocity BC. In practice, the pressure satisfies a mixed BC that is given by the continuity of the normal component of equation (2.3a) on ∂V (see the ESM).

(b) Mixed hyperbolic-elliptic problem

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

[v, π, ξ] (r, t) = [u, Φ, ζ] (r) exp(iωt), (2.4)
where [u, Φ, ζ] are complex-valued amplitudes depending on space, and with ω ∈ C.

-ω 2 u + 2Ω × (iωu) + N 2 (z)(u • 1z)1z = -iω∇Φ, ∇ • u = 0, (2.5a,b)
Equations (2.5a,b) and the no-penetration BC u • n| ∂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 [START_REF] Friedlander | Internal oscillations in the Earth's fluid core[END_REF]. On the contrary, non-oscillatory modes with ω = 0 require a specific consideration. For unstratified fluids with N = 0, they are called geostrophic modes because they obey the geostrophic balance 2Ω × u = -∇Φ. When N (z) > 0, the steady modes are two-dimensional (2D) such that u • 1z = 0, and the density perturbation obeys the diagnostic equation ρ * (2Ω • ∇u) = ∇ζ × 1z (resulting from the thermal wind balance).

Because of historical reasons (dating back to Poincaré [START_REF] Poincaré | Sur l'équilibre d'une masse fluide animée d'un mouvement de rotation[END_REF]), the problem is generally recast into a scalar equation for the pressure Φ [START_REF] Friedlander | Internal waves in a rotating stratified fluid in an arbitrary gravitational field[END_REF][START_REF] Friedlander | Internal waves in a contained rotating stratified fluid[END_REF]. As shown in the ESM, the pressure is governed in V by a second-order equation whose higher-order terms are given by

[ω 2 -N 2 (z)]∇ 2 Φ + N 2 (z) (1z • ∇) 2 Φ -(2Ω • ∇) 2 Φ + • • • = 0 (2.6) when |ω| ̸ = ω ± (z) with 2ω 2 ± (z) = [N 2 (z) + 4|Ω| 2 ] ± [N 2 (z) + 4|Ω| 2 ] 2 -16N 2 (z)(Ω • 1z) 2 . (2.7)
The lower-order terms in equation (2.6), which vanish when N is constant, are here omitted for concision but are given in the ESM. When |ω| = ω ± (z), the pressure instead obeys a first-order equation (see the ESM). Finally, the pressure obeys a mixed BC on ∂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 [START_REF] Grubb | Distributions and operators[END_REF]. 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 [START_REF] Faure | Manifestation of the topological index formula in quantum waves and geophysical waves[END_REF][START_REF] Venaille | From ray tracing to waves of topological origin in continuous media[END_REF]. In particular, the nature of equation (2.6) is governed by a scalar quantity that is called the principal symbol. The latter is obtained from the equation by freezing the variable coefficients and keeping only the highest-degree terms. In practice, it amounts to transforming the derivatives as ∇ ↔ ik, where k = (k 1 , k 2 , k 3 ) ⊤ is a real-valued wave vector at the position (x, y, z). The principal symbol is then a homogeneous polynomial in k, which encapsulates many properties of the spectral problem [START_REF] De Verdière | Attractors for two-dimensional waves with homogeneous Hamiltonians of degree 0[END_REF][START_REF] De Verdière | Spectral theory of pseudodifferential operators of degree 0 and an application to forced linear waves[END_REF]. 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

|k| 2 ω 2 -N 2 (z)|1z × k| 2 + (2Ω • k) 2 = 0, (2.8)
which is the dispersion relation of IGWs in unbounded fluids when N is constant [3]. The wavelike equation is thus hyperbolic when equation (2.8) is satisfied for some non-zero wave vectors k ∈ R 3 . Otherwise, it is elliptic (except when ω = 0). Next, the bounds for the hyperbolic domain can be obtained from equation (2.8). The latter defines a quadric cone given by (g ⋆ ) ij k i k j = 0 (using Einstein notation) when ω ̸ = 0, where g ⋆ is the associated metric tensor. The nature of this quadric cone depends on the eigenvalues of g ⋆ , which are given by µ 0 = ω 2 -N 2 (z) and µ ± = ω 2 -ω 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 [START_REF] Friedlander | Internal waves in a rotating stratified fluid in an arbitrary gravitational field[END_REF]. The wave-like equation is thus hyperbolic in V when ω -(z)

< |ω| < ω + (z), whereas it is elliptic in V when 0 < |ω| < ω -(z) or |ω| > ω + (z).

(c) Spectral problem for the velocity

Further wave properties can be obtained by considering the velocity equation. We introduce the Hilbert space 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 ⟨u, v⟩ :=

V u † • v dV, (2.9)
where † is the complex conjugate. The associated norm is ||u|| = ⟨u, u⟩ 1/2 . Then, we denote by V 0 ⊂ V the closed subspace that is orthogonal, with respect to inner product (2.9), to the space spanned by vector fields made of gradients of smooth functions. A smooth element in V 0 is divergenceless and tangent to the boundary [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF], that is

V 0 := {u ∈ V | ∇ • u = 0 in V, u • n = 0 on ∂V }. (2.10)
Finally, we introduce the orthogonal projector L from V onto V 0 (called Leray projector after [START_REF] Leray | Sur le mouvement d'un liquide visqueux emplissant l'espace[END_REF]).

For any vector e ∈ V written as e = ∇Ψ + ∇ × A (using a Helmholtz decomposition) for some potentials [Ψ, A], the projected vector is defined by L(e) := e -∇Ψ . This projector is used in the mathematical analysis of incompressible flows [START_REF] Foias | Navier-Stokes equations and turbulence[END_REF], but also in some numerical studies [START_REF] Teed | Solenoidal force balances in numerical dynamos[END_REF]. Equipped with the above definitions, we can seek solutions of equation (2.5) with u ∈ V 0 in ellipsoids. We apply the orthogonal projector to the equation to remove the pressure term, and find that (ω, u) obey a quadratic eigenvalue problem (QEP). The latter is given by

Q ω (u) = 0, Q ω := -ω 2 I + ωiC + K, (2.11a,b)
where I is the identity operator, and with the two operators

iC(u) := iL[2Ω × u], K(u) := L[N 2 (z)(u • 1z)1z]. (2.12a,b)
The operator iC, called the Poincaré operator [START_REF] De Verdière | The spectrum of the Poincaré operator in an ellipsoid[END_REF], is a bounded self-adjoint operator in V 0 [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF].

The buoyancy operator K has properties given by Proposition 2.1.

Proposition 2.1. The operator K is bounded, self-adjoint and positive in V 0 .

Proof. We must show that K is symmetric, bounded, and positive in V 0 . We verify that K is symmetric since ⟨a, K(b)⟩ = ⟨K(a), b⟩ for any (a, b) ∈ V 0 × V 0 . Moreover, we have

||K(u)|| 2 = ⟨K(u), N 2 (u • 1z)1z⟩ -⟨K(u), ∇Φ⟩ = ⟨K(u), N 2 (u • 1z)1z⟩
since ⟨K(u), ∇Φ⟩ = 0 using the definition of L. Using Cauchy-Schwarz inequality, we then obtain

||K(u)|| 2 ≤ ||K(u)|| ||N 2 (u • 1z)1z|| ≤ max(N 2 )||K(u)|| ||u||, showing that K is bounded. Finally, K is positive because ⟨u, K(u)⟩ ≥ 0 if N 2 ≥ 0. ■
We want to characterise the spectrum σ of the quadratic spectral form Q ω in ellipsoids, which is the set of complex numbers ω for which Q ω is not continuously invertible. Actually, the spectral problem Q ω (u) = 0 can be converted into a standard eigenvalue problem (SEP) for a bounded self-adjoint operator in

V 0 × V 0 .
Proof. This is a general result for spectral families of the form Q ω := -ω 2 + ωA + B acting on a Hilbert space H, where [A, B] are bounded self-adjoint operators and with B ≥ 0 [START_REF] Barston | Eigenvalue Problem for Lagrangian Systems[END_REF][START_REF] Valette | Étude d'une classe de problèmes spectraux[END_REF]. When ω ̸ = 0, the QEP Q ω (u) = 0 can be recast as a SEP for the operators L or T given by

L := 0 I B A , T := A B 1/2 B 1/2 0 , (2.13a,b)
which are both acting on the Hilbert space H × H equipped with the inner product ⟨⟨z,

ζ⟩⟩ = ⟨z 1 , ζ 1 ⟩ + ⟨z 2 , ζ 2 ⟩ with z = (z 1 , z 2 ) ⊤ and ζ = (ζ 1 , ζ 2 ) ⊤ .
Here, T is a bounded self-adjoint operator (where B 1/2 is the square root of B, which is also self-adjoint [START_REF] Wouk | A note on square roots of positive operators[END_REF]), whereas L has no particular symmetries. The spectra of Q ω , L and T are identical outside ω = 0. ■ The spectrum σ(Q ω ) is thus the disjoint union of the point spectrum given by

σ P (Q ω ) := {ω ∈ C | ∃ u ̸ = 0 Q ω (u) = 0}, (2.14)
which is the subset spanned by the eigenvalues of QEP (2.11), and the remaining continuous spectrum (i.e. the set of ω for which T -ω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 [START_REF] De Verdière | Attractors for two-dimensional waves with homogeneous Hamiltonians of degree 0[END_REF]. These geometrical structures, called attractors [START_REF] De Verdière | Attractors for two-dimensional waves with homogeneous Hamiltonians of degree 0[END_REF], are associated with nearly singular (i.e. non-square-integrable) velocity fields when diffusion is vanishingly small [START_REF] Dintrans | Gravito-inertial waves in a rotating stratified sphere or spherical shell[END_REF]. Given the symmetries of the different operators, some general properties of the spectrum are given in Proposition 2.2.

Proposition 2.2. The spectrum of QEP (2.11) is real-valued and even. Moreover, the spectrum is bounded by |ω| ≤ ωmax with ωmax ≤ Ωs + Ω 2 s + max(N 2 ). Finally, two eigen-pairs (ω 1 , u 1 ) and (ω 2 , u 2 ) with ω 1 ̸ = ω 2 are not orthogonal with respect to inner product (2.9), but satisfy the two integral properties

ω 1 ω 2 ⟨u 2 , u 1 ⟩ = -⟨u 2 , K(u 1 )⟩, (ω 2 + ω 1 )⟨u 2 , u 1 ⟩ = ⟨u 2 , iC(u 1 )⟩.
(2.15a,b)

Proof. ω ∈ R result from the fact that the operator T defined in equation (2.13b) is self-adjoint when K ≥ 0. Then, the complex-conjugate of QEP (2.11) gives

-ω 2 u † -ωiC(u † ) + K(u † ) = 0,
which shows that (-ω, u † ) is an eigen-pair of the QEP. The upper bound of ωmax can be obtained using a min-max principle [START_REF] Barston | Eigenvalue Problem for Lagrangian Systems[END_REF]. Finally, the orthogonality conditions can be found using vector manipulations of the QEP [START_REF] Barston | Eigenvalue problem for Lagrangian systems[END_REF]. ■

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 ∂V is a triaxial ellipsoid of semi-axes [a, b, c], where the c-axis is aligned with gravity (figure 2b). The use of ellipsoids has a long-standing history in the study of geophysical vortices [14], and is supported by recent fluid dynamic experiments [10,11]. In the rotating frame, the boundary is expressed using the Cartesian coordinates as F (x, y, z) = 1 with F (x, y, z) := (x/a) 2 + (y/b) 2 + (z/c) 2 , such that the unit normal is n := ∇F/||∇F ||. Moreover, geophysical vortices with N ≫ 2Ωs are usually strongly flattened with a, b ≫ c (e.g. Jupiter's GRS is expected to have a depth of only a few hundreds of km [11]). In such vortices, the ambient gravity is expected to weakly varies with depth. Thus, we also assume that g = -g1z is uniform and consider a background density ρ * + ρ 0 (z) that varies linearly with z as

ρ 0 (z)/ρ * = -(N 2 /g) z, (3.1)
where N is a constant BV frequency. With these assumptions, the eigensolutions can be characterised using finite-dimensional arguments in ellipsoids.

(a) Polynomial description

Since an ellipsoid is a quadratic surface and expression (3.1) is a polynomial, we can seek eigenvectors using polynomial functions in the Cartesian coordinates. We denote by Pn : R 3 → 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 ≥ 1. Its dimension is N (Pn) = (n + 1)(n + 2)(n + 3)/6. Then, we denote by Vn the space of vectors in ellipsoids with coefficients in Pn and define the subspace of dimension N (V 0 n ) = n(n + 1)(2n + 7)/6, whose elements are divergenceless vector fields that are tangent to the boundary ∂V . Finally, we introduce the subspaces W 0

V 0 n := Vn ∩ V 0
n := (V 0 n-1 ) ⊥ ∩ V 0 n of dimension N (W 0 n ) = n(n + 2)
, which are spanned by vectors in V 0 n of degree n that are orthogonal to polynomial vectors of degree ≤ n -1. These polynomial spaces are key to characterise the spectrum because of Proposition 3.1.

Proposition 3.1. The spaces V 0 n with n ≥ 1 are invariant by the action of C and K in triaxial ellipsoids, that is

C |V 0 n ⊆ V 0 n and K |V 0 n ⊆ V 0 n .
Proof. The invariance of V 0 n under the action of C has been demonstrated elsewhere [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF][START_REF] De Verdière | The spectrum of the Poincaré operator in an ellipsoid[END_REF]. Here, we can use similar arguments to show that it is invariant by the action of K. Briefly, we must prove that the Leray projector L is well-defined on V 0 n . This is ensured by Theorem 2.4 in [START_REF] Axler | The Neumann problem on ellipsoids[END_REF], which shows that there is a unique scalar Ψ ∈ P n+1 in the definition of L. ■

We deduce from Proposition 3.1 that V 0 n is invariant by the action of Q ω . Moreover, since

⊕ n∈N * V 0 n is dense in V 0 [40]
, it shows that σ(Q ω ) is pure point in ellipsoids (i.e. σ = σ P ). The continuous spectrum is thus empty, and equation (2.5) admits square-integrable eigenvectors in V 0 n . 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 [START_REF] Pacary | Observation of inertia-gravity wave attractors in an axisymmetric enclosed basin[END_REF]). 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 [START_REF] Rabitti | Inertial wave rays in rotating spherical fluid domains[END_REF] but, since the continuous spectrum is also empty in a full sphere [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF][START_REF] Ivers | Enumeration, orthogonality and completeness of the incompressible Coriolis modes in a sphere[END_REF], these structures are not associated with singular global modes for vanishingly small diffusion. Here, because the spectrum is pure point, we can solve the eigenvalue problem by restricting u to V 0 n in ellipsoids.

(b) Galerkin algorithm

An accurate numerical description can be developed by seeking u ∈ V 0 n as u = j≥1 α j e j , where {α j } j≥1 are complex-valued coefficients, and {e j } j≥1 are real-valued basis vector elements of V 0 n [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF] normalised such that ||e j || = 1. Then, we substitute the above expansion into the QEP Q ω (u) = 0, and project the entire equation on every e i using inner-product (2.9) to cancel out the spatial dependence (Galerkin method [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF]). This gives the finite-dimensional QEP

-ω 2 M n + ωiCn + Kn α = 0 (3.2)
where α = (α 1 , α 2 , . . . ) ⊤ is the eigenvector, and where [M n, Cn, Kn] are three real-valued matrices of size

N (V 0 n ) × N (V 0 n ) and of elements M n,ij = ⟨e i , e j ⟩, C n,ij = ⟨e i , 2Ω × e j ⟩ and K n,ij = ⟨e i , N 2 (e j • 1z)1z⟩.
The matrix M n is positive definite but non-diagonal (because the chosen basis elements are not mutually orthogonal by construction in triaxial ellipsoids [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF]). The matrix Cn is antisymmetric, and Kn is a real-valued positive matrix. Further properties of [Cn, Kn] are given below in Proposition 3.2. Proposition 3.2. dim(ker Cn) = ⌈n/2⌉ and dim(ker Kn) = n(n + 1)(n + 2)/6.

Proof. When N = 0, ker(Cn) is associated with geostrophic modes of degree n. Its dimension was proved in [START_REF] De Verdière | The spectrum of the Poincaré operator in an ellipsoid[END_REF]. When N ̸ = 0, ker(Kn) is associated with 2D motions such that u • 1z = 0. The latter can be sought as u = ∇ × (Ψ 1z) where Ψ = p[1 -F (x, y, z)] is a stream function and with p ∈ P n-1 . The dimension of ker(Kn) is thus dim(P n-1 ) = n(n + 1)(n + 2)/6. ■

For numerical convenience, we convert the finite-dimensional QEP into an equivalent finitedimensional SEP when ω ̸ = 0. The finite-dimensional counterpart of the self-adjoint operator defined in equation (2.13b) is not well suited for numerical computations (because of the square root). Instead, we consider the SEP Lnz = ωz with z = (α, β) ⊤ and the matrix Ln given by

Ln := 0 I Bn iAn , An := M -1 n Cn, Bn := M -1 n Kn, (3.3a-c)
where I is the identity matrix. As illustrated in figure 3, the matrices An and Bn 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 Anα = ωα, can indeed be computed separately for every space W 0 n by seeking the eigenvalues of every diagonal block of An. Yet, since the matrix Ln has not an upper triangular block structure, we cannot compute separately the eigenvalues for every degree n when N ̸ = 0.

We have modified the bespoke numerical code SHINE, initially developed for compressible fluids [START_REF] Vidal | Compressible fluid modes in rigid ellipsoids: towards modal acoustic velocimetry[END_REF][START_REF] Vidal | Acoustic and inertial modes in planetary-like rotating ellipsoids[END_REF], to implement the above algorithm for Boussinesq fluids. The matrices [M n, Cn, Kn] are first built analytically. Since the basis vector elements are not mutually orthogonal in an ellipsoid [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF], the matrices become numerically ill-conditioned when the polynomial degree is large. Hence, we perform the computations using extended-precision algorithms. We can illustrate the Galerkin algorithm by considering the six polynomial solutions u ∈ V 0 1 . The latter are sought in the form of uniform-vorticity flows as [START_REF] Noir | Precession-driven flows in non-axisymmetric ellipsoids[END_REF]

u = α × r + ∇ψ, ψ = αx c 2 -b 2 b 2 + c 2 yz + αy a 2 -c 2 a 2 + c 2 xz + αz b 2 -a 2 a 2 + b 2 xy, (3.4a-b)
and with α = (αx, αy, αz) ⊤ . Next, the block matrices of L 1 are given by

B 1 =       N 2 c 2 b 2 + c 2 0 0 0 N 2 c 2 a 2 + c 2 0 0 0 0       , A 1 =        0 - 2Ωza 2 a 2 + c 2 2Ωya 2 a 2 + b 2 2Ωzb 2 b 2 + c 2 0 - 2Ωxb 2 a 2 + b 2 - 2Ωyc 2 b 2 +c 2 2Ωxc 2 a 2 + c 2 0        , (3.5a,b)
and with Ω = (Ωx, Ωy, Ωz) ⊤ . Finally, the eigenvalue-eigenvector solutions of L 1 can be solved analytically to extend the linear solutions obtained by Vantieghem [START_REF] Vantieghem | Inertial modes in a rotating triaxial ellipsoid[END_REF] when N = 0 and 1 Ω = 1z. However, we must rely on numerical computations to obtain Ln for higher degrees n.

(c) Enumeration of the eigenvalues

We can explicitly count the number of non-zero and zero eigenvalues of the QEP when u ∈ V 0 n . Since the spectra of the primitive equations, the one of Q ω and the one of L are identical when ω ̸ = 0, their finite-dimensional counterparts have the same number of non-zero eigenvalues. Proof. The number of non-zero eigenvalues can be estimated from the primitive equations

iωu + 2Ω × u = -∇π -ζN 1z, ∇ • u = 0, iωζ = (u • 1z) N.
Solutions of the latter equations belong to the space spanned by u ∈ V 0 n and ζ ∈ Pn, whose dimension is

N (V 0 n ) + N (Pn) = (n + 1)(n 2 + 4n + 2)/2.
Next, we count the number of zero eigenvalues. From the density equation, we deduce that ω = 0 is associated with 2D flows such that u • 1z = 0. Such 2D flows u ∈ V 0 n can be sought as u = ∇ × (Ψ 1z) where Ψ = p[1 -F (x, y, z)] is a stream function and with p ∈ P n-1 . Assuming that Ω = (Ωx, Ωy, Ωz) ⊤ , we obtain 2Ωz∂xΨ = -∂xπ and 2Ωz∂yΨ = -∂yπ from the horizontal components of the momentum equation. Then, we deduce that π = -2ΩzΨ + h(z), where h(z) is an arbitrary polynomial function of z and of degree n + 1 (with n ≥ 1). Finally, the density perturbation is obtained from the vertical component of the momentum equation, which gives ζN = 2(Ωy∂yΨ + Ωx∂xΨ ) -∂zπ. Thus, ω = 0 is associated with an eigenspace of dimension dim(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 Ln, which has 2N (V 0 n ) eigenvalues. Hence, the algebraic multiplicity of the zero eigenvalue ω = 0 is given by n(n + 1)(n + 2)/3. However, ω = 0 is not only associated with standard eigenvectors. The geometric multiplicity m L of ω = 0 is the dimension of ker Ln. The null space of the equation Ln(α, β) ⊤ = 0 is spanned by the vector elements such that β = 0 and Knα = 0. Therefore, we have ker Ln = ker Kn {0} and ω = 0 has the geometric multiplicity m L = n(n + 1)(n + 2)/6 according to Proposition 3.2. However, since Ln is not Hermitian, the eigenvalue ω = 0 has an algebraic multiplicity r L > m L . The missing r L -m L = n(n + 1)(n + 2)/6 polynomial solutions are associated with generalized eigenvectors, which belong to ker(L 2 n ).

Proposition 3.4. The eigenvalue ω = 0 of the matrix Ln 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 ω = 0 is thus partly spanned by generalized eigenvectors, which here satisfy Lnz ̸ = 0 and L 2 n z = 0. Back in the temporal domain, they correspond to solutions that are not steady but linear in t. However, these admissible solutions are mathematical artefacts due to the formulation of the quadratic problem in terms of the matrix Ln. Indeed, QEP (3.2) can be converted into an equivalent self-ajoint SEP [START_REF] Barston | Eigenvalue problem for Lagrangian systems[END_REF], which is the finite-dimensional counterpart of the self-adjoint operator defined in equation (2.13b). This shows that, from a physical viewpoint, ω = 0 is only associated with standard eigenvectors. Therefore, we will not further elaborate on the (artificial) generalized eigenvectors of Ln below. Following our prior investigation of the non-stratified problem [START_REF] De Verdière | The spectrum of the Poincaré operator in an ellipsoid[END_REF], we can combine polynomial computations for fixed values of n and spectral theory for n → ∞ to study the eigensolutions in ellipsoids.

Point spectrum in ellipsoids

We have shown that σ(Q ω ) 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 → ∞), it belongs to the discrete spectrum σ D . Otherwise, it belongs to the essential spectrum σ E := σ\σ D (i.e. the complement to the discrete spectrum). Since ⊕ n∈N * V 0 n is dense in V 0 , we conclude from Proposition 3.4 that the eigenvalue ω = 0 has an infinite multiplicity when n → ∞ (i.e. ω = 0 belongs to σ E ). Now, it remains to characterise the other eigensolutions.

When N is constant, the frequencies ω ± defined in equation (2.7) are constant in the volume. Moreover, it is possible to find a better estimate of ωmax (compared to the one given in Proposition 2.2). In the aligned case, it was shown in [START_REF] Friedlander | Internal waves in a contained rotating stratified fluid[END_REF] that ωmax = ω + with ω + = max(N, 2Ωs) . Thus, σ 3 is empty in ellipsoids when Ω and 1z are collinear. We also show below in Proposition 4.1 that ωmax = ω + , even when Ω and 1z are not collinear. Proposition 4.1. When N is constant, the pure point spectrum is bounded by ω + in an ellipsoid.

Proof. We introduce the projector operator P acting on a pair z = (u, ζ) ⊤ by P(z) := (Lu, ζ) ⊤ , where L is the Leray projector introduced in §2. Then, we can reformulate the spectral problem as the spectral properties of the linear operator P defined by P(z) := P(EP(z)) with

E = i      0 -2Ω 3 2Ω 2 0 2Ω 3 0 -2Ω 1 0 -2Ω 2 2Ω 1 0 N 0 0 -N 0     
.

P is a self-adjoint operator for square-integrable fields. We have ||P|| = 1, and the eigenvalues of the Hermitian matrix E are ±ω ± such that ||E|| ≤ ω + . Thus, the spectrum is bounded by ω + . ■ Therefore, we can divide the allowable frequency range into disjoint intervals with

I 1 = [-ω + , -ω -] ∪ [ω -, ω + ] and I 2 =] -ω -, ω -[ \ {0}
. The velocity equation is hyperbolic in V when ω ∈ I 1 , and elliptic in V when ω ∈ I 2 . We also define the associated spectra σ i = σ P ∩ I i for i = {1, 2} and, when u is restricted to V 0 n , the subsets are denoted by [σ 1,n , σ 2,n ] such that σ i,n-1 ⊆ σ i,n . In unbounded fluids, a continuum of IGWs exists with angular frequencies ω ∈ I 1 given by dispersion relation (2.8). We will show below that the spectrum σ 1 ⊂ σ E is dense in I 1 within ellipsoids. On the contrary, σ 2 is known to be empty in unbounded fluids (which agrees with prior experimental observations [START_REF] Peacock | The effect of rotation on conical wave beams in a stratified fluid[END_REF]). This is consistent with the fact that the velocity satisfies an elliptic-type equation when ω ∈ 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 |ω| < ω - [START_REF] Allen | Some aspects of the initial-value problem for the inviscid motion of a contained, rotating, weakly-stratified fluid[END_REF][START_REF] Friedlander | Internal waves in a contained rotating stratified fluid[END_REF]. The eigensolutions with ω ∈ σ 1 are referred as class-I solutions, and those with ω ∈ σ 2 are called class-II solutions in [START_REF] Allen | Some aspects of the initial-value problem for the inviscid motion of a contained, rotating, weakly-stratified fluid[END_REF][START_REF] Friedlander | Internal waves in a contained rotating stratified fluid[END_REF]. We describe the two spectra below.

(a) Low-frequency (surface) modes (i) Illustrative examples

Contrary to the free-space case, it is known that σ 2 is not empty when Ω ∝ 1z in some bounded geometries (e.g. in cylinders or spheres [START_REF] Kerswell | Elliptical instabilities of stratified, hydromagnetic waves[END_REF][START_REF] Friedlander | Internal waves in a contained rotating stratified fluid[END_REF]). In triaxial ellipsoids, we can also show that σ 2 ̸ = ∅ when Ω • 1z ̸ = 0 (otherwise, ω -= 0 and σ 1 fills the entire interval 0 < |ω| < ω + ). To this end, we first consider the eigenvalues for u ∈ V 0 1 , which are the roots of the characteristic equation det(L 1 -ωI) = 0. The six roots ω j of the characteristic polynomial are the two degenerate solutions ω 1,2 = 0 and the four solutions

ω 3,4 = ± β 4 2∆ + (β 4 /∆) 2 -4β 2 /∆ 2 , ω 5,6 = ± β 4 2∆ - (β 4 /∆) 2 -4β 2 /∆ 2 , (4.1a-b)
where we have introduced the positive coefficients

β 4 = N 2 c 2 (a 2 + b 2 )(a 2 + b 2 + 2c 2 ) + 4Ω 2 x b 2 c 2 (b 2 + c 2 ) + 4Ω 2 y a 2 c 2 (a 2 + c 2 ) + 4Ω 2 z a 2 b 2 (a 2 + b 2 ), β 2 = c 4 N 2 [4(Ω 2 x b 2 + Ω 2 y a 2 ) + N 2 (a 2 + b 2 )], ∆ = 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 ω 3,4 = ± β 4 /∆, which reduce to the eigenvalues of the pure inertial modes in V 0 1 obtained by Vantieghem [START_REF] Vantieghem | Inertial modes in a rotating triaxial ellipsoid[END_REF] when Ωx = Ωy = 0. Solutions (4.1) are illustrated in figure 5(a) for a sphere and a flattened (triaxial) ellipsoid. Several points are worth commenting on in the figure. First, ω 3 and ω 5 smoothly vary as a function of N . We observe that ω 3 ∈ σ 1,1 for any value of N but, on the contrary, ω 5 can either belong to σ 1,1 or σ 2,1 when N is varied. For instance, the transition occurs here at N/Ωs > 4 in a sphere and at a much larger value N/Ωs ≳ 25 in the pancake-like ellipsoid.

Next, we show in figure 5(b) the evolution of dim(σ 2,n ) as a function of n for some illustrative cases. We see that we have always σ 2,n ̸ = ∅ for a sufficiently large degree n, and that dim(σ 2,n ) increases when n is increased. This suggests that σ 2 has an infinite (but countable) number of eigenvalues when n → ∞. However, dim(σ 2,n ) is found to strongly depend on the geometry and the value of N/Ωs. More specifically, dim(σ 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 ≪ 1. Finally, the number of eigenvalues seems to be bounded by n(n + 1), which corresponds to the number of spherical harmonics Y m l of degree 2 ≤ l ≤ n + 1 and order 0 < |m| < l. This shows that σ 2 is related to solid (spherical) harmonics.

Further insight into the link between σ 2 and the theory of spherical harmonics can be gained by considering the aligned case 1 Ω = 1z with N = 2Ωs. For this particular configuration, we have ω -= ω + from equation (2.7). Then, pressure equation (2.6) and its associated BC for eigenvalues ω ∈ σ 2 reduce to

∇ 2 Φ = 0, n • ∇Φ = - 1 iω (2Ω × n) • ∇Φ on ∂V . (4.2a,b)
Laplace equation (4.2a) can then be solved using harmonic functions in spherical or ellipsoidal geometries, and the eigenvalues ω are the roots of a transcendental equation obtained from BC (4.2b). In a sphere, the pressure associated with u ∈ V 0 n is given in spherical coordinates (r, θ, ϕ) by the solid harmonics Φ(r) = r l Y m l (θ, ϕ) with ω = ±2Ωsm/l, where Y m l is the spherical harmonic of degree l ≤ n + 1 and azimuthal order m < l. The modes with m > 0 have a negative frequency ω < 0 when Ωs > 0 (respectively positive ω > 0 if m < 0). Given our phase convention ∝ exp[i(ωt + mϕ)], these modes have a positive azimuthal phase velocity c ϕ = -ω/m > 0 when Ωs > 0 (i.e. prograde direction). The explicit solutions show that the spectrum σ 2,n is here isomorph to the set of rational numbers Q, such that dim(σ 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 [START_REF] Friedlander | Internal waves in a contained rotating stratified fluid[END_REF]. Since Q is dense in R, σ 2 belongs to the essential spectrum σ E when 1 Ω = 1z and N = 2Ωs.

(ii) General properties

The above examples strongly suggest that σ 2 ⊂ σ E . This can be shown by further looking at the pressure problem. Assuming that |ω| ̸ = ω ± , the pressure is governed by a boundary-value problem given by equation (2.6) in V and a mixed-type BC on ∂V . This boundary-value scalar problem is elliptic if the equation is elliptic in V (which occurs here when 0 < |ω| < ω -), but also if some invertibility conditions are satisfied by the BC on ∂V [START_REF] Grubb | Distributions and operators[END_REF] (see the ESM). If both requirements were met in V and everywhere on ∂V in a given interval [ω 1 , ω 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 σ 2 belongs to the essential spectrum in an ellipsoid when Ω and 1z are aligned (see the ESM). We also conjecture that σ 2 ∈ σ E in the misaligned case. Because the spectrum is pure point in ellipsoids, an infinite number of modes would thus exist with 0 < |ω| < ω -when n → ∞.

Finally, we explore the properties of the corresponding eigenvectors. We remind the reader that these modes only exist when ω -̸ = 0. Otherwise, pure inertial modes (when N = 0), pure internal (gravity) modes (when Ω = 0) or pure IGMs (when Ω • 1z = 0) fill the entire interval 0 < |ω| < ω + . Two modal families are shown in figure 6. When N/Ωs ≤ 2, the energy of the lowest-frequency modes with |ω| → 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/Ωs ≪ 1. As investigated by Allen [START_REF] Allen | Some aspects of the initial-value problem for the inviscid motion of a contained, rotating, weakly-stratified fluid[END_REF] when Ω ∝ 1z, 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 |ω| ≲ ω -have an energy that is maximum near the boundary and decays more or less significantly away from the boundary. These lowfrequency (surface) modes seem similar to the (trapped) Kelvin waves in shallow-water models [START_REF] Thomson | On gravitational oscillations of rotating water[END_REF]. 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 ∂V (see figure 6). Therefore, the presence of a boundary introduces a new subset of low-frequency (surface) modes in a rotating stratified fluid (which differ from the higher-frequency IGMs).

(b) High-frequency IGMs

We investigate IGMs with angular frequencies belonging to σ 1 . A few illustrative modes are shown in figure 7. IGMs are mainly inertial modes modified by gravity when N ≪ 2Ωs, whereas they are gravity modes more or less affected by rotation when N ≫ 2Ω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 ω/Ωs → 2 and 1 Ω = 1z. From dispersion relation (2.8), the angle ϑ between 1z and k is such that ϑ → 0 when ω/Ωs → 2. Since the group velocity is orthogonal to k for IGWs [3], the energy propagates in the horizontal direction (in agreement with the observed pancake-like structure). On the contrary, we have θ → π/2 for high-frequency IGWs with ω → N , such that the energy varies in the vertical direction. IGMs with ω ∈ σ 1 are thus the direct counterparts in ellipsoids of IGWs in unbounded fluids.

Such as IGWs without boundaries, IGMs are dense in this frequency interval when n → ∞. 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 σ 1 can be obtained using mathematical analysis when n → ∞, and then compared to the density obtained from the polynomial computations at fixed values of n. We define the probability measure πn of the eigenvalues πn := (1/dn) dn j=1 δ(ω j ) with dn = 2N (W 0 n ), where δ is the Dirac function. As explained above, the eigenvalues of SEP (3.3) cannot be easily separated in different subsets for every degree n. Therefore, we rather consider the joint repartition of the eigenvalues of Q ω | V 0 n . Naturally, the two repartitions have the same asymptotic distribution when n → ∞. 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.

{ω j } of Q ω | W 0 n given by
Proposition 4.2. The measure πn converges weakly to the asymptotic measure π∞ when n → ∞. The asymptotic measure is characterised by a cumulative distribution function given by

λ -∞ dπ∞ = 1 8π Area(S λ ∩ S 2 )
where S 2 is the surface of the unit sphere, S λ is the surface in R 3 defined by

S λ := { k ∈ R 3 | ω( k) ≤ λ}, k = ( A 1 kx, A 2 ky, A 3 kz) ⊤ , with A 1 = 1/a 2 , A 2 = 1/b 2 and A 3 = 1/c 2 ,
and where ω( k) is the plane-wave dispersion relation given by equation (2.8) for the rescaled wave vector k. Hence, the asymptotic measure has a probability density function f∞ := dπ∞/dω that is even and non-vanishing when

ω 2 -≤ ω 2 ≤ ω 2 + .
Proof. This is the most technical part of the paper. The proof closely follows the strategy presented in Colin de Verdière & Vidal [START_REF] De Verdière | The spectrum of the Poincaré operator in an ellipsoid[END_REF] 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 T defined in equation (2.13b). Thus, the starting point is to consider the matrix operator L defined in equation (2.13a). The fact that this operator is not self-adjoint is not a problem in the analysis (i.e. in the trace formula). Then, we apply the tools of microlocal analysis in the presence of boundary effects. They allow us to get the asymptotics of dn j=1 f (ω j ) when n → ∞ for any function f ∈ C 0 (R \ 0) (i.e. the set of real-valued continuous functions vanishing near zero). The analysis heavily relies on the principal symbol of the operator L when ω ̸ = 0. This gives the asymptotic measure in the interval I 1 (where the velocity equation is hyperbolic). Note that we avoid the value ω = 0 because ker 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 σ 1 is dense in I 1 when n → ∞. 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 ∞ -∞ dπ∞ = 1 but, since Proposition 4.2 heavily relies on the dispersion relation of IGWs, it is assumed that dπ∞/dω = 0 when ω ∈ I 2 . To have a correct comparison, we have thus removed the spectra σ 2,n before normalising the density. An excellent agreement is then obtained with the asymptotic density (even for the moderate polynomial degree n = 35 considered here). We have also represented a particular configuration for which 1 Ω • 1z = 0, such that ω -= 0 from expression (2.7). In this case, IGMs are dense in the entire interval ]0, ω + [. Note that a similar agreement is found for other configurations (not shown). Hence, Proposition 4.2 can be used to obtain the density of σ 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Ωs, 2Ωs[ for a sphere [START_REF] De Verdière | The spectrum of the Poincaré operator in an ellipsoid[END_REF]. 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 lowfrequency modes). Between these two limits, the density is very close to the non-rotating case when N/Ωs ≫ 1. On the contrary, the density of the high-frequency IGMs with |ω|/ωmax → 1 approaches the density of the non-stratified modes when N/Ωs ≤ 2. We also observe that the geometry does modify the spectral density. When N/Ωs < 2, the density is maximal near |ω| → ω -in elongated (i.e. prolate) ellipsoids with c/a ≫ 1, whereas it is maximal near |ω| → ω + in flattened (i.e. oblate) ellipsoids with c/a ≪ 1. Instead, low-frequency IGMs with |ω| → ω -are favoured in flattened ellipsoids when N/Ωs > 2.

Finally, the situation ω + = ω -occurring when 1 Ω = 1z and N/Ωs = 2 deserves a specific consideration. The spectrum σ 1 is then reduced to the degenerate eigenvalue ω/Ωs = ±2, whose multiplicity can be determined. Indeed, we know that dim(σ 2,n ) = n(n + 1) from the explicit solutions of equations (4.2a,b). Because the total number of non-zero eigenvalues is fixed for a degree n (Proposition 3.3), we have dim(σ 1,n ) = n(n + 1)(n + 2)/3 when ω + = ω -.

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 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Ω • 1z = 2Ωz cos(θ), where θ the mean colatitude of the vortex (figure 2).

f [rad.s -1 ] N/f a [km] b [km] c [km]
Jupiter's GRS ∼ 1.4 × 10 -4 O( 102 ) ∼ 9000 ∼ 5000 ∼ 100 Jupiter's Oval BA ∼ 1.9 × 10 -4 O( 102 ) ∼ 3000 ∼ 3000 ∼ 50 Mediterranean eddies ∼ 8.5 × 10 -5 10 -50 10 -10 2 10 -10 2 0.5 -1 with N/Ωs ≫ 1 and, consequently, they usually have a very small aspect ratio c/a ≤ 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 [START_REF] Lanchon | Internal wave turbulence in a stratified fluid with and without eigenmodes of the experimental domain[END_REF]). The bounds given by formula (2.7), which depend on the orientation of Ω and on the ratio N/Ωs, can be simplified in the relevant regime N/Ωs ≫ 1. It gives

ω 2 + ≃ N 2 + 4 |Ω| 2 -(Ω • 1z) 2 , ω 2 -≃ 4|Ω • 1z| 2 . (5.1a,b)
We see that the upper bound approaches ω 2 + → N 2 for realistic values N/Ωs ≫ 1. Detecting IGMs within a vortex may thus yield an estimate of the vortex's stratification.

A striking property of the global modes is that they admit, under our assumptions, polynomial solutions. This behaviour directly results from Proposition 3.1, and IGMs in other configurations can be very different. To illustrate this point, we consider a sphere of radius R that is stratified under the radial gravity g = -(g/R) r. The background hydrostatic density is then ∇ρ 0 = -(β/R) r (where β > 0 is a constant). The BV frequency is thus given by N 2 = N 2 0 (r/R) 2 using the Boussinesq approximation, where N 0 = βg/ρ * is the constant value of N at r = R. Then, the IGMs obey a mixed hyperbolic-elliptic equation that changes of type when [START_REF] Friedlander | Internal waves in a rotating stratified fluid in an arbitrary gravitational field[END_REF][START_REF] Dintrans | Gravito-inertial waves in a rotating stratified sphere or spherical shell[END_REF]

s 2 /(ω 2 -4Ω 2 s ) + z 2 /ω 2 = 1/N 2 0 , (5.2)
where s = x 2 + y 2 is the cylindrical radius measured from the rotation axis Ω = Ωs1z. For a given angular frequency ω, equation (5.2) defines a turning surface that separates regions where the equation is hyperbolic and elliptic in V . Such turning surfaces strongly impact the spatial structure of IGMs, which are mainly trapped in the hyperbolic region [START_REF] Dintrans | Gravito-inertial waves in a rotating stratified sphere or spherical shell[END_REF]. Moreover, IGMs with |ω| → 2Ωs < N become trapped in the equatorial waveguide when N/Ωs ≫ 1 [START_REF] Stewartson | On waves in a thin shell of stratified rotating fluid[END_REF]. Finally, such turning surfaces are also likely responsible for a non-empty continuous spectrum without diffusion [START_REF] Dintrans | Gravito-inertial waves in a rotating stratified sphere or spherical shell[END_REF]. On the contrary, the eigenmodes in a stratified ellipsoid with a constant N are always smooth and can penetrate deep inside the volume.

We have also shown the existence of low-frequency (surface) modes when 0 < |ω| < ω -. Their properties are very similar to those of the (trapped) Kelvin waves in rotating shallow-water models [START_REF] Thomson | On gravitational oscillations of rotating water[END_REF], which owe their existence to topological features [START_REF] Tauber | A bulk-interface correspondence for equatorial waves[END_REF][START_REF] Venaille | Wave topology brought to the coast[END_REF]. 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,[START_REF] Bars | Numerical study of the McIntyre instability around Gaussian floating vortices in thermal wind balance[END_REF] τν = (a 2 /ν)(2Ωs/N ) 2 , where ν 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 [START_REF] Rieutord | On the analogy between gravity modes and inertial modes in spherical geometry[END_REF][START_REF] Dintrans | Gravito-inertial waves in a rotating stratified sphere or spherical shell[END_REF]. Hence, we could wonder how much energy could be dissipated if some IGMs were forced inside an elliptical vortex. Without stratification, a

V u † • ∇ 2 u dV = V u † • L(∇ 2 u) dV = 0.
(5.3)

This could be interpreted as a vanishing viscous dissipation of the inertial modes in an ellipsoid if there were no boundary layers. However, when viscous BCs are properly taken into account on the boundary (either with no-slip or stress-free BCs), the inertial modes have a non-zero viscous dissipation [2,[START_REF] Vidal | Precession-driven flows in stress-free ellipsoids[END_REF]. From a mathematical viewpoint, formula (5.3) is valid in an ellipsoid because (i) L(∇ 2 u) ∈ V 0 n-2 when u ∈ V 0 n and (ii) two inertial eigenvectors (u i , u j ) are orthogonal with respect to inner product (2.9) such that ⟨u i , u j ⟩ = δ ij (when properly normalised). Here, the IGMs do not satisfy property (5.3) because they are not mutually orthogonal with respect to inner product (2.9) as explained in Proposition 2.2. This shows that the IGMs are (at least) subject to viscous dissipation in the bulk on the slow time scale 1/E in a stratified vortex. In addition, they will also be affected by thermal diffusion (if the stratification has a thermal origin, otherwise by molecular diffusion). Since the IGMs are here smooth polynomial solutions, it is possible to account for the leading-order dissipative effects using boundary-layer theory (see the ESM). This also shows that the IGMs would be dissipated on a typical time scale that is of the same order as τν 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 [START_REF] Kerswell | Elliptical instability[END_REF][START_REF] Lanchon | Internal wave turbulence in a stratified fluid with and without eigenmodes of the experimental domain[END_REF][START_REF] Boury | Triadic resonant instability in confined and unconfined axisymmetric geometries[END_REF]. Therefore, IGMs may play a role in the transition towards turbulence and the long-term evolution of stratified vortices.

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 [START_REF] De Verdière | The spectrum of the Poincaré operator in an ellipsoid[END_REF]. The normal modes consist of high-frequency IGMs and low-frequency (surface) modes. Using microlocal analysis, we have shown that the boundary plays a key role in sustaining such low-frequency (surface) modes. Finally, we have discussed our results in the light of geophysical applications, arguing that accounting for IGMs in the models could be useful to (better) understand the dynamics of stratified vortices.

Several questions have remained unanswered by our work (even in the linear theory). A more in-depth investigation of the properties of the low-frequency (surface) modes will be presented in a mathematical work. From a physical viewpoint, the model could also be extended to account for additional ingredients (while keeping polynomial solutions). For instance, a uniform-vorticity background flow U 0 will be taken into account to study the effects of the vortex's differential rotation (as measured by the Rossby number Ro). Large-scale vortices (e.g. Jupiter's vortices or Mediterranean eddies) are expected to be in the low-Rossby regime |Ro| ≲ O(10 -1 ), whereas smaller-scale vortices are in the opposite regime |Ro| ≳ O(1) (e.g. vortices at sub-mesoscales ≲ 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 [START_REF] Kerswell | Elliptical instabilities of stratified, hydromagnetic waves[END_REF][START_REF] Friedlander | Hydromagnetic waves in the Earth's fluid core[END_REF]. Moreover, there is a class of magnetic fields leaving invariant the three-dimensional polynomial vector spaces in an ellipsoid [START_REF] Malkus | Hydromagnetic planetary waves[END_REF][START_REF] Gerick | Pressure torque of torsional Alfvén modes acting on an ellipsoidal mantle[END_REF]. 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 [START_REF] Cébron | Tidal instability in a rotating and differentially heated ellipsoidal shell[END_REF][START_REF] Reun | Parametric instability and wave turbulence driven by tidal excitation of internal waves[END_REF][START_REF] Vidal | Magnetic fields driven by tidal mixing in radiative stars[END_REF][START_REF] Onuki | Breaking of internal waves parametrically excited by ageostrophic anticyclonic instability[END_REF]. The polynomial Galerkin algorithm presented in this work can be extended in the linear theory to determine the onset of the instabilities upon an elliptical flow. Therefore, elliptical instabilities will be considered in future work to explore, as a function of the Rossby number, the transition towards turbulence in stratified vortices. Data Accessibility. The paper has an Electronic Supplementary Material (ESM). The source code is available at https://bitbucket.org/vidalje/.

Authors' Contributions. This work is an idea of JV to pursue the interdisciplinary collaboration with YCdV after their joint work [START_REF] De Verdière | The spectrum of the Poincaré operator in an ellipsoid[END_REF]. JV designed the study and performed the Galerkin computations, whereas YCdV conducted the microlocal analysis. Both authors discussed the results presented in the paper and drafted the manuscript before giving final approval for submission.

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

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

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

Proof. We recast Euler equation (1a) and buoyancy equation (2) as a linear system given by

M     u x u y u z ζ/ρ *     = -     ∂ x Φ ∂ y Φ ∂ z Φ 0     , M :=     iω -2Ω z 2Ω y 0 2Ω z iω -2Ω x 0 -2Ω y 2Ω x iω g(z) 0 0 -N 2 (z)/g(z) iω     , (8a,b)
and with u = (u x , u y , u z ) ⊤ . The system admits non-zero solutions when det(M ) := D ω (z) ̸ = 0. In this case, we can invert the linear system and the result follows. Note that it reduces to equation (2.7.2) in Greenspan [1] when N = 0. ■ Proposition 1.2. When D ω (z) ̸ = 0, the pressure obeys the second-order equation given by

∇ • (f /D ω ) = 0 ⇐⇒ ∇ • f -(1/D ω ) f z ∂ z D ω = 0, (9)
with

∂ z D ω = -(ω 2 -ω 2 -)∂ z ω 2 + -(ω 2 -ω 2 + )∂ z ω 2 -.
When N (z) is variable and Ω = Ω z 1 z , equation ( 9) reduces to

∇ 2 Φ + ∂ ∂z N 2 -4Ω 2 z ω 2 -N 2 ∂Φ ∂z = 0, (10)
which is identical to equation (5.1) in Friedlander & Siegmann [2]. For oscillatory motions with ω ̸ = 0, equation ( 11) is thus equivalent to equation (3) when |ω| ̸ = ω ± (z). However, the pressure must obey another equation D ω (z) = 0 (see below).

Proposition 1.3. When N is constant, equation ( 9) reduces to P ω (Φ) = 0 with

P ω := -∇ 2 - N 2 ω 2 -N 2 (1 z • ∇) 2 + 4 ω 2 -N 2 (Ω • ∇) 2 . ( 11
)
The operator P ω is elliptic when 0 < |ω| < |ω -, and its principal symbol p is given by

p := Ak 2 x + Bk 2 y + Ck 2 z + 2Dk y k z + 2Ek x k z + 2F k x k y (12)
with

A = ω 2 -N 2 -4Ω 2 x ω 2 -N 2 , B = ω 2 -N 2 -4Ω 2 y ω 2 -N 2 , C = ω 2 -4Ω 2 z ω 2 -N 2 , D = -4Ω y Ω z ω 2 -N 2 , E = -4Ω x Ω z ω 2 -N 2 , F = -4Ω x Ω y ω 2 -N 2 .
Finally, it is worth explicitly writing the Cartesian components of f = (f x , f y , f z ) ⊤ when N is constant (as we will reuse them below). They are given by

f x ω 2 -N 2 = iω[A∂ x Φ + F ∂ y Φ + E∂ z Φ] - 2ω 2 Ω y ω 2 -N 2 ∂ z Φ + 2Ω z ∂ y Φ, (13a) f y ω 2 -N 2 = iω[F ∂ x Φ + B∂ y Φ + D∂ z Φ] + 2ω 2 Ω x ω 2 -N 2 ∂ z Φ -2Ω z ∂ x Φ, (13b) f z ω 2 -N 2 = iω[E∂ x Φ + D∂ y Φ + C∂ z Φ] + 2ω 2 ω 2 -N 2 (Ω y ∂ x Φ -Ω x ∂ y Φ). ( 13c
)
where the coefficients (A, B, C, D, E, F ) appear in the principal symbol of P ω (Proposition 1.3).

Proposition 1.5. When D ω (z) = 0 (i.e. on turning surfaces), the pressure is given by the firstorder equation

(2Ω • 1 z )(2Ω • ∇Φ) + 2iω Ω • (1 z × ∇Φ) -ω 2 (1 z • ∇Φ) = 0, (16)
which is a particular case of equation (3.5) in Friedlander & Siegmann [2].

Proof. We apply 2Ω× to equation (1a), which gives

iω(2Ω × u) = -2Ω × (2Ω × u) -2Ω × ∇Φ -g(z)(ζ/ρ * )(2Ω × 1 z ), = -4 (Ω • u)Ω -|Ω| 2 u -2Ω × ∇Φ -g(z)(ζ/ρ * )(2Ω × 1 z ).
The latter equation can be rearranged to give

4|Ω| 2 -ω 2 u = -iω∇Φ + 4(Ω • u)Ω + 2Ω × ∇Φ + N 2 (z) iω (1 z • u) [2Ω × 1 z -iω1 z ] , (17)
where we have used equation (1a) to simplify 2Ω × u and used equation (2a) to simplify the buoyancy term. Thus, it only remains to express Ω • u and 1 z • u as a function of Φ only. We take the scalar product of equation (1a) with iωΩ, which gives

ω 2 (Ω • u) = iω(Ω • ∇Φ) + N 2 (z)(Ω • 1 z )(1 z • u). (18)
Then, we take the scalar product of equation (1a) with iω1 z to obtain

N 2 (z) -ω 2 (1 z • u) = iω 2Ω • (1 z × u) -iω(1 z • ∇Φ). (19)
Since we have also from equation (1a)

iω(1 z × u) = -1 z × (2Ω × u) -1 z × ∇Φ = -2 [(1 z • u)Ω -(1 z • Ω)u] -1 z × ∇Φ, (20)
we combine equations ( 18) to (20) to obtain

-D ω (z)(1 z • u) = iω(2Ω • 1 z )(2Ω • ∇Φ) -2ω 2 Ω • (1 z × ∇Φ) -iω 3 (1 z • ∇Φ). (21)
From equation ( 21), we see that the pressure must satisfy equation ( 16) when D ω (z) = 0. ■

Microlocal analysis of boundary-value problems

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

Elliptic boundary-value problem

Remark. In this subsection, the Cartesian coordinates are denoted by (x 1 , x 2 , x 3 ) to make use of index notations.

We consider a smooth bounded domain V in R d , and a symmetric elliptic differential operator P ω of order 2 with smooth coefficients and of principal symbol p. We also assume that the Dirichlet extension of P ω is > 0. The scalar Φ is governed by the boundary-value problem given by P ω (Φ) = 0 in V , B(Φ) := n ⋆ • ∇Φ + V(Φ) = 0 on ∂V , (22a,b)

where n ⋆ is the unit vector normal computed with respect to the metric g induced on ∂V by the principal symbol of P ω (see Box 2.1), and V is a differential operator of degree 1 tangent1 to ∂V .

To solve system [START_REF] Bars | Thermo-elliptical instability in a rotating cylindrical shell[END_REF], we introduce the Dirichlet-to-Neumann operator N in Proposition 2.1. Note that N is a self-adjoint operator acting on functions defined in the space L2 (∂V, dσ) for the area dσ defined by the metric g restricted on ∂V [5]. To solve problem [START_REF] Bars | Thermo-elliptical instability in a rotating cylindrical shell[END_REF], we first solve the boundary problem N (Ψ)+V(Ψ) = 0 on ∂V , and then consider the Dirichlet problem in V given the boundary field Ψ. In practice, this amounts to solving the problem on the boundary ∂V and, then, propagating the boundary solution in the interior V by using the Dirichlet problem. Actually, the nature of the boundary-value problem crucially depends on the nature of the boundary equation N (Ψ) + V(Ψ) = 0 on ∂V . This can be analysed using microlocal analysis. Indeed, N is a nice pseudo-differential operator of principal symbol is given in Proposition 2.2. The ellipticity of boundary-value problem [START_REF] Bars | Thermo-elliptical instability in a rotating cylindrical shell[END_REF] is then given by Proposition 2.3. Finally, the nature is given by Proposition 2.4 below.

Proposition 2.2. The principal symbol of the Dirichlet-to-Neumann operator N , defined for an elliptic operator with a positive principal symbol p > 0, is given by g ⋆ | ∂V , where g ⋆ | ∂V is the dual of the restriction on ∂V of the metric g induced by p in V .

Proof. It is given in the Appendices of Taylor's book [6]. Note that the principal symbol is defined as --g| ∂V in [6]. The first minus sign comes from the fact that an inward vector n ⋆ was considered, and the second from the fact that the other convention p ≤ 0 was used. ■ Proposition 2.3. Let us denote by v the principal symbol of V. Then, boundary-value problem [START_REF] Bars | Thermo-elliptical instability in a rotating cylindrical shell[END_REF] is elliptic if and only if the principal symbol g ⋆ | ∂V + v is elliptic on ∂V , that is if it does not vanish for non-zero covectors k ∈ R 2 on the boundary ∂V .

Proof. See Chapter 9 in [7]. This is also known as the Shapiro-Lopatinskii condition [7]. ■ Proposition 2.4. If boundary-value problem [START_REF] Bars | Thermo-elliptical instability in a rotating cylindrical shell[END_REF] is elliptic in the interval [ω 1 , ω 2 ], then its spectrum is discrete (or empty) in this interval. Otherwise, the spectrum is essential in [ω 1 , ω 2 ].

Proof. If the boundary-value pressure problem is elliptic, the operator P ω admits a parametrix L (i.e. an approximate inverse) such that P ω • L = Id + R, where R is a pseudo-differential operator of degree -1. Hence, R is a compact operator. Moreover, P ω and R are also analytic functions. Then, the result follows from the Fredholm analytic theorem (see Theorem VI.14 in [9]). Otherwise, the values ω ∈ [ω 1 , ω 2 ] belong to the essential spectrum [10]. ■ BOX 2.1: Normal vector induced by the metric on ∂V

The normal vector n defined from the metric on ∂V is computed as follows. The principal symbol p defines a metric given by g ⋆ := (g ⋆ ) ij k i k j for the covectors k ∈ R 3 in the cotangent space, whose inverse g := (g ⋆ ) -1 induces a metric in V defined by g := g ij dx i dx j . Then, the normal vector n ⋆ computed with respect to g ⋆ is given by

n ⋆ ∝ g ⋆ n, g ij n ⋆ i n ⋆ j = 1,
where n = (n 1 , n 2 , n 3 ) ⊤ is the unit vector normal to ∂V in Cartesian coordinates with n 2 1 + n 2 2 + n 2 3 = 1. The constant of proportionality is such that n ⋆ has a unit norm with respect to the metric g. The vector n ⋆ is also called the conormal vector (with respect to g ⋆ ) in the study of partial differential equations [8]. Note that g ⋆ is a contravariant metric tensor associated with the operator (g ⋆ ) ij ∂ x i ∂ x j . Thus, the conormal n ⋆ is the contravariant description (in the space with the metric defined by g) of the covariant normal vector n to ∂V (in the space with the Euclidean metric).

In practice, the difficult step is to evaluate the principal symbol of the Dirichlet-to-Neumann operator (Proposition 2.2). In practice, this can be done as follows. The principal symbol p induces a metric in V given by g := g ij dx i dx j . Let us assume that x i = x i (q j ) with the curvilinear coordinates (q 1 , q 2 , q 3 ) such that the boundary ∂V is given by q 1 = constant (for instance). Then, we can rewrite the two-dimensional metric restricted on ∂V as g| ∂V = g ij dq i dq j . Then, the dual of the metric induced by g on ∂V is given by g ⋆ | ∂V = ( g ⋆ ) ij k i k j with g ⋆ = g -1 .

Back to the pressure equation

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

g ⋆ =   A F E F B D E D C   , g = 1 det(g ⋆ )   BC -D 2 -CF + DE -BE + DF -CF + DE AC -E 2 -AD + EF -BE + DF -AD + EF AB -F 2   , (23a,b)
with det(g ⋆ ) = ABC -AD 2 -BE 2 -CF 2 + 2DEF . For a given normal vector n = (n 1 , n 2 , n 3 ) ⊤ in Cartesian coordinates, the conormal vector n ⋆ is thus given by

n ⋆ = 1 α   n 1   A F E   + n 2   F B D   + n 3   E D C     , (24)
with α = [An 2 1 + Bn 2 2 + Cn 2 3 + 2(Dn 2 n 3 + En 1 n 3 + F n 1 n 2 )] 1/2 such that g ij n ⋆ i n ⋆ j = 1. Thanks to expression [START_REF] Cartan | Sur les petites oscillations d'une masse fluide[END_REF], we clearly see that the conormal naturally appears in equations (13a-c). Thus, we can apply Proposition 2.3 to determine the nature of the boundary-value pressure problem.

Let us consider here the aligned case with Ω = (0, 0, Ω z ) ⊤ , for which the equations are then greatly simplified (see Box 1.2). In this case, the two metrics reduce to

g ⋆ = k 2
x + k 2 y + Ck 2 z , g = dx 2 + dy 2 + C -1 dz 2 , (25a,b)

Figure 1 .

 1 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 Ω.

Figure 2 .

 2 Figure 2. Sketch of the mathematical model for geophysical vortices (not to scale). (a) Meridional cross-section of a vortex centred at colatitude θ, which is embedded within an ocean or an atmosphere (blue region). (b) Front view of a flattened ellipsoidal vortex of semi-axes a ≫ b ≫ c, which is stratified under the uniform planet's gravity g (i.e. the background density ρ 0 decreases with increasing z).

  royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 .

 3 Figure 3. Sparse upper triangular block matrices An (left) and Bn (right) for n = 5 and 1 Ω = 1z.

Proposition 3 . 3 .

 33 The number of non-zero eigenvalues of QEP (3.2) is n(n + 1)(n + 5)/3.

  royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Elliptic in V Hyperbolic in V Low-frequency (surface) modes Inertial-Gravity Modes (IGMs)

Figure 4 .

 4 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 |ω| ≤ ω + (Proposition 4.1). The velocity equation is either hyperbolic or elliptic in volume. The spectrum is spanned by IGMs and low-frequency (surface) modes.

Figure 5 . 1 (

 51 Figure 5. Spectrum σ 2 in the aligned case 1 Ω = 1z. (a) Eigenvalues ω 3,5 ≥ 0 in σ 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 σ 2 in which the problem is hyperbolic. (b) Number of eigenvalues in σ 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/Ωs.

  royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 6 .

 6 Figure 6. Spatial structure of the low-frequency (surface) modes for u ∈ V 0 10 in the aligned case 1 Ω = 1z as a function of N/Ω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 |ζ| (in logarithmic scale) in the (Oxz) plane (with rescaled axes to have a 1:1 ratio).

  royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 7 .

 7 Figure 7. Spatial structure of IGMs for u ∈ V 0 10 in the aligned case 1 Ω = 1z as a function of N/Ω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 |ζ| (in logarithmic scale) in the (Oxz) plane (with rescaled axes to have a 1:1 ratio).

  royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 8 .

 8 Figure 8. Density of IGMs in σ 1 . Comparison between asymptotic density dπ∞/dω given in Proposition 4.2 and σ 2,35 obtained from polynomial computations. Left: Flattened ellipsoid b/a = 1 and c/a = 0.1 in the aligned case 1 Ω = 1z with N/Ωs = 5. Right: Sphere b/a = c/a = 1 in the perpendicular case 1 Ω • 1z = 0 with N/Ωs = 10.

1 c/a = 10 Figure 9 .

 1109 Figure 9. Asymptotic density dπ∞/d ω of IGMs with normalised eigenvalues ω := ω/ω + , given by Proposition 4.2. Left: Sphere b/a = c/a = 1 in the aligned case 1 Ω = Right: Axisymmetric ellipsoids with b/a = 1 and varying c/a ratios. Top panel computed for N/Ωs = 0.5, and bottom one for N/Ωs = 5.

  royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table

  royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .pure inertial mode (ω, u) satisfies in an ellipsoid the intriguing integral property[START_REF] Vantieghem | Inertial modes in a rotating triaxial ellipsoid[END_REF][START_REF] Ivers | Enumeration, orthogonality and completeness of the incompressible Coriolis modes in a tri-axial ellipsoid[END_REF]

Proposition 2 . 1 .

 21 If Ψ : ∂V → C is a field defined on the boundary, then the Dirichlet-to-Neumann operator is defined by N (Ψ) := ∂ n ⋆ Φ where ∂ n ⋆ = n ⋆ • ∇ is the normal derivative with respect to the vector n ⋆ , and where Φ is the solution of the Dirichlet problem in V given by P ω (Φ) = 0 in V and Φ| ∂V = Ψ on ∂V .

  equations, it is worth considering a master equation that governs the evolution of either u or Φ. If we multiply Euler equation (2.3a) by iω and replace iωζ by buoyancy equation (2.3c), we obtain for ω ̸ = 0

The wave properties are deeply tied to the nature of the mathematical problem. Instead of solving the royalsocietypublishing.org/journal/rspa Proc. R. Soc. A . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . primitive

tangent is here understood in the language of differential geometry[4]. For instance, if ∂V is described by q 1 = constant for a set of curvilinear coordinates (q 1 , q

, q

), then V only involves derivatives (∂ q2 , ∂ q3 ) on the tangent plane at a given point.

Acknowledgements. JV would like to dedicate this work to the memory of Norman R. Lebovitz (1935Lebovitz ( -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.

Electronic Supplementary Material of "Inertia-gravity waves in geophysical vortices"

Jérémie Vidal * ID & Yves Colin de Verdière † February 14, 2024

Abstract

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

Boundary-value pressure problem

We consider a diffusionless and incompressible fluid, which is enclosed in a smooth bounded domain V (whose boundary is denoted by ∂V ). The fluid is rotating at the angular velocity Ω = (Ω x , Ω y , Ω z ) ⊤ with respect to an inertial frame, and is stratified in density under the imposed gravity field g = -g(z)1 z within the Boussinesq approximation. At rest, the pressure P 0 and the density satisfy the hydrostatic balance ∇P 0 = [ρ * + ρ 0 (z)]g, where ρ * is the typical (mean) density of the fluid within the Boussinesq approximation. Then, we seek small-amplitude oscillatory perturbations (where ω ∈ R is the angular frequency) for the velocity u exp(iωt), the pressure Φ exp(iωt), and the density ζ exp(iωt) upon the hydrostatic basic state. These perturbations are governed in the rotating frame by the incompressible Euler equations

and the buoyancy equation

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

(3) * Université Grenoble Alpes, CNRS, ISTerre, 38000 Grenoble, France † Université Grenoble Alpes, CNRS, Institut Fourier, 38000 Grenoble, France

1: Poincaré problem

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

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

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

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

Finally, equation ( 3) is supplemented with the no-penetration boundary condition (BC) given by u • n| ∂V = 0 (i.e. u is tangent to ∂V ), where n is the unit vector normal to the boundary.

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

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

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

Second-order equation in the volume

To find the equation satisfied by the pressure Φ in V , we first find the equation relating u and Φ (Proposition 1.1). Then, we enforce the incompressibility condition to obtain the pressure equation (Proposition 1.2). Finally, we give in Proposition 1.3 the pressure equation when N is constant.

It will be of prime importance for the microlocal analysis presented in Section 2.

Proposition 1.1. The velocity u is related to the pressure Φ by the equation

and with

BOX 1.2: Aligned case

When Ω = Ω z 1 z and N is constant, pressure equation ( 10) simplifies into

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

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

which simplifies into

Boundary condition

Proposition 1.4. The pressure BC is given by f • n| ∂V = 0, where the components of f are given in equations (13a-c).

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

Since we have

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

1.3 First-order equation when D ω (z) = 0

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

and the conormal vector is given by n

, where the Euclidian normal vector is n = (n 1 , n 2 , n 3 ) ⊤ with n 2 1 + n 2 2 + n 2 3 = 1. To demonstrate that the pressure problem is not elliptic for a given frequency 0 < |ω| < ω -, it is sufficient to find one location on ∂V where Proposition 2.3 is not valid. We look at the behaviour at the equator, where the tangent plane is vertical with n = 1 x (for instance) such that the pressure BC is simply f x = 0. The restriction of the metrics on ∂V gives

and the conormal vector is simply given by n ⋆ = 1 x . The principal symbol of the pressure BC is then given by

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

Boundary-layer analysis

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

with the (dimensionless) Ekman number E = ν/(Ω s L 2 ) and E κ = E/P r where P r = ν/κ is the Prandtl number, and with the normalised BV frequency N = N/Ω s . Multiplying equation (29a) by u † and the complex-conjugate of equation (29b) by N 2 ζ, we obtain after volume integration

where the scalar product between two scalar quantities is defined as [START_REF] Veronis | The analogy between rotating and stratified fluids[END_REF] is also valid in a sphere with a radial gravity [11]. Finally, the problem is supplemented with BC. We assume that the velocity field satisfies the stress-free BCs

where ϵ(u) = (1/2) ∇u + (∇u) ⊤ is the strain-rate tensor. SF-BCs (31a,b) allow the tangential velocity to slip on ∂V . For the density perturbation, we enforce a Neumann BC ∇ζ • n| ∂V = 0. A density jump is thus possible between the vortex and the ambient fluid (as considered in idealised models [12,13,14]), but without exchange of mass between the two fluids. We can simplify equation ( 31) using BLT when E ≪ 1. Formally, classical BLT for rotating fluids is not expected to be strongly modified by stratification when N is not too large (e.g. for rigid boundaries [15]). Moreover, it is known that the boundary-layer flow is E 1/2 smaller than the bulk flow for SF-BCs [16,17]. This considerably simplifies the BLT, because explicit expressions for the boundary-layer solutions are no longer required to estimate the decay rate of the modes. We seek the variables at the leading order in E for the SF-BCs as

where [ω 0 , u 0 ] is the eigenvalue-eigenvector pair of a diffusionless IGM, λ 1 ∈ C is the first-order correction of the eigenvalue, and [u 1 , ζ 1 ] are the first-order corrections within the boundary layer such that u 0 + E 1/2 u 1 satisfies SF-BCs (31a,b) and ζ 0 + E 1/2 ζ 1 satisfies the Neumann BC for the density. We substitute the above asymptotic expansions into equation ( 30) and we obtain

at the order E, where we have used the SF-BCs and the Neumann BC to simplify the volume integrals on the right-hand side. Since the latter is real-valued and negative, there is no frequency correction at the leading order in E due to SF-BCs such that ℑ m (λ 1 ) = 0 and ℜ e (λ 1 ) = τ 1 ≤ 0. We recover from equation ( 33) the viscous decay rate of pure inertial modes with SF-BCs [17]. Finally, we can crudely estimate whether forced IGMs would quickly decay or not in largescale stratified vortices by plugging geophysical estimates into equation [START_REF] Brunet | Linear and nonlinear regimes of an inertial wave attractor[END_REF]. Typical values for the kinematic viscosity are ν ∼ 10 -6 m 2 .s -1 for water and ν ∼ 4 × 10 -7 m 2 .s -1 for gas giants (according to ab-initio simulations [18]). Using Table 1 in the main text, typical values for the Ekman number are thus E ∼ 10 -8 -10 -12 for Mediterranean eddies and E ∼ 10 -12 -10 -17 for Jovian vortices. Estimating the Prandtl number depends on whether the stratification is due to thermal effects or compositional ones (for which P r ≫ 1 in both water and gas giants). Heat diffusion in Mediterranean eddies is characterised by P r ≃ 0.7, and typical values P r ∼ 10 -2 -1 are expected for thermal diffusion in gas giants [18]. Hence, we have N E κ ∼ 10 -4 -10 -8 for Mediterranean eddies and N E κ ∼ 10 -6 -10 -11 for Jovian vortices. Consequently, the largest-scale IGMs are not expected to be strongly damped by diffusion in stratified vortices.