Jérémie Vidal

Sylvie Su

David Cébron

Compressible fluid modes in rigid ellipsoids: towards modal acoustic velocimetry

Keywords: waves in rotating fluids, acoustics, geophysical and geological flows

Motivated by planetary-driven applications and experiments in non-spherical geometries, we study compressible fluid modes in rotating rigid ellipsoids. Such modes are also required for modal acoustic velocimetry (MAV), a promising non-invasive method to track the velocity field components in laboratory experiments. To calculate them, we develop a general spectral method in rigid triaxial ellipsoids. The description is based on an expansion onto global polynomial vector elements, satisfying the non-penetration condition on the boundary. Then, we investigate the diffusionless compressible modes in rotating (and magnetised) rigid ellipsoids. The spectral description is successfully benchmarked against three-dimensional finite-element computations and analytical predictions. A spectral convergence is obtained. Our results have direct implications for MAV in experiments, for instance in the ZoRo experiment (gas-filled rigid spheroid). So far, deformation and rotational effects have been theoretically considered separately, as small perturbations of the solutions in non-rotating spheres. We carefully compare the perturbation approach, in this illustrative geometry, to the polynomial solutions. We show that second-order ellipticity effects are often present, even in weakly deformed ellipsoids. Moreover, high-order effects due to rotation and/or ellipticity should be observed for some acoustic modes in experimental conditions. Thus, perturbation theory should be used with care in MAV. Instead, the spectral polynomial method paves the way for future MAV applications in fluid experiments with rigid ellipsoids.

Introduction

Planetary context

The dynamics of planetary fluid layers often defies our physical knowledge. A fundamental understanding can be obtained with massive numerical simulations (e.g. [START_REF] Schaeffer | Turbulent geodynamo simulations: a leap towards Earth's core[END_REF]. The success of this approach is somewhat surprising, considering the gap between the parameters of the simulations and the expected planetary values. Realistic numerical simulations will remain for a long time inaccessible to the computational power. Moreover, a long-term endeavour is to go beyond the spherical geometry in these models. Indeed, planetary fluid bodies are closer to ellipsoids than to spheres, for instance due to tidal effects (e.g. [START_REF] Chandrasekhar | Ellipsoidal Figures of Equilibrium[END_REF][START_REF] Nduka | The Roche problem in an eccentric orbit[END_REF]). However, simulations in deformed spheres are very challenging. Numerical methods usually used in non-spherical geometries, such as finite elements (e.g. [START_REF] Cébron | Elliptical instability in hot Jupiter systems[END_REF], are not very efficient for this problem. Spectral numerical methods have been developed in (weakly) non-spherical containers (e.g. [START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF][START_REF] Rieutord | An algorithm for computing the 2D structure of fast rotating stars[END_REF]), but they rely on non-orthogonal spherical-like coordinates [START_REF] Bonazzola | Numerical approach for high precision 3D relativistic star models[END_REF]) that strongly complicate the numerics. Furthermore, accounting for density variations (compressibility) is another computational burden. Thus, mainly incompressible flows have been simulated in deformed geometries (e.g. [START_REF] Cébron | A systematic numerical study of the tidal instability in a rotating triaxial ellipsoid[END_REF][START_REF] Favier | Generation and maintenance of bulk turbulence by libration-driven elliptical instability[END_REF].

To gain physical insight, we can develop reduced global models in ellipsoids. Reduced models have proven useful to approach more realistic values in spherical convection (e.g. [START_REF] Guervilly | Turbulent convective length scale in planetary cores[END_REF]. In ellipsoids, a first step would be to investigate the global modes of rotating compressible fluids. Indeed, the inertial modes (sustained by global rotation) play a fundamental role in rotating flows (e.g. [START_REF] Zhang | Theory and Modelling of Rotating Fluids: Convection, Inertial Waves and Precession[END_REF], for incompressible flows). They can be triggered by orbital forcings (Le [START_REF] Bars | Flows driven by libration, precession, and tides[END_REF] and coupled nonlinearly to yield flow instabilities (e.g. [START_REF] Kerswell | The instability of precessing flow[END_REF][START_REF] Kerswell | Elliptical instability[END_REF][START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF]. Thus, reduced compressible models in non-spherical rotating domains are worth developing.

Developing such models was largely inhibited by the mathematical complexity of the ellipsoidal system (e.g. [START_REF] Cartan | Sur les petites oscillations d'une masse de fluide[END_REF]. Analytical expressions for the acoustic modes do exist in non-rotating spheroids [START_REF] Chang | Natural resonant frequency of a prolate acoustical resonator[END_REF][START_REF] Chang | Natural resonant frequencies of an oblate acoustical resonator[END_REF] and ellipsoids [START_REF] Willatzen | Eigenmodes of triaxial ellipsoidal acoustical cavities with mixed boundary conditions[END_REF]. Similarly, analytical solutions of the incompressible inertial modes [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF] exist in spheroids [START_REF] Zhang | Theory and Modelling of Rotating Fluids: Convection, Inertial Waves and Precession[END_REF], but not in ellipsoids (except for a few, see [START_REF] Vantieghem | Inertial modes in a rotating triaxial ellipsoid[END_REF]. Unfortunately, these analytical solutions cannot be extended to account for additional ingredients, notably (i) global rotation for the acoustic modes and (ii) compressibility for the inertial modes. In the latter case, extensions towards compressible models have been attempted. Yet, they were restricted to neutrally buoyant fluids, by considering specific isentropic (polytropic) states within the anelastic approximation (Clausen & Tilgner 2014, in ellipsoids). A consistent mathematical treatment of the acoustic and inertial modes may appear as a prerequisite, to pave the way for future planetary-driven reduced models.

Experimental context

A complementary physical understanding can be obtained with experimental analogues. Despite they lack some ingredients compared to simulations, experiments can often probe more turbulent regimes (as measured by the Reynolds number). Motivated by planetary-driven applications, experiments have been conducted in spherical geometries (e.g. [START_REF] Kelley | Inertial waves driven by differential rotation in a planetary geometry[END_REF][START_REF] Triana | Helioseismology in a bottle: modal acoustic velocimetry[END_REF][START_REF] Tigrine | Torsional Alfvén waves in a dipolar magnetic field: experiments and simulations[END_REF]. Similarly, planetarydriven ellipsoidal experiments have been built to mimic tidal deformations, for instance with water (e.g. [START_REF] Noir | Experimental evidence of inertial waves in a precessing spheroidal cavity[END_REF][START_REF] Noir | Experimental study of libration-driven zonal flows in non-axisymmetric containers[END_REF][START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF][START_REF] Lemasquerier | Libration-driven flows in ellipsoidal shells[END_REF][START_REF] Le Reun | Experimental study of the nonlinear saturation of the elliptical instability: inertial wave turbulence versus geostrophic turbulence[END_REF] or the gas-filled spheroidal experiment ZoRo [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]. A major difficulty in the experimental approach is to reconstruct, from a sparse set of measurements, the velocity field within the fluid domain.

Conventional velocimetry methods, such as particle-imaging techniques, give often partial information on the flow components. In order to sample larger portions of the fluid domain, a large and flexible imaging apparatus is required. This is quite challenging to develop for (rapidly) rotating experiments. The fluid must be also seeded with neutrally buoyant tracer particles, acting as scatterers (of light or sound). Yet, buoyant tracers can float, sink or be quickly centrifuged (especially in rapidly rotating gas-filled experiments). Moreover, some methods work only in non-opaque fluids (for the part of the light spectrum probed by the sensor), such as particle image velocimetry and laser Doppler velocimetry, and alternatives to optical techniques must be employed for liquid metal experiments (e.g. [START_REF] Kelley | Inertial waves driven by differential rotation in a planetary geometry[END_REF][START_REF] Tigrine | Torsional Alfvén waves in a dipolar magnetic field: experiments and simulations[END_REF]. Consequently, it is worth developing flexible imaging techniques for transparent or opaque fluids, with a large sampling domain and using non-invasive sensors.

Modal acoustic velocimetry (MAV) does match these requirements. This technique consists in the identification of the normal modes of oscillations of the fluid. Since the acoustic modes are highly sensitive to the flow properties, the observed distribution of acoustic resonant frequencies can feed an inverse problem, to gain a passive access to the (hidden) dynamics. This technique is routinely used in helioseismology and asteroseismology (e.g. [START_REF] Aerts | Asteroseismology[END_REF], to remotely probe the interior conditions of stars. Thus, there are strong motivations to use MAV as a complementary imaging technique for rotating fluid experiments [START_REF] Triana | Helioseismology in a bottle: modal acoustic velocimetry[END_REF][START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]. It is capable of imaging remotely the entire medium on relatively short time scales (compared to the flow time scale), since the global modes are sensitive to the properties of the whole fluid domain. Moreover, only non-invasive sensors are required on the boundary of the container. Consequently, MAV is a very promising technique. This would allow the experimentalists to recover the three components of the (rather large-scale) flow, simply by using non-invasive acoustic probes.

Modal acoustic velocimetry

MAV aims at providing flow reconstructions, by measuring an experimental acoustic signal. Synthetic predictions of the acoustic signal are thus required, to determine the flow that reproduces at best the experimental data. A high accuracy on the synthetic solutions is desirable, to restrict the survey of the parameter space. Numerical computations with local numerical methods (such as finite elements) could be used to determine the (visco-thermal) acoustic resonances in any bounded geometry (e.g. [START_REF] Berggren | Acoustic boundary layers as boundary conditions[END_REF][START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]). Yet, the computations are slow and often limited to a few solutions. Thus, they cannot be efficiently coupled with inverse schemes (so far), thereby limiting their practical use for MAV. Alternatively, the analytical theory of the acoustic modes has been used, since analytic solutions are available in various idealised geometries. For instance, [START_REF] Triana | Helioseismology in a bottle: modal acoustic velocimetry[END_REF] used the spherical modes for MAV, to retrieve the mean azimuthal velocity in a rotating spherical experiment. Despite this proof-of-concept validation, the identification of the spherical acoustic modes in the experimental data is often difficult, due to the azimuthal degeneracy of the spherical modes. Then, a misleading physical identification would strongly affect the reconstruction of the velocity field. For this reason, [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF] excited the acoustic modes of a spheroid in the ZoRo experiment, in order to (partially) lift the azimuthal degeneracy of the modes. For this experimental set-up, the ellipsoidal acoustic modes should be employed.

So far, only small departures from the standard spherical solutions of the acoustic equation were considered with perturbation theory (to overcome the difficulty of the ellipsoidal coordinates). Indeed, the angular eigenfrequency of a global mode is shifted by a small amount when the cavity is squashed (e.g. [START_REF] Mehl | Acoustic eigenvalues of a quasispherical resonator: second order shape perturbation theory for arbitrary modes[END_REF][START_REF] Guianvarch | Acoustic field in a quasispherical resonator: unified perturbation model[END_REF]. Similarly, frequency shifts are expected due to global rotation (e.g. [START_REF] Dahlen | Theoretical Global Seismology[END_REF], the velocity field (e.g. [START_REF] Aerts | Asteroseismology[END_REF]) and diffusive effects [START_REF] Moldover | Gas-filled spherical resonators: Theory and experiment[END_REF]. Previous MAV applications [START_REF] Triana | Helioseismology in a bottle: modal acoustic velocimetry[END_REF][START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]) combined separate perturbation theories, one for each effect, which have different mathematical expressions. A similar perturbation framework has proven accurate enough in metrology, in which either (i) non-rotating and slightly deformed spherical resonators (e.g. [START_REF] Moldover | Gas-filled spherical resonators: Theory and experiment[END_REF][START_REF] Guianvarch | Acoustic field in a quasispherical resonator: unified perturbation model[END_REF] or (ii) rotating gyroscopes [START_REF] Bruneau | A rate gyro based on acoustic mode coupling[END_REF][START_REF] Ecotiere | Inertial coupling of resonant normal modes in rotating cavities: Acoustic gyrometers for high rotation rates[END_REF] were considered. Undertaking such a (difficult) task may not be satisfactory in moderately deformed ellipsoids.

The perturbation approach must be carefully assessed against consistent solutions in the experimental conditions. Indeed, possible cross-effects have been hitherto neglected, such as the interplay between rotation and ellipticity. In typical experimental conditions, the effective frequency resolution for MAV is 10 -1 -1 Hz (for typical observational temporal windows between 10 and 100 seconds, see [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]. For most of the acoustic modes that can be excited (with frequencies up to several kHz), even slight modelling errors could be misleading for the velocity reconstruction. Therefore, perturbation theory should be used with care (without any prior validation).

Motivations

In the present study, we focus on the theory and computations of the rotating compressible modes. To simplify the physical problem, we consider full ellipsoids (no inner core) with rigid boundaries, on which the flows satisfy the non-penetration condition. This is a reasonable starting assumption for planetary-driven (reduced) models, which also matches the geometry of many fluid experiments. Firstly, we aim at developing a new theoretical description of the diffusionless modes in compressible and rotating rigid ellipsoids. The diffusionless approximation filters out the thin boundary layers, which are (often) unimportant for the bulk dynamics [START_REF] Berggren | Acoustic boundary layers as boundary conditions[END_REF]. This assumption also greatly simplifies the mathematical analysis. In ellipsoids, the acoustic [START_REF] Chang | Natural resonant frequency of a prolate acoustical resonator[END_REF][START_REF] Chang | Natural resonant frequencies of an oblate acoustical resonator[END_REF][START_REF] Willatzen | Eigenmodes of triaxial ellipsoidal acoustical cavities with mixed boundary conditions[END_REF] and inertial [START_REF] Vantieghem | Inertial modes in a rotating triaxial ellipsoid[END_REF][START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF] modes are known to be smooth and differentiable. Thus, we can develop polynomial approximations of the compressible modes in Cartesian coordinates. This will allow us to overcome the mathematical complexity of both ellipsoidal coordinates (e.g. [START_REF] Cartan | Sur les petites oscillations d'une masse de fluide[END_REF]) and non-orthogonal coordinates (e.g. [START_REF] Bonazzola | Numerical approach for high precision 3D relativistic star models[END_REF][START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF]. Such a path has been followed in the pioneering work of [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF]. Yet, he only considered ellipsoids with free-surface boundaries, which are not relevant for experiments.

Secondly, we will revisit MAV in rigid ellipsoids. Previous works have mainly employed perturbation theory [START_REF] Triana | Helioseismology in a bottle: modal acoustic velocimetry[END_REF][START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF], to deal with small departures from non-rotating spherical containers. Perturbation calculations are often valid for most planetary [START_REF] Dahlen | Theoretical Global Seismology[END_REF]) and stellar applications [START_REF] Lignières | Acoustic oscillations of rapidly rotating polytropic stars-I. Effects of the centrifugal distortion[END_REF][START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF]). Yet, the validity of perturbation theory remains elusive for rotating experiments in (moderately) deformed spheres (e.g. [START_REF] Noir | Experimental evidence of inertial waves in a precessing spheroidal cavity[END_REF][START_REF] Noir | Experimental study of libration-driven zonal flows in non-axisymmetric containers[END_REF][START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF][START_REF] Lemasquerier | Libration-driven flows in ellipsoidal shells[END_REF][START_REF] Le Reun | Experimental study of the nonlinear saturation of the elliptical instability: inertial wave turbulence versus geostrophic turbulence[END_REF]. This currently prevents us from using perturbation theory for MAV in rigid ellipsoids.

To sum up, the present paper has a interdisciplinary twofold purpose. First, we aim at introducing a new polynomial spectral method in rigid ellipsoids. Second, we will apply it to solve the magneto-acoustic problem in rotating rigid ellipsoids (as a proof of concept). The paper is organised as follows. We describe the full physical model in §2. Then, we introduce in §3 the new polynomial spectral method, valid for any vector field satisfying the non-penetration condition in rigid ellipsoids. Next, results for the magnetoacoustic modes are presented and validated in §4 and §5. Finally, we end the paper with a discussion and outline some perspectives in §6.

Description of the model

Assumptions

Taking into account all the physical ingredients in the mathematical model is unnecessary for MAV. We build a forward model accounting only for the key physical ingredients, namely adiabatic compressibility, rotation and the leading-order magnetic field component. As depicted in figure 1, we consider a fluid-filled rigid ellipsoid, of semi-

* = Ω * s 1 Ω
, where 1 Ω is a unit vector with respect to an inertial frame. In the following, we work exclusively in the rotating frame where the boundary is steady and employ the Cartesian coordinates (x * , y * , z * ). In the rotating frame, the ellipsoidal boundary ∂V is expressed by (x * /a * ) 2 +(y * /b * ) 2 +(z * /c * ) 2 = 1. Moreover, the fluid can have a possible (steady) background mean flow U * 0 (r * ) (such as a differential rotation with respect to the solid-body rotation of the container, accompanied with a meridional circulation due to viscous effects), with r * = (x * , y * , z * ) ⊤ the position vector. This mean flow is sub-sonic for experimental conditions (yielding ∇ • U * 0 = 0). For the sake of the numerical validation, we assume that the fluid has a spatially uniform dynamic (shear) viscosity η, bulk viscosity η B and magnetic diffusivity ν m . In the theory, diffusive effects will be entirely neglected since they are often significant only in the boundary layers [START_REF] Berggren | Acoustic boundary layers as boundary conditions[END_REF]. We also discard non-adiabatic effects to consider isentropic fluids, characterised in the reference state by the homogeneous background density ρ * 0 . Within our idealised framework, we neglect gravitational effects, as well as the dynamical pressure generated by the background velocity. They are small in rotating laboratory experiments (when the flow rotation is small compared to the speed of sound). For the same reason, we also neglect centrifugal effects which are likely negligible for moderate rotation (although they could be included within our framework). Hence, the reference state has the background pressure P * 0 and the speed of sound C * 0 = K * 0 /ρ * 0 , with K * 0 the isentropic bulk modulus, that are spatially uniform. For an isentropic gas, we have K * 0 = γP * 0 with γ the adiabatic index. The fluid is also pervaded by a background magnetic field B * 0 (r * ). The leading-order spatial component of the background magnetic field is usually the (aligned) dipole in fluid experiments, such as in the Maryland spherical Couette flow experiment (e.g. [START_REF] Kelley | Inertial waves driven by differential rotation in a planetary geometry[END_REF] or the Derviche Tourneur Sodium (DTS) experiment (e.g. [START_REF] Tigrine | Torsional Alfvén waves in a dipolar magnetic field: experiments and simulations[END_REF]. Thus, we only retain for the background magnetic field its spatially uniform component, that is B * 0 = B * 0 1 B with B * 0 the strength of the magnetic field and 1 B its unit direction.

We expand the velocity v * , the density ρ * , the pressure P * and the magnetic field B * as isentropic perturbations around the (steady) isentropic background state. This reads

v * (r * , t * ) = U * 0 (r * ) + u * 1 (r * , t * ), (2.1a) [ρ * , P * , B * ] (r * , t * ) = [ρ * 0 , P * 0 , B * 0 ] + [ρ * 1 , p * 1 , b * 1 ] (r * , t * ), (2.1b) ∇ • U * 0 = ∇ • B * 0 = ∇ • b * 1 = 0. (2.1c)
Then, linearising the governing equations around the background state leads to the governing (dimensional) equations

MΩ 8 × 10 -3 10 -2 2 × 10 -1 ≪ 10 -2 MA 2 × 10 -4 9 × 10 -5 -4 × 10 -3 0 0 Le < ∞ 0.01 -0.3 0 0 |Ro| 0.03 -400 0.1 0.1 10 -1 Ek > 10 -8 > 10 -7 10 -6 > 10 -6
∂u * 1 ∂t * + 2 Ω * × u * 1 = -∇(p * 1 /ρ * 0 ) + 1 ρ * 0 µ 0 (∇ × b * 1 ) × B * 0 -(U * 0 • ∇) u * 1 -(u * 1 • ∇) U * 0 + ν f * v (u * 1 ), (2.2a) ∂b * 1 ∂t * -∇ × (U * 0 × b * 1 ) = ∇ × (u * 1 × B * 0 ) + ν m ∇ 2 b * 1 , (2.2b) ∂ρ * 1 ∂t * + U * 0 • ∇ρ * 1 = -ρ * 0 ∇ • u * 1 , (2.2c)
with µ 0 the magnetic permeability of the vacuum, ν = η/ρ * 0 the kinematic viscosity, the viscous force

f * v (u * 1 ) = ∇ 2 u * 1 + 1 3 + ν B ν ∇(∇ • u * 1 ) (2.3)
and the bulk kinematic viscosity ν B = η B /ρ * 0 . The governing equation for the pressure is obtained by using the equation of state for an isentropic fluid. We get

∂p * 1 ∂t * + U * 0 • ∇p * 1 = -ρ * 0 C * 0 2 ∇ • u * 1 .
(2.4) Finally, equations (2.2) are supplemented with boundary conditions. The velocity satisfies the non-penetration (Dirichlet) condition on the rigid ellipsoidal wall u * 1 • 1 n = 0, where 1 n is the unit vector normal to the boundary. In addition to the non-penetration condition, diffusive and pressure boundary conditions will be enforced in the diffusive numerical computations that will be used for the validation (see appendix A). However, no additional boundary condition is required in the diffusionless theory.

Dimensionless variables

For the numerical convenience, we work in dimensionless variables. We use the semimajor axis a * as length scale, a * /C * 0 as time scale, ρ * 0 as density scale and ρ * 0 C * 0 2 as pressure scale. We choose the typical amplitude B * 0 of the dimensional background magnetic field as magnetic scale. We also introduce the typical amplitude of the back-ground velocity field U * 0 . The dimensionless fields are written without an asterisk, to distinguish them from their dimensional counterparts. Then, equations (2.2)-(2.4) read in dimensionless form

∂u 1 ∂t + 2M Ω 1 Ω × u 1 = -∇p 1 + M 2 A (∇ × b 1 ) × 1 B (2.5a) -M S [(U 0 • ∇) u 1 + (u 1 • ∇) U 0 ] + Re -1 S f v (u 1 ), ∂p 1 ∂t + M S (U 0 • ∇p 1 ) = -∇ • u 1 , (2.5b) ∂b 1 ∂t -M S ∇ × (U 0 × b 1 ) = ∇ × (u 1 × 1 B ) + Rm -1 S ∇ 2 b 1 . (2.5c)
We have introduced in equations (2.5) the sonic, rotational and the Alfvénic Mach dimensionless numbers

M S = U * 0 C * 0 , M Ω = a * Ω * s C * 0 and M A = V * A C * 0 , (2.6a-c)
where V * A = B * 0 / ρ * 0 µ 0 is the Alfvén velocity. They compare the typical time scales for the flow, rotation and the magnetic field with the sonic time scale. Note that M S can be negative, if the background flow is retrograde compared to global rotation. We can also define (for rotating flows only) the Rossby number Ro = M S /M Ω , which measures the strength of the background mean flow compared to global rotation. The diffusive effects are quantified by the sonic Reynolds and magnetic Reynolds numbers (based on the speed of sound), defined as

Re S = C * 0 a * ν and Rm S = C * 0 a * ν m . (2.7a,b)
In the presence of global rotation, we can also introduce the Ekman number Ek = ν/(Ω * s a * 2 ). Finally, M A is related in rotating magnetohydrodynamics to the Lehnert number Le = M A /M Ω . The typical values of these numbers in experiments are given in table 1. Their magnitudes will allow us to simplify the theory of the compressible modes. Finally, given the high accuracy of typical acoustic measurements (e.g. [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF], note that we must be able to resolve frequency variations as small as 10 -4 -10 -3 in dimensionless units.

Master wave-like equation

We are now in a position to build the wave-like equation of the rotating magnetoacoustic modes in rigid ellipsoids. We leave aside diffusive and non-adiabatic effects. Indeed, (rapidly) rotating fluid experiments are usually characterised by small diffusive effects (except possibly in the thin boundary layers, see [START_REF] Berggren | Acoustic boundary layers as boundary conditions[END_REF], that is Re S ≫ 1 (or Ek ≪ 1 with rotation) in table 1. They can (often) be considered as small perturbations of the diffusionless modes (e.g. [START_REF] Moldover | Gas-filled spherical resonators: Theory and experiment[END_REF]. We can recast diffusionless primitive equations (2.5) into a wave-like equation for the fluid particle displacement vector ξ 1 (r, t). This is be a prerequisite to obtain the frequency shift of the modes generated by the background mean flow (e.g. [START_REF] Aerts | Asteroseismology[END_REF]. The Eulerian velocity perturbation is related to ξ 1 by (e.g. [START_REF] Lynden-Bell | On the stability of differentially rotating bodies[END_REF])

u 1 = ∂ξ 1 ∂t + (U 0 • ∇) ξ 1 -(ξ 1 • ∇) U 0 .
(2.8)

We also obtain from pressure and induction equations in dimensionless units (e.g. [START_REF] Bernstein | An energy principle for hydromagnetic stability problems[END_REF][START_REF] Frieman | On hydromagnetic stability of stationary equilibria[END_REF])

p 1 = -∇ • ξ 1 , b 1 = ∇ × (ξ 1 × 1 B ). (2.9a,b)
Now, we substitute (2.8)-(2.9) into momentum equation (2.2a), yielding the (dimensionless) master wave-like equation for the displacement vector

∂ 2 ξ 1 ∂t 2 + C ∂ξ 1 ∂t + K(ξ 1 ) = 0, ξ 1 • 1 n = 0 on ∂V, (2.10)
with the linear operators

C(ξ 1 ) = 2 1 Ω × ξ 1 + 2M S (U 0 • ∇) ξ 1 , (2.11a) K(ξ 1 ) = M 2 S (U 0 • ∇) 2 ξ 1 + 2M Ω M S 1 Ω × (U 0 • ∇ξ 1 ) -∇ [∇ • ξ 1 ] (2.11b) -ξ 1 • ∇ 2M Ω M S 1 Ω × U 0 + M 2 S (U 0 • ∇) U 0 -M 2 A [∇ × ∇ × (ξ 1 × 1 B )] × 1 B . (2.11c)
Operator C is skew-Hermitian (Lynden-Bell & Ostriker 1967), whereas operator K is Hermitian without considering our chosen background magnetic field. Finally, we seek modal solutions

ξ 1 (r, t) = ξ(r) exp(λt), ξ • 1 n = 0 on ∂V, (2.12)
where λ = σ + iω is the (complex-valued) eigenvalue with the decay (or growth) rate σ ∈ R and the angular frequency ω ∈ R. Then, problem (2.11) reduces to the quadratic eigenvalue problem

λ 2 ξ + λ C(ξ) + K(ξ) = 0, ξ • 1 n = 0 on ∂V.
(2.13)

MAV consists in reconciling predictions obtained by solving (2.13) with the observed acoustic spectrum, to retrieve the (unknown) background flow.

Polynomial spectral method

We introduce the mathematical result of the paper, that is the new polynomial spectral method in rigid ellipsoids. We will apply it to solve magneto-acoustic problem (2.13). Yet, a wide class of linear problems in (rigid) ellipsoids is readily amenable to the attack presented below. For instance, the hydromagnetic modes (e.g. Vidal et al. 2019a) or various orbitally driven instabilities (e.g. [START_REF] Kerswell | The instability of precessing flow[END_REF][START_REF] Kerswell | Elliptical instability[END_REF][START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF]Vidal & Cébron 2017) can be considered for incompressible fluids. Thus, the polynomial method for compressible flows is worth presenting in general terms. In the rotating frame, we consider the canonical linear problem for a vector field u bounded in a rigid ellipsoid

M ∂ 2 u ∂t 2 + C ∂u ∂t + K(u) = 0 in V, u • 1 n = 0 on ∂V, (3.1)
where [M, C, K] are three differential operators. Problem (3.1) is equipped with the usual inner product

u, v = V u † • v dV, (3.2)
where † denotes the complex conjugate. We can reasonably approximate solutions of (3.1) by using finite-dimensional functional spaces, made of polynomial elements which are square integrable and infinitely differentiable. A similar approach has been pursed in the pioneering work of [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF] for free-surface flows. We present the spectral decomposition for rigid ellipsoids in §3.1, its polynomial expansion in §3.2 and discuss its application to boundary-value problem (3.1) in §3.3.

Spectral decomposition

Assuming that u is square integrable and (at least) two-times continuously differentiable in (3.1), we expand u by using the Helmholtz-Hodge decomposition. This reads

u = ∇Φ + ∇ × Ψ , u • 1 n = 0 on ∂V, (3.3)
with Φ the scalar potential and Ψ the vector potential. The two components ∇Φ and ∇ × Ψ are orthogonal with respect to inner product (3.2). The potentials are found by solving the following problems

∇ 2 Φ = ∇ • u in V, ∇Φ • 1 n = 0 on ∂V, ∇ × (∇ × Ψ ) = ∇ × u in V, (∇ × Ψ ) • 1 n = 0 on ∂V. (3.4)
Solutions of Neumann problem (3.4) for the scalar potential are usually written as the sum of homogeneous and particular potentials Φ = Φ u + Φ w , such that

∇ 2 Φ u = 0, ∇ 2 Φ w = ∇ • u. (3.5a,b)
Non-trivial solutions Φ u = 0 only exist when the normal component of u does not vanish on some part of the boundary [START_REF] Kellog | Foundations of Potential Theory[END_REF]. Hence, we have Φ u = 0 in rigid ellipsoids (e.g. [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF]. Therefore, any vector u bounded in rigid ellipsoids lies in the union of two orthogonal linear vector spaces V W, such that

V : {e = ∇ × Ψ , ∇ • e = 0, e • 1 n = 0 on ∂V }, (3.6a)
W : {e = ∇Φ w , e • 1 n = 0 on ∂V }.

(3.6b)

Note that, for ellipsoids with a free surface (not considered here), any vector lies instead in the union of three orthogonal vector spaces (because Φ u = 0, see [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF]. The associated decomposition for Φ u has been employed in acoustics by Vidal et al. (2019b).

Polynomial vector spaces

We introduce the finite-dimensional space P [l, l], spanned by all Cartesian monomials x i y j z k of degree i + j + k = l. Its dimension is dim P [l, l] = (l + 1)(l + 2)/2 (e.g. [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF]. We denote P [0, n] the finite-dimensional space of all polynomial scalars of degree i + j + k n. We now define the finite-dimensional space P [0, n] of all polynomial vectors for which each component, belonging to P [0, n], is made of Cartesian monomials x i y j z k (of maximum degree i + j + k n). Its dimension is dim P [0, n] = (n + 1)(n + 2)(n + 3)/2 (e.g. [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF][START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF]. It is known that any two-times continuously differentiable field appearing in equation (3.1) can be uniformly approximated by Cartesian polynomials belonging to P [0, n] (Weierstrass approximation theorem). This polynomial decomposition has proven useful for freesurface flows in astrophysics (e.g. [START_REF] Chandrasekhar | Ellipsoidal Figures of Equilibrium[END_REF], and has been rediscovered in elasticity (e.g. [START_REF] Visscher | On the normal modes of free vibration of inhomogeneous and anisotropic elastic objects[END_REF][START_REF] Saviot | Acoustic vibrations of anisotropic nanoparticles[END_REF]. However, these polynomial elements do not naturally satisfy the rigid boundary condition (contrary to the ones introduced below).

Instead, we restrict the vector space V W to the finite-dimensional space P [0, n], by considering basis elements in P [0, n] that satisfy the non-penetration condition. The quest for a polynomial description of V W hinges on the facts that (i) the ellipsoid is a smooth quadratic surface in Cartesian coordinates (1 n is a polynomial vector) and that (ii) some modes supported by equation (2.10) admit polynomial solutions (e.g. [START_REF] Vantieghem | Inertial modes in a rotating triaxial ellipsoid[END_REF][START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF]. We define the restrictions of [V, W] to the finite-dimensional polynomial vector space

P [0, n], denoted V [0, n] ⊂ P [0, n] and W [0, n] ⊂ P [0, n]. We have dim V [0, n] < dim P [0, n] and dim W [0, n] < dim P [0, n].
We build polynomial basis vectors for V [0, n] and W [0, n], involving linear combinations of Cartesian monomials x i y j z k (of maximum degree i + j + k n) that satisfy the nonpenetration condition, as follows.

3.2.1. Basis for V [0, n]
Divergenceless fields represented by the potential vector Ψ in decomposition (3.6a) can be described by two scalar fields, denoted [A, B] in the following. Various descriptions for these scalars have been be proposed, for instance based on different forms of poloidal/toroidal-like decompositions (e.g. [START_REF] Gledzer | Instability of bounded flows with elliptical streamlines[END_REF][START_REF] Wu | High order instabilities of the Poincaré solution for precessionally driven flow[END_REF][START_REF] Ivers | Enumeration, orthogonality and completeness of the incompressible Coriolis modes in a tri-axial ellipsoid[END_REF]. Instead, we define [A, B] as Clebsch (or Euler) scalars. Then, the spectral decomposition for e ∈ V simplifies into e = ∇ × (A∇B) = ∇A × ∇B.

(3.7)

This decomposition has also been introduced to describe quasi-geostrophic motions, almost invariant along the rotation axis (e.g. in spheres [START_REF] Labbé | On magnetostrophic inertia-less waves in quasigeostrophic models of planetary cores[END_REF]. Decomposition (3.7) is of practical interest, since the non-penetration condition is automatically satisfied if either A or B is constant on the boundary ∂V . [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF] found admissible polynomial bases for the scalars [A, B] as follows.

We consider the linearly independent Cartesian monomials in P [0, n-1], that is x i y j z k with i + i + k n -1. Their number is N 2 = n(n + 1)(n + 2)/6. Among them, there are N 1 = n(n + 1)/2 monomials that are independent of z, denoted g i . The other monomials, denoted h i , contain z as factor. We index the set of these polynomials as

{g i } = 1, x, y, x 2 , xy, y 2 , . . . , x n-1 , y n-1 , i ∈ [1, N 1 ],
(3.8a)

{h i } = z, xz, yz, z 2 , . . . , z n-1 , i ∈ [N 1 + 1, N 2 ], (3.8b) such that {p i } = {g i } {h i } with i ∈ [1, N 2 ]
. Then, we define the three pairs of Clebsch scalars (e.g. [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF])

{A = p i F , B = x} with i ∈ [1, N 2 ] ⇒ e i = ∇[p i F ] × 1 x , (3.9a) {A = p i F , B = y} with i ∈ [1, N 2 ] ⇒ e N2+i = ∇[p i F ] × 1 y , (3.9b) {A = g i F , B = z} with i ∈ [1, N 1 ] ⇒ e 2N2+i = ∇[g i F ] × 1 z , (3.9c)
with the shape function

F = 1 -(x/a) 2 -(y/b) 2 -(z/c) 2 . The enumeration of these elements gives the dimension of V [0, n] (see also Backus & Rieutord 2017) dim V [0, n] = N 1 + 2N 2 = n(n + 1)(2n + 7)/6. (3.10)
Hence, V [0, 0] is empty and the polynomial elements of V [0, n] are at least linear in the Cartesian coordinates. Polynomial elements (3.9) are linearly independent [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF]. They are neither normalised nor (fully) orthogonal, but they can be orthonormalised with the (modified) Gram-Schmidt algorithm.

Basis for W [0, n]

In investigating the properties of the inertial modes, [START_REF] Cartan | Sur les petites oscillations d'une masse de fluide[END_REF] outlined an algorithm to build the elements e = ∇Φ w ∈ W [0, n] involving ellipsoidal harmonics.

n = 2, φ = π/2 n = 3, φ = π/2 Figure 2.
Meridional sections of Φw for the potential basis elements ∇Φw ∈ W [0, 3] in an oblate spheroid (a = b = 1, c = 0.9). Each section is taken in a meridional plane containing the z-axis, at the longitude φ measured from the long x-axis. Horizontal axis shows x 2 + y 2 at the longitude φ.

Colour bar shows the scalar potential Φw. Solid grey lines are positive iso-contours, whereas dashed grey lines represent negative iso-contours.

Yet, they admit simple explicit expressions in Cartesian coordinates. We introduce the (non-normalised) normal gradient operator N , defined by

N = n • ∇ = x a 2 ∂ ∂x + y b 2 ∂ ∂y + z c 2 ∂ ∂z , (3.11)
with n = (x/a 2 , y/b 2 , z/c 2 ) ⊤ . For a given monomial x i y j z k ∈ P [1, n + 1], we have

N (x i y j z k ) = i a 2 + j b 2 + k c 2 x i y j z k . (3.12) Since i + j + k 1, we can obtain the inverse operator N -1 such that N -1 N (x i y j z k ) = x i y j z k , yielding N -1 = 1 (i/a 2 + j/b 2 + k/c 2 ) 2 N (3.13) and N -1 (x i y j z k ) = x i y j z k / [i/a 2 + j/b 2 + k/c 2 ] for each Cartesian monomial.
Any polynomial element e = ∇Φ w ∈ W [0, n] is then expressed as a function of the polynomial scalar potential Φ w defined by

Φ w = N -1 (F Ψ w ), F = 1 -(x/a) 2 -(y/b) 2 -(z/c) 2 , (3.14) with Ψ w ∈ P [1, n -1]. We have dim W [0, 0] = dim W [0, 1] = 0 and, for n 2, the one of P [1, n -1]. Thus, we get dim W [0, n] = n(n + 1)(n + 2)/6 -1. (3.15)
Elements e ∈ W [0, n] involve at least quadratic products in the Cartesian coordinates. These elements are neither normalised nor orthogonal, but they can be orthonormalised with the (modified) Gram-Schmidt algorithm. Two elements in W [0, 3] are illustrated in figure 2.

Method of weighted residuals

Equipped with the previous mathematical results, we develop the new global spectral method in ellipsoids for generic boundary-value problem (3.1). We seek a trial solution for the vector field u, projected onto the trial space V [0, n] W [0, n] (by considering polynomial vectors up to the truncation degree n). The trial solution has the form [START_REF] Chang | Natural resonant frequency of a prolate acoustical resonator[END_REF][START_REF] Chang | Natural resonant frequencies of an oblate acoustical resonator[END_REF][START_REF] Willatzen | Eigenmodes of triaxial ellipsoidal acoustical cavities with mixed boundary conditions[END_REF]. In spheres, the solutions combine Bessel functions and spherical harmonics.

u(r, t) = N j=1 α j (t) e j (r), e j • 1 n = 0 on ∂V, (3.16) with N = dim V [0, n]+dim W [0, n], α = (α 1 , α 2 , . . . ) T the (complex-valued) modal coef- ficients and {e j (r)} the real-valued (trial) elements in V [0, n] W [0, n].
Then, we substitute expansion (3.16) in boundary-value problem (3.1). This leads to a set of equations for the modal coefficients, depending on the position (x, y, z) in space. The spatial dependence is removed by using the method of weighted residuals. To do so, we define the residual R(u) of equation (3.1), which measures to what extent the trial expansion is a solution of the original boundary-value problem. When the polynomial truncation degree n is increased in expansion (3.16), the residual becomes smaller to converge towards R(u) = 0 when the true solution u is approached. The latter condition is approximated by projecting, with respect to inner product (3.2), the residual onto a set of test functions {w i (r)}. Note that the choice of test functions is problem dependent. We obtain a set of weighted residual integrals that are set to zero, that is w i , R(u) = 0. This gives the matrix system for the modal coefficients

M d 2 α dt 2 + C dα dt + K α = 0, (3.17)
with the three matrices [M , C, K]. Their elements are given by

M ij = V w † i •M(e j ) dV, C ij = V w † i •C(e j ) dV, K ij = V w † i •K(e j ) dV. (3.18a-c)
Finally, enforcing the (Dirichlet) non-penetration condition in problem (3.17) is not necessary because the boundary condition is automatically satisfied by the basis elements.

Diffusionless magneto-acoustic modes

We are now in a position to solve the magneto-acoustic problem in rigid ellipsoids. We survey illustrative numerical solutions of the eigenvalue problem when M S = 0. The treatment of the physical situation M S = 0, that is required for MAV, is only outlined in appendix B. Although our method can calculate magneto-acoustic modes of arbitrary ellipsoids, we mostly consider spheres in this section (for the sake of simplicity). The polynomial solutions will be also compared and cross-validated with finite-element computations, by using the commercial software COMSOL (see appendix A). Numerical applications related to MAV in ellipsoids will be undertaken in §5. We present the numerical implementation in §4.1. Then, we separately investigate the solutions belonging to the finite-dimensional spaces V [0, n] in §4.2 and W [0, n] in §4.3. The two spaces are combined to study the rotational compressible modes in §4.4 and the magneto-sonic modes in §4.5.

Numerical implementation

We expand ξ onto the trial set of polynomial vectors {e j (r)} of maximum polynomial degree n, as in expansion (3.16), and we apply the method of weighted residuals described in §3.3. We choose the set of test functions {w i (r)} to be equal to the set of trial functions {e i (r)} satisfying the boundary condition. Hence, our variational method is a Galerkin procedure. Problem (2.13) then reduces to the finite-dimensional quadratic eigenvalue problem (QEP) for α

λ 2 M + λ C + K α = 0, (4.1)
with the matrices [M , C, K] obtained from operators (2.11). Their (real) elements are

M ij = V e i • e j dV, C ij = V e i • C(e j ) dV, K ij = V e i • K(e j ) dV. (4.2a-c)
Since the set {e i (r)} is made of Cartesian monomials, volume integrals (4.2) can be evaluated analytically in ellipsoids (see formula 50 in [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF]).

We have implemented the polynomial algorithm within a bespoke numerical code. The matrices [M , C, K] are first computed symbolically, and then are converted to Fortran subroutines (for fast matrix computations from Python with f2py). To reduce the conditioning of the matrices, which affects the numerical accuracy, we have normalised the basis elements with respect to inner product (3.2). This enhances the numerical convergence of the eigenvalue solutions. Yet, we have not orthogonalised the basis elements because the symbolic Gram-Schmidt algorithm is slow (for the large number of elements we have considered). Moreover, it would introduce additional monomials in the polynomial description of each basis element, thereby strongly reducing the numerical efficiency of the Galerkin projections. We have used double-precision arithmetic in the eigenvalue computations, which was found to be sufficient for the modes we are interested in. We have always found that the eigenvalues are purely imaginary λ = iω (up to the machine precision). This is expected in the diffusionless theory and confirmed by the numerics. Thus, we do not discuss the real part of the eigenvalues in the following.

In practice, we convert QEP (4.1) into a generalised eigenvalue problem (GEP) of size 2N . This process is called linearisation [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]. To do so, we recast matrix quadratic problem (4.1) as

λ I 0 0 M α λα = 0 I -K -C α λα , (4.3)
where I is the identity matrix. An obvious difference with QEP (4.1) is that GEP (4. we have used a built-in direct solver for dense matrices based on the generalised Schur decomposition (since we are interested in the full spectrum).

Incompressible modes V [0, n]

By considering the solenoidal space V [0, n] alone, the compressible term vanishes in (2.11) to yield incompressible modes. Among them, the inertial modes play a fundamental role in rotating flows (e.g. Le [START_REF] Bars | Flows driven by libration, precession, and tides[END_REF][START_REF] Zhang | Theory and Modelling of Rotating Fluids: Convection, Inertial Waves and Precession[END_REF]. The inertial modes of a given spatial complexity are exactly described by the finite-dimensional space V [0, n] (Vantieghem 2014). Moreover, [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF] and [START_REF] Ivers | Enumeration, orthogonality and completeness of the incompressible Coriolis modes in a tri-axial ellipsoid[END_REF] proved mathematically that any solenoidal flow in rigid ellipsoids can be expanded onto a linear combination of inertial modes. This rigorously ensures the completeness of the polynomial space V [0, n] in the limit n → ∞ (for square-integrable vector fields).

The spectrum of the inertial modes is illustrated in figure 4. The eigenvalue is purely imaginary λ = iω, with an angular frequency ω which is bounded by |ω|/M Ω < 2 (e.g. [START_REF] Greenspan | The Theory of Rotating Fluids[END_REF]). In compressible units, the frequency grows linearly with M Ω as shown in figure 4(a). Yet, the probability density function of the inertial modes, within the interval |ω|/M Ω < 2, has not been reported yet. We show in figure 4(b) the (sampled) probability density function, obtained by increasing the polynomial degree for values as large as n 55. This suggests that the probability density function of the inertial modes is nearly uniform in spheres. This behaviour starts emerging when considering large enough polynomial degrees (n 55).

Finally, various modes exist in the presence of magnetic fields. The hydrodynamic modes have been already discussed in spheres (e.g. [START_REF] Malkus | Hydromagnetic planetary waves[END_REF][START_REF] Friedlander | Hydromagnetic waves in the Earth's fluid core[END_REF][START_REF] Zhang | Nonaxisymmetric instabilities of a toroidal magnetic field in a rotating sphere[END_REF][START_REF] Labbé | On magnetostrophic inertia-less waves in quasigeostrophic models of planetary cores[END_REF], and do survive in ellipsoids (Vidal et al. 2019a).

Compressible modes W [0, n]

We seek compressible solutions ξ = ∇Φ w ∈ W, by considering non-rotating and nonmagnetic configurations. Then, QEP (4.1) simply reduces to the Helmholtz equation with a Neumann boundary condition, that is

λ 2 Φ w = ∇ 2 Φ w , ∇Φ w • 1 n = 0 on ∂V, (4.4)
with λ = iω the eigenvalue and ω ∈ R the angular frequency. The decay rate is rigorously zero for the non-penetration condition in the absence of diffusive effects [START_REF] Willatzen | Eigenmodes of triaxial ellipsoidal acoustical cavities with mixed boundary conditions[END_REF]. We could seek polynomial solutions of equation (4.4) belonging to W [0, n], by directly applying the method of weighted residuals (see Vidal et al. 2019b, but with a Dirichlet boundary condition on Φ w ). However, this approach is less general than solving (4.1), because it cannot be consistently extended to account for rotation (see appendix C). Nonetheless, discussing the solutions of equation (4.4) is worthy of interest to validate the polynomial solutions against known analytical solutions. Analytical solutions of equation (4.4) can be obtained in spheroids [START_REF] Chang | Natural resonant frequency of a prolate acoustical resonator[END_REF][START_REF] Chang | Natural resonant frequencies of an oblate acoustical resonator[END_REF]) and ellipsoids [START_REF] Willatzen | Eigenmodes of triaxial ellipsoidal acoustical cavities with mixed boundary conditions[END_REF], by using separation of variables in ellipsoidal coordinates. To benchmark the polynomial description, we consider axisymmetric spheroids (a = b = c) for which the eigenfrequencies have been tabulated (only for a few configurations in [START_REF] Chang | Natural resonant frequency of a prolate acoustical resonator[END_REF][START_REF] Chang | Natural resonant frequencies of an oblate acoustical resonator[END_REF]. They are characterised by a triplet (i, l, m) with i an index characterising the radial-like complexity, l the meridional wavenumber (e.g. the spherical harmonic degree in spheres) and m ∈ Z the azimuthal wavenumber (with an azimuthal Fourier decomposition in exp(imφ)). Because of the symmetry of revolution, the modes characterised by the same doublet (i, l) but with different m are degenerate, that is, have the same angular frequency. For a given doublet (l, m), we denote ω i,l,m the angular frequency of the i th spheroidal acoustic mode. This frequency is given by the i th zero of the first derivative (with respect to the radial-like variable) of the prolate (c a) or oblate (a c) spheroidal function of the first kind for a Neumann condition on Φ w . The solutions i = 1 are the fundamental modes, the ones with i = 2 the first harmonics and so on. Contrary to the inertial modes in figure 4, the acoustic modes have a spectrum with a non-uniform (discrete) distribution.

We show in table 2 the comparison between theoretical predictions and polynomial solutions, for a few fundamental modes in spheroids. All the eigenvalues are imaginary, that is λ = iω with ω ∈ R. The real part is numerically zero (up to the machine precision, not shown). Then, we vary the truncation degree n to outline how the numerical

a/c Theory n = 5 n = 10 n = 15 n = 20 (i, l, |m|) ω i,l,m ω ∆ω (%) ω ∆ω (%) ∆ω (%) ∆ω (%)
1 (1, 1, 0) 2.08158 2.08159 8 × 10 -4 2.08158 0 0 0 1 (1, 2, 0) 3.34209 3.34234 7 × 10 -3 3.34209 1 × 10 -8 5 × 10 -13 0 1 (1, 3, 0) 4.51410 4.60590 2 × 10 0 4.51410 2 × 10 -7 0 0 1 (1, 4, 0) 5.64670 ∅ ∞ 5.64673 4 × 10 -4 3 × 10 -12 7 × 10 -13 1 (1, 5, 0) 6.75646 6.48753 4 × 10 0 6.75652 1 × 10 -3 2 × 10 -8 1 × 10 -12 1 (1, 6, 0) 7.85108 ∅ ∞ 7.86149 1 × 10 -1 7 × 10 -8 1 × 10 -11

2 (1, 0, 0) 7.23768 7.24343 8 × 10 -2 7.23769 1 × 10 -4 6 × 10 -11 6 × 10 -11 2 (1, 1, 0) 3.91366 3.91366 3 × 10 -2 3.91254 3 × 10 -10 2 × 10 -10 2 × 10 -10 2 (1, 1, 1) 2.10698 2.10699 7 × 10 -2 2.10698 7 × 10 -10 7 × 10 -10 7 × 10 -10 2 (1, 2, 0) 3.98939 3.98939 4 × 10 -2 3.98799 2 × 10 -7 2 × 10 -10 2 × 10 -10 2 (1, 2, 1) 5.02171 5.02171 8 × 10 -2 5.01764 2 × 10 -6 1 × 10 -10 1 × 10 -10 2 (1, 2, 2) 3.39289 3.39289 7 × 10 -2 3.39264 1 × 10 -8 5 × 10 -10 5 × 10 -10 (up to the truncation degree n = 20). The relative error ∆ω is given by expression (4.5). The symbol ∅ means that the modes do not appear in the spectrum for this polynomial degree.

convergence is affected. We cannot expect the polynomial decomposition to be as efficient as a separable spectral decomposition (for instance using spherical harmonics), which diagonalises exactly Helmholtz equation (4.4). We compute the relative error (between numerics and theory) on the angular frequency

∆ω = |ω -ω i,l,m |/ω i,l,m , (4.5)
as a function of the truncation degree n. Expression (4.5) gives a good proxy of the convergence for eigenvalue computations (e.g. [START_REF] Valdettaro | Convergence and round-off errors in a two-dimensional eigenvalue problem using spectral methods and Arnoldi-Chebyshev algorithm[END_REF]. We find that only the first modes are approximatively described by the polynomial expansion at small degrees (e.g. at n = 5 in spheres), with a relative error smaller than 10%. Modes with higher spatial complexity are not captured at this polynomial degree, for instance the spherical mode (1, 6, 0). Increasing n is required for the modes to appear in the spectrum (here at n = 10). When n is further increased, we observe that the relative errors become extremely small for this subset of modes. The azimuthal degeneracy of the modes is also well recovered, up to the machine precision (when n is large enough, not shown).

Even if the polynomial description is not exact, polynomial solutions of (4.4) are accurate when n is large enough. Its accuracy is rooted in the Taylor expansion of the (known) analytical solutions. In spheres (a = c), spheroidal modes reduce to the textbook solutions in spherical coordinates (r, θ, φ)

Φ w (r, θ, φ) ∝ j l (rω i,l,m ) Y m l (θ, φ), (4.6)
with Y m l (θ, φ) the spherical harmonic (of degree l and azimuthal order m) and the associated spherical Bessel function j l (rω i,l,m ). The latter admits a convergent series in powers of r 2 = x 2 + y 2 + z 2 (e.g. §10 in [START_REF] Abramovitz | Handbook of Special Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. Solution (4.6) has then the infinite polynomial expansion Φ w ∝ r l Y m l T (r 2 ), where r l Y m l is a homogeneous Cartesian polynomial of degree l (i.e. a solid spherical harmonic) and T (r 2 ) a power ). We have used for ω i,l,m either high-precision analytical predictions [START_REF] Chang | Natural resonant frequencies of an oblate acoustical resonator[END_REF] or high-resolution computations at n = 25.

series. Thus, the (almost) arbitrary accuracy of the polynomial description is equivalent to vary the truncation in the power series T (r 2 ). Similar expansions could be obtained in ellipsoids (albeit often not explicitly), due to the Cartesian form of the ellipsoidal harmonics (e.g. [START_REF] Dassios | Ellipsoidal Harmonics: Theory and Applications[END_REF]).

As outlined in table 2, the convergence towards the expected analytical eigenvalues is very fast for the modes with the largest spatial scales. We also quantify the convergence of the eigenvalues with the spatial resolution in figure 5. We show the evolution of relative error (4.5), as a function of the truncation degree n. The polynomial description strikingly exhibits an exponential convergence, which is typical of accurate spectral methods [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF]. The errors decrease very fast, with the numerical scaling ∆ω ∝ exp(-2.4n), until the computations are limited by the round-off errors of the eigenvalue solver (here around 10 -12 -10 -13 ). The errors then depart from an exponential decay, to become almost insensitive to the polynomial degree (as previously reported in table 2 for n = 20). This plateau appears for large enough values of n (depending on the spatial complexity of the modes). No better approximation to the eigenvalue can be obtained by increasing further the truncation degree n with double-precision arithmetic. Yet, the numerical accuracy could be increased with computations using quadruple precision (as considered in [START_REF] Rieutord | Axisymmetric inertial modes in a spherical shell at low Ekman numbers[END_REF]. Finally, the spectral polynomial method outperforms the (standard) finite-element computations performed with COMSOL in terms of numerical accuracy, as shown in appendix A. Indeed, the errors decrease much slower as a function of the spatial resolution, following only a power-law scaling that is characteristic of finiteelement methods using piecewise continuous polynomials [START_REF] Boyd | Chebyshev and Fourier Spectral Methods[END_REF].

Therefore, we have shown that the (spectral) polynomial method is very efficient, due to its spectral convergence. In the following, we have considered truncation degrees n 20, which ensure an excellent convergence for all the modes we are interested in.

Rotational compressible modes

We reintroduce global rotation in the system, by assuming M Ω = 0 and 1 Ω = 1 z . When the system is rotating, the acoustic modes cannot be sought separately from the inertial modes in the polynomial space W [0, n]. Indeed, the incompressible and compressible spaces are now coupled through the Coriolis operator, as illustrated by the non-zero coupled entries of the matrix C in figure 3. Thus, we must seek a priori solutions in

V [0, n] W [0, n].
Yet, in the relevant regime M Ω ≪ 1 (see table 1), we may still expect the acoustic modes to be mainly compressible (i.e. being mostly described by W [0, n]).

We first survey in figure 6 the whole spectrum in rotating spheres, computed by considering solely W [0, 20] in (a) and V [0, n] W [0, 20] in (b). We have also superimposed targeted diffusive computations performed with COMSOL (see appendix A). They are in excellent agreement with the diffusionless polynomial solutions, thereby cross-validating the general results. Several points of figure 6 are worth commenting. We have not shown the range M Ω 10, where we get three families of modes. Two of them are in good agreement with a local analysis (not shown). The missing one results from the rotational splitting of the acoustic modes, which exists only in bounded geometries. Yet, the results for such large rotation rates are unlikely physical. Indeed, the centrifugal force has been entirely neglected in the analysis, although it may affect some acoustic modes for rapid rotation [START_REF] Lignières | Acoustic oscillations of rapidly rotating polytropic stars-I. Effects of the centrifugal distortion[END_REF][START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF]. Moreover, the range M Ω 1 is unreachable in laboratory experiments (see table 1).

When M Ω 1, the two spectra obviously differ by the presence of the inertial modes, associated with the space V [0, 20] in figure 6(b). Their frequency lies in the interval |ω| < 2 M Ω (Valette 1989a,b), as for incompressible fluids (e.g. [START_REF] Backus | Completeness of inertial modes of an incompressible inviscid fluid in a corotating ellipsoid[END_REF]. They are only weakly affected by compressibility when M Ω 1 (not shown). By contrast, the acoustic modes, located above the line |ω| 2 M Ω , are significantly affected by rotation. The latter lifts the azimuthal degeneracy of the spherical acoustic modes. More precisely, a non-rotating spherical solution ω i,l,m splits into 2l + 1 distinct frequencies [START_REF] Gough | The effect of rotation and a buried magnetic field on stellar oscillations[END_REF]. For instance, the fundamental mode with l = 1 splits into three branches (as observed), the second mode with l = 2 into five branches and so on. This splitting is not predicted by a plane-wave analysis (not shown), such that Coriolis effects only occur in bounded geometries.

Magneto-acoustic modes

We finally consider non-rotating magnetic configurations (M A = 0). We investigate the regime M A ≪ 1, which is relevant for experiments (see table 1). Before solving numerically the bounded problem in rigid ellipsoids, it is worth pointing out the properties of the various waves with an unbounded (local) analysis. To do so, we assume that 1 B = 1 z and M Ω = M S = 0. Then, we seek plane-wave solutions ξ ∝ exp(ik • r), with k the (local) wave vector. The typical plane-wave spectrum is shown in figure 7(a). We observe a change of regime, depending on the strength of the background magnetic field (measured by M A ). The transition occurs around M A |k| ≃ 1. In the strong field regime (M A |k| 1), we get three branches of magneto-acoustic waves. The fast branch is made of magnetic (Alfvén) waves, with a dimensionless angular frequency |ω| ∝ M A . The slow waves are almost pure acoustic waves, constrained to move along magnetic field lines (e.g. [START_REF] Campos | On waves in gases. Part II: Interaction of sound with magnetic and internal modes[END_REF]). The intermediate branch represents hybrid magneto-acoustic waves. In the weak field regime (M A |k| 1), these three branches coalesce to yield two branches. The acoustic waves become the fast waves and the Alfvén waves the slow ones.

We illustrate in figure 7(b) the modal spectrum in non-rotating magnetised spheres, as a function of M A . We have also compared our diffusionless polynomial modes with threedimensional (3-D) diffusive computations obtained with COMSOL. Computing magnetic solutions with COMSOL is numerically challenging in three dimensions. Thus, because of numerical constraints, we have fixed the diffusive numbers Re S = 10 3 and Rm S = 10. We have used 416 456 degrees of freedom for these computations. This yields a resolution of approximately 50 points in each spatial direction. We refer the reader to appendix A for additional details on the numerics. We find a broad agreement between the polynomial solutions and the 3-D diffusive solutions (for both the Alfvén modes and the acoustic modes, see below). The discrepancies, usually smaller than 10 -1 , are due to the large diffusive effects (absent in our diffusionless polynomial method) and to the lack of spectral decomposition in the azimuthal direction (in the 3-D finite-element computations). We have checked that we slowly approach the diffusionless solutions by increasing Re S and Rm S (not shown). Then, several points are worth commenting in figure 7(b). We obtain two families of modes. On the one hand, the slow family (when M A 1) disappears by considering only the compressible space W [0, n] (not shown), and their angular frequency scales as ω ∝ M A . This clearly shows that these modes are the magnetic Alfvén modes, mainly associated with the incompressible space V [0, n]. Their eigenfrequencies scale as in the local theory (e.g. [START_REF] Campos | On waves in gases. Part II: Interaction of sound with magnetic and internal modes[END_REF]). On the other hand, the fast modes do survive when we suppress the divergenceless space (not shown). They are acoustic modes, only slightly perturbed by the magnetic field (when M A 1). Moreover, the magnetic field lifts the azimuthal degeneracy of the acoustic part of the spectrum. This was also found by [START_REF] Gough | The effect of rotation and a buried magnetic field on stellar oscillations[END_REF] with perturbation theory, but they considered different background magnetic fields.

Therefore, we have recovered the two families of magneto-acoustic modes in this regime, as predicted by the local theory. We do not study further the modes in rotating magnetised spheres. Indeed, the properties of the incompressible modes should remain largely valid for their compressible counterparts. The Alfvén modes have also been observed and discussed in the experimental context (e.g. in rotating fluids [START_REF] Schmitt | Rotating spherical Couette flow in a dipolar magnetic field: experimental study of magneto-inertial waves[END_REF][START_REF] Tigrine | Torsional Alfvén waves in a dipolar magnetic field: experiments and simulations[END_REF]. Moreover, we will show that magnetic effects are negligible on the acoustic spectrum, such that the interplay between rotation and magnetic field is not worth investigating for MAV.

Acoustic splittings for MAV

We revisit the properties of the acoustic modes in relation to MAV for experiments. The small values of M A in experiments (see table 1) suggest that magnetic fields would have very tiny influence in experiments. We confirm this assumption in appendix B, such that we discard magnetic effects in the following. We investigate quantitatively how the ellipsoidal flattening and rotation along 1 Ω = 1 z alter the (initially degenerate) spherical acoustic modes. Dimensionless variations as small as 10 -4 -10 -3 can be measured in air experiments [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]) and, thus, should be accurately computed. Similarly, we leave aside the effects of a background flow by setting M S = 0. The general method is only presented in appendix B. Quantitative results with M S = 0 are problem dependent and beyond the scope of the present study. Applications for flow inversions will be provided related to the ZoRo experiment [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]). The full ellipsoid is a canonical situation, in which perturbation theory can be carefully assessed against non-perturbative solutions (due to the existence of the polynomial solutions). To quantify ellipticity effects, we introduce the polar flattening ǫ = (ac)/a.

(5.1)

Oblate containers are characterised by 0 ǫ 1 and prolate ones by -1 ǫ 0. We mainly focus on spheroidal geometries (i.e. been discarded. Note that the domain of validity of perturbation computations has been assessed in astrophysics up to ǫ 0.16 [START_REF] Lignières | Acoustic oscillations of rapidly rotating polytropic stars-I. Effects of the centrifugal distortion[END_REF][START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF], by increasing global rotation (which is responsible for spheroidal shapes due to centrifugal effects with a free-surface condition). Yet, these works cannot be used to disentangle easily the effects due to the ellipticity and rotation in rigid ellipsoids. Indeed, the rotation rate and the flattening are two independent parameters in experimental conditions. Hence, their effects remain to be carefully evaluated for MAV.

We denote in the following ω(ǫ, M Ω ) the angular eigenfrequency of the consistent solutions (without approximation) and ω(0, 0) the ones of the acoustics modes in nonrotating spheres (ǫ = 0, M Ω = 0). Then, ω(0, 0) is shifted by small amounts due to flattening δ ǫ and rotation δ Ω along 1 Ω = 1 z . In the perturbation framework, these effects are often treated separately. If this assumption were consistent for experimental conditions, we would obtain

ω(ǫ, M Ω ) ≃ ω(0, 0) + δ ǫ + δ Ω .
(5.2)

The residual of expression (5.2), denoted δ ǫΩ , is due to any cross-effect between ellipticity and rotation. The frequency shifts would be given by

δ ǫ = ω(ǫ, 0) -ω(0, 0), δ Ω = ω(0, M Ω ) -ω(0, 0). (5.3a,b)
We assess separately ellipticity effects δ ǫ in §5.1 and rotational ones δ Ω in §5.2. They are then combined in §5.3 to assess approximation (5.2), which neglects any possible cross-effect δ ǫΩ between rotation and flattening.

Flattening effects

We investigate quantitatively the ellipticity splitting of the compressible modes belonging to W [0, n], to assess the validity of perturbation theory in ellipsoids. Ellipticity effects are usually modelled with first-order perturbation theory (see in metrology [START_REF] Guianvarch | Acoustic field in a quasispherical resonator: unified perturbation model[END_REF]), but [START_REF] Mehl | Acoustic eigenvalues of a quasispherical resonator: second order shape perturbation theory for arbitrary modes[END_REF] extended the theory up to the second order. However, moderate ellipsoidal deformations are often encountered in fluid experiments, for which the perturbation approximation may be inaccurate. The various modes given in table 2 for spheroidal (axisymmetric) geometries have illustrated that an ellipsoidal boundary lifts the azimuthal degeneracy in |m| of the spherical modes, for a given doublet (i, l). Indeed, the degeneracy of the modes (1, 2, |m|), with |m| ∈ [1, 2, 3], is lifted as the spherical cavity is squashed. Yet, the degeneracy of positive and negative m is not removed in spheroids (not shown).

We show in figure 8(a) how ellipticity affects several low-frequency acoustic modes, as a function of the polar flattening ǫ in spheroids. We clearly observe that the degeneracy of the spherical modes (ǫ = 0) is lifted by ellipsoidal flattening. Yet, the different branches, each associated with a given azimuthal number m in the spheroidal geometry, are not modified in the same way. A naive first-order theory would predict that the frequency shift is linear in ǫ. This is obviously not the case for some modes at moderate deformations |ǫ| 0.2. Some curves are closer to parabolas than to straight lines, such that they could be fairly described by second-order perturbation theory [START_REF] Mehl | Acoustic eigenvalues of a quasispherical resonator: second order shape perturbation theory for arbitrary modes[END_REF]). In the strong deformation limit, the acoustic branches tend to diverge and crossings can occur between branches associated with different azimuthal modes in the initial spherical geometry (e.g. near ω = 4.5 at ǫ ≃ 0.4 in figure 8a). Finally, we illustrate some triaxial configurations in figure 8(b) (for the first low-frequency modes), by assuming b = a(1 + ǫ) and c = a(1ǫ) with a = 1 (in dimensionless units). The curves are symmetric with respect to ǫ = 0, since the acoustic problem is invariant under an exchange of b and c. As naively expected for these small deformations (|ǫ| 0.05), the branches are close to straight lines, in agreement with first-order theory for the ellipticity splitting.

We have shown that ellipticity effects do not affect similarly all the acoustic branches in figure 8. The branches tend to diverge by increasing the flattening. Thus, accidental degeneracies of acoustic modes (with initially distinct frequencies in spheres) can occur for large deformations. Unexpectedly, accidental degeneracy also occurs for initially identical branches, due to high-order ellipticity effects. This striking feature is clearly illustrated in figure 9, for two low-frequency acoustic modes in spheroidal geometries. The m = 1 branches are straight lines, that is they are mainly affected by first-order ellipticity effects. However, the m = 0 modes undergo second-order effects (i.e. in ǫ 2 ), such that the m = 0 branches cross the m = 1 branches at small values of ǫ. This is not an artefact of the polynomial description, since we also get this behaviour with the finiteelement computations in COMSOL. The numerical agreement is excellent, because we can perform diffusionless computations without rotation in COMSOL (see appendix A). Thus, figure 9 completely cross-validates the diffusionless results, which are recovered by using the two methods. These points are not avoided crossings that may exist in spheroidal geometries [START_REF] Lignières | Acoustic oscillations of rapidly rotating polytropic stars-I. Effects of the centrifugal distortion[END_REF], since each branch has its own azimuthal symmetry. Moreover, it turns out that the crossing point in figure 9 (left panel) occurs near ǫ ≃ 0.05, that is the spheroidal deformation of the ZoRo experiment [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF].

We illustrate in figure 10 the modes which are close to the accidental degeneracies observed in figure 9 (left panel). We have represented the dimensionless acoustic pressure p 1 ∝ -∇ • ξ. The two modes have large spatial scales, but different meridional structures and azimuthal numbers (m = 0 and m = 1). Note that these modes do not undergo avoided crossing. Indeed, they keep their own azimuthal (and longitudinal) structure from each side of the crossing point (not shown). Finally, we point out that large-scale modes are likely sensitive to the large-scale components of the background velocity field, which are (a priori) easier to determine with MAV. Thus, a poor identification of these two particular low-frequency modes in the experimental data would be certainly misleading for the acoustic inversion.

This unexpected phenomenon clearly shows that an accurate description of the acoustic modes is a prerequisite to any application of MAV (at least in moderately deformed ellipsoids). In particular, if one would like still to use domain perturbation theory in ellipsoids, we advocate to employ (at least) a second-order theory [START_REF] Mehl | Acoustic eigenvalues of a quasispherical resonator: second order shape perturbation theory for arbitrary modes[END_REF]. , νB/ν = 0.6). The experimental range of parameters is MΩ 2 × 10 -1 (see table 1).

10 -2 10 -1 M Ω 1.6 1.8 2.0 2.2 2.4 ω m = -1 m = 1 m = 0 V W W COMSOL 10 -2 10 -1 M Ω 3.0 3.2 3.4 3.6 ω m = -2 m = 2 m = -1 m = 0 m = 1 V W W COMSOL

Rotational splitting

We now focus on the frequency shift due to rotation δ Ω in the absence of deformation. Coming back to figure 6, we observe that the rotational shift seems to be accurately described solely by the compressible space W [0, n] for small enough rotations. Indeed, the curves seem poorly distinguishable when M Ω ≪ 10 -1 , whereas we do notice some differences when M Ω 10 -1 (for the intermediate branch of the first mode). This behaviour is confirmed in figure 11, in which we have superimposed the computations with W [0, 20] and V [0, 20] W [0, 20]. We illustrate only the two first modes, but we have checked that this behaviour remains valid for modes with higher frequencies (if their polynomial descriptions are well converged). The polynomial results are also in very good agreement with the 3-D computations in COMSOL (even if the latter necessarily include diffusion). This confirms the robustness of the results, even if the diffusive angular frequencies are slightly smaller (in absolute value) than the diffusionless frequencies. This is consistent with the diffusive theory [START_REF] Moldover | Gas-filled spherical resonators: Theory and experiment[END_REF]). Thus, acoustic modes could be fairly described in the range M Ω ≪ 10 -1 by considering solely the space W [0, n].

The frequency shift δ Ω is often estimated by first-order perturbation theory (e.g. [START_REF] Gough | The effect of rotation and a buried magnetic field on stellar oscillations[END_REF][START_REF] Dahlen | Theoretical Global Seismology[END_REF]. Without centrifugal effects and when M Ω ≪ 1, we may expect the acoustic frequencies to scale linearly with rotation (i.e. δ Ω ∝ M Ω ). Observing linear or quadratic effects (δ Ω ∝ M 2 Ω ) is not easy in figure 11, due to the log-linear representation. To disentangle the two scalings, we show in figure 12 the evolution of |δ Ω |/M Ω , as a function of M Ω in the range 0 M Ω 10 -1 . For a dominant linear scaling, |δ Ω |/M Ω should be roughly constant. We have illustrated the observed behaviour for some low-frequency modes, which are representative of the results. We do obtain a linear variation |δ Ω | ∝ M Ω when M Ω 10 -2 , in agreement with first-order theory. When M Ω 10 -2 , the typical error on the ratio |δ Ω |/M Ω is smaller than 10 -3 , yielding an upper-bound error on the frequency shift of 10 -5 . Then, second-order Coriolis effects are already present when M Ω 10 -2 . At M Ω = 10 -1 , the typical error is 3 × 10 -2 , yielding as typical estimate for the error on the frequency shift 3 × 10 -3 . Finally, we recover that the rotational shift is fairly described by considering space W [0, n] when M Ω ≪ 10 -1 . Moreover, high-order modes are less affected by the Coriolis force (not shown), as found experimentally [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]) and in agreement with astrophysical predictions [START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF]. We thus expect small Coriolis effects in the high-frequency regime. Consequently, the linear approximation seems accurate enough for typical experimental purposes when M Ω 10 -2 . However, the linear approximation may not be sufficient for some large-scale modes and for large values of M Ω , which can be considered in gas experiments (e.g. [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF].

Interplay between rotation and flattening

We have shown that, for typical experimental conditions (see table 1), ellipsoidal flattening is quickly responsible for second-order variations in ǫ, even for values as small as ǫ 0.025 (see figure 9). In addition, second-order Coriolis effects in M Ω might be also observed, notably in gas experiments (for which M Ω 2 × 10 -1 ). This leads us to naturally question the relevance of formal perturbation expansion (5.2), which treats independently ellipticity and rotational effects. Indeed, this approach rules out any possible cross-effect between rotation and flattening δ ǫΩ , which would scale as ǫM Ω at leading order. These effects may not be vanishingly small for experimental conditions, in which we can measure slight frequency changes of order 10 -4 -10 -3 [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF].

Accurate estimates of δ ǫ and δ Ω are required to assess generic formula (5.2). We have shown that the polynomial method describes precisely δ ǫ and δ Ω at any order. This is strongly different in perturbation theory. For instance, δ ǫ is usually limited to the second order in the deformation [START_REF] Mehl | Acoustic eigenvalues of a quasispherical resonator: second order shape perturbation theory for arbitrary modes[END_REF][START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]. However, second-order ellipticity effects are often not accurate enough (e.g. when ǫ = 0.3). We clearly illustrate this point in figure 13. We have first determined the rotating spherical solutions ω(0, M Ω ) in V [0, 20] W [0, 20]. Second, we have computed the frequency shift δ ǫ , in the absence of rotation, with formula (5.3). Each acoustic branch, characterised here by a given azimuthal number m, has its own frequency shift δ ǫ . We have only considered the secondorder ellipticity effects in δ ǫ , which would be obtained with perturbation theory. To isolate these second-order effects, we have used a quadratic polynomial fit in powers , 20]. Black squares: perturbation solutions (5.2), obtained by adding the frequency shift δǫ (due to ellipticity) and the rotating spherical solutions ω(0, MΩ). Frequency shift δǫ for each branch has been computed with formula (5.3), by using the non-rotating polynomial solutions shown in (a). Rotating spherical modes ω(0, MΩ)

computed in V [0, 20] W [0, 20].
of ǫ (dashed line in figure 13a), to describe the polynomial solutions (when ǫ ≪ 0.1). On the one hand, the fits are close to the true solutions when ǫ 0.1 (as naively expected). On the other hand, the fitted curves are far from the true branches for the upper acoustic branches when ǫ 0.1. It becomes obvious that second-order theory is largely inaccurate for the flattening ǫ = 0.3, even for the low-frequency modes. Computing δ ǫ from the quadratic law yields large errors in formula (5.2), as observed in figure 13(a). For the upper branches, the accurate estimates (containing all ellipticity effects) are δ ǫ (ǫ, 0) = [0.4991, 0.7121], whereas the fitted values are δ ǫ (ǫ, 0) ≃ [0.5816, 0.6518]. These discrepancies would be directly responsible for erroneous predictions for the acoustic branches, as illustrated by the large offsets on the upper curves in figure 13(b).

Then, the goal is to evaluate the independence of ellipticity and Coriolis effects in expression (5.3). We assume that δ ǫ contains all the effects due to ellipticity (without approximation) in the absence of rotation, whereas δ Ω accounts for all the Coriolis effects (i.e. without the centrifugal effects, which have been entirely neglected) in the absence of ellipsoidal deformation. We assess formula (5.2) in figure 14, by considering two polar flattenings ǫ = 0.05 (e.g. [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]) and ǫ = 0.3 (e.g. [START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF][START_REF] Lemasquerier | Libration-driven flows in ellipsoidal shells[END_REF]. We have computed the frequency shift δ ǫ in the absence of rotation with formula (5.3). Each acoustic branch, characterised here by a given azimuthal number m, has its own frequency shift δ ǫ . To avoid neglecting high-order ellipticity effects, the frequency shift δ ǫ is computed for each branch from accurate computations of ω(ǫ, 0) and ω(0, 0) in W [0, 20]. We present illustrative results for some low-frequency modes.

We observe in figure 14(a) that perturbation solutions (5.2) seem in good agreement with the consistent solutions ω(ǫ, M Ω ) when ǫ = 0.05. However, this is largely Blue triangles: COMSOL diffusive computations (separate m, ReS = 4 × 10 6 , νB/ν = 0.6). The experimental range of parameters is MΩ 2 × 10 -1 (see table 1).

10 -2 10 -1 M Ω 3.1 3.2 3.3 3.4 3.5 3.6 ω m = -2 m = 2 m = -1 m = 0 m = 1 V W Perturbation COMSOL 10 -2 10 -1 M Ω 3.2 3.4 3.6 3.8 4.0 4.2 ω m = 0 m = 2 m = -2 m = 1 m = -1 V W Perturbation COMSOL (
a consequence of the chosen graphical representation. Quantitatively, we do obtain small discrepancies between the two approaches. They are due to cross-effects δ ǫMΩ in this figure, because the frequency shift δ Ω contains any order of the Coriolis effects (independently of the ellipticity). Note that the error amplitude is also mode-sensitive. Their typical magnitude is δ ǫMΩ ≃ c 1 (M Ω ) ǫM Ω , where the numerical prefactor is c 1 (M Ω ) 0.5 when M Ω 10 -1 . For instance, we get for the mode with m = ±1 at M Ω = 10 -1 in figure 14(a) a frequency shift around 1 Hz (dimensional units), that is a dimensionless error of 4 × 10 -3 on the angular frequency (not shown). Hence, we theoretically predict that cross-effects δ ǫΩ are indeed observable in experiments. The situation is even much clearer for the ellipsoidal configuration ǫ = 0.3 illustrated in figure 14(b). For instance, the upper curves around ω = 4.1 differ in angular frequency by already 10 -2 at M Ω = 0.1, which is one order of magnitude larger than for the lower branches. The frequency shift is then typically 10 -3 when M Ω ≃ 10 -2 for the upper branches and 10 -4 for the lower ones in figure 14(b). Thus, we support that cross-effects between ellipticity and rotation are observable given the experimental precision (at least for some modes), even if they are discarded with perturbation theory [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF].

Experimental implications

Previous results have direct implications for experiments. As shown in appendix B, magnetic splittings are far too small and can be safely ignored in the analysis. Yet, ellipticity and Coriolis effects can be significant and must be considered with care. For the ZoRo set-up (ǫ = 0.05), [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF] conducted a preliminary experimental survey with rotation, by using air at ambient temperature (C * 0 ≃ 343 m/s). They confronted their experimental observations to a mixed perturbation theory similar to equation (5.2), taking into account second-order ellipticity effects and first-order Coriolis effects (without ellipticity). They considered (dimensional) rotation rate frequencies Ω * s /(2π) 30 Hz, that is M Ω 10 -1 in dimensionless units. The quantitative results outlined in figure 14(a) show that cross-effects are observable at M Ω = 10 -1 . In particular, as also done in astrophysics (e.g. [START_REF] Aerts | Asteroseismology[END_REF], [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF] measured the absolute acoustic splitting δ m ω. The latter is the difference between the angular frequency of the acoustic branches with azimuthal wavenumbers m and -m, that is δ m ω = |ω(m)ω(-m)|. Considering first-order perturbation theory for the Coriolis effects (without ellipticity), they observed small differences with their experimental data for the |m| = 1 and |m| = 2 branches shown in figure 14(a), and these discrepancies do quantitatively agree with our predictions for the cross-effects δ ǫΩ .

In addition, we also predict, for this experiment, that second-order Coriolis effects should be observable in the experimental data for M Ω 10 -1 . They are illustrated in figure 15(a). Second-order Coriolis effects could be detected by directly measuring the angular frequency ω of the modes. We have superimposed, on the exact solutions in figure 15(a), the linear fits (i.e. in M Ω ) of ω that model both the first-order Coriolis effects and cross-effects (because the latter scale as ǫM Ω ). Second-order effects in M Ω , that is secondorder Coriolis effects at leading order, are responsible for the departures between the exact and the fitted curves. Yet, we point out that second-order Coriolis are mainly filtered out by measuring the absolute splitting δ m ω, as considered in [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]. As shown in figure 15(b), the splitting δ m ω obtained with the fitted curves nearly coincides with the ones computed from the exact solutions. Indeed, by looking at the Taylor expansion in M Ω of ω, the prefactor of the term in M 2 Ω (i.e. the second-order Coriolis effects) has the same sign and roughly the same amplitude between the two split branches with |m| (not shown). Consequently, the second-order Coriolis effects are filtered out by considering the absolute splitting δ m ω. However, they could be isolated from the cross-effects by considering instead the sum of the acoustic branches (not shown). For these reasons, the spectral polynomial method should be preferred for future MAV applications in rigid ellipsoids. This would avoid the cumbersome development of high-order perturbation theory, to describe rotational and ellipticity effects.

In moderately deformed ellipsoids (e.g. [START_REF] Noir | Experimental evidence of inertial waves in a precessing spheroidal cavity[END_REF][START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF][START_REF] Lemasquerier | Libration-driven flows in ellipsoidal shells[END_REF][START_REF] Le Reun | Experimental study of the nonlinear saturation of the elliptical instability: inertial wave turbulence versus geostrophic turbulence[END_REF], the situation is even worse. The illustrative comparison between figures 14 and 15 clearly demonstrates that second-order domain perturbation theory [START_REF] Mehl | Acoustic eigenvalues of a quasispherical resonator: second order shape perturbation theory for arbitrary modes[END_REF] is not often sufficient to model ellipticity effects. The predictions for some acoustic modes would be largely inaccurate and would pollute the velocity reconstruction in MAV. Instead, the polynomial method is perfectly suited to predict the resonant frequencies in the presence of global rotation and any deformation. These examples confirm the relevance of the polynomial method (compared to perturbation method) in acoustic studies of rigid ellipsoids. Indeed, its accuracy outperforms perturbation theory to determine the acoustic modes. This would prevent from introducing (avoidable) modelling errors, which may pollute the interpretation of the experimental data in MAV.

In addition to figures 14(a) and 15, already computed for the ZoRo geometry, we show in figure 16 theoretical acoustic spectra for some acoustic modes that could be detected during the forthcoming experimental surveys with the ZoRo apparatus. In the range M Ω ≪ 1, we predict accidental degeneracies for several acoustic branches (with different azimuthal wavenumbers) due to global rotation in the flattened geometry (blue squares). These crossings are worth tracking in the experimental data. Note that avoided crossings may also occur, for high-frequency acoustic modes sharing the same azimuthal symmetry [START_REF] Lignières | Acoustic oscillations of rapidly rotating polytropic stars-I. Effects of the centrifugal distortion[END_REF][START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF]. The spatial complexity of two of these degenerate modes is illustrated in figure 17. If these crossings were overlooked, this would alter the results from MAV technique (which is sensitive to the large-scale modes). Moreover, the reconstruction of the flow components with MAV is much harder in the presence of accidental degeneracy of the acoustic modes. Indeed, the usual perturbation formula employed to take a background flow into account (e.g. [START_REF] Aerts | Asteroseismology[END_REF] ought to be modified in the presence of these degenerate cases (as explained in appendix B).

Conclusion

Summary

We have addressed the problem of compressible modes in rotating rigid ellipsoids. We have conducted an interdisciplinary theoretical and numerical study, related to planetary-driven flow experiments [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]. Indeed, the compressible modes could be used in modal acoustic velocimetry, to reconstruct passively the three-dimensional components of the (large-scale) velocity field in ellipsoidal experiments. To this end, an accurate description of the resonant acoustic frequencies of the fluid cavity is required. Indeed, perturbation theories are often mathematically involved (e.g. [START_REF] Moldover | Gas-filled spherical resonators: Theory and experiment[END_REF][START_REF] Gough | The effect of rotation and a buried magnetic field on stellar oscillations[END_REF][START_REF] Mehl | Acoustic eigenvalues of a quasispherical resonator: second order shape perturbation theory for arbitrary modes[END_REF] and have also (very) restricted domains of validity (see [START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF], for astrophysical flows).

We have developed a new spectral method in rigid ellipsoids. It is valid for any vector field, satisfying the non-penetration condition. This relies on the orthogonal expansion onto a divergenceless vector space V and an irrotational vector space W. We have constructed admissible square-integrable fields out of global Cartesian polynomials, which are infinitely differentiable and exactly satisfy the non-penetration boundary condition. We have emphasised the completeness of the orthogonal decomposition, ensuring that no admissible field is suppressed. We have favoured a systematic exposure of the procedure, since we believe this approach to be useful beyond experimental problems. Indeed, the method only involves clever mathematics in Cartesian coordinates to attack physical problems that are cumbersome to solve in ellipsoidal (e.g. [START_REF] Cartan | Sur les petites oscillations d'une masse de fluide[END_REF] or non-orthogonal coordinates (e.g. [START_REF] Bonazzola | Numerical approach for high precision 3D relativistic star models[END_REF]. Then, we have applied this new method to compute the magneto-acoustic modes. They are solution of a quadratic eigenvalue problem, formulated for the fluid particle displacement vector. The numerical validation (against finite-element computations) has shown that the polynomial spectral method has an exponential convergence. This is a desirable property for accurate predictions in MAV, which is sensitive to slight acoustic variations.

Finally, we have carefully assessed perturbation theory for MAV in ellipsoids, as employed in the preliminary experimental study of [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF] with the ZoRo apparatus. Indeed, the ellipsoidal configuration offers a canonical situation in which diffusionless polynomial solutions exist to be used as benchmarks. We have investigated how global rotation and the ellipsoidal flattening disturb the acoustic spectrum, since magnetic effects were entirely negligible. They lift (even partially) the azimuthal degeneracy of the acoustic modes, which exist in non-rotating spheres. We have shown that ellipticity effects cannot be predicted by first-order perturbation theory [START_REF] Guianvarch | Acoustic field in a quasispherical resonator: unified perturbation model[END_REF], even in weakly deformed contained. Indeed, second-order effects quickly appear by increasing the polar flattening ǫ, even for values as small as ǫ = 0.05 as encountered in the ZoRo experiment [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]. They are also responsible for accidental degeneracies of some low-frequency acoustic modes for typical experimental deformations. Then, we have investigated the frequency shift due to rotation, which is a (diffusionless) mechanism that exists only in bounded geometries. For most experimental regimes, with small enough fluid rotation rates compared to the speed of sound (i.e. M Ω ≪ 10 -1 ), rotational effects vary mostly linearly with the rotation rate. Yet, we have shown that second-order Coriolis effects could be observed (for some modes) in more extreme experimental conditions, in agreement with preliminary experimental results [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]. Similarly, we have quantified the cross-effects between rotation and ellipticity for experimental conditions. They appear to be negligible for (i) small rotations (typically M Ω ≃ 10 -2 ) and (ii) small enough deformations. However, we have shown that they are responsible for the observed small variations between perturbation theory and the experimental observations reported in [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF], for rotation rates 10 -2 M Ω 10 -1 .

Perspectives

Despite being rather idealised, full ellipsoids are simple analogues of planetary fluid cores. They offer a rich flow dynamics, which can be simulated in fluid experiments and could be probed with MAV. In addition to spheroidal geometries (e.g. [START_REF] Noir | Experimental evidence of inertial waves in a precessing spheroidal cavity[END_REF][START_REF] Noir | Experimental study of libration-driven zonal flows in non-axisymmetric containers[END_REF][START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF], future experimental applications should cover further triaxial geometries to mimic tidally deformed fluid cores (e.g. [START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF][START_REF] Lemasquerier | Libration-driven flows in ellipsoidal shells[END_REF][START_REF] Le Reun | Experimental study of the nonlinear saturation of the elliptical instability: inertial wave turbulence versus geostrophic turbulence[END_REF]. The triaxial geometry would naturally lift the azimuthal degeneracy of the spheroidal acoustic modes. It would also introduce additional accidental degeneracies, for instance due to cross-effects between rotation and ellipticity. This would complicate further the interpretation of observed acoustic frequencies. The comparison with the polynomial solutions has showed that the perturbation framework (e.g. [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF] is not always accurate enough to match the experimental conditions in deformed geometries. Therefore, we hope the present study will be an impetus in favour of the polynomial method for future applications of MAV in full (triaxial) ellipsoids.

Additional ingredients would be worth including in the physical model, to improve the accuracy of MAV. Since the full acoustic problem is rather complex, we have only considered without approximation the Coriolis and ellipticity effects in ellipsoids. However, geometrical imperfections are hardly avoidable in experiments [START_REF] Moldover | Gas-filled spherical resonators: Theory and experiment[END_REF][START_REF] Mehl | Acoustic eigenvalues of a quasispherical resonator: second order shape perturbation theory for arbitrary modes[END_REF][START_REF] Guianvarch | Acoustic field in a quasispherical resonator: unified perturbation model[END_REF]. Small departures from ellipsoidal boundaries could be modelled in our framework, by using domain perturbation theory [START_REF] Lebovitz | Perturbation expansions on perturbed domains[END_REF]. We have also neglected all diffusive effects, which are usually small in experiments (as well as in celestial objects) except in the boundary layers [START_REF] Berggren | Acoustic boundary layers as boundary conditions[END_REF]. However, diffusive effects may be comparable in amplitude with the cross-effects between ellipticity and rotation. Thus, it would be worth investigating how diffusion alters the acoustic spectrum. Viscous and non-adiabatic effects at the boundary can be approximated with perturbation theory in non-rotating spheres [START_REF] Moldover | Gas-filled spherical resonators: Theory and experiment[END_REF]). Hence, we could estimate the diffusive effects in weakly deformed ellipsoids by using the diffusive spherical solutions. Alternatively, we could directly solve, in non-rotating ellipsoids, the diffusive scalar equation for the temperature (see equation 6.4.22 in [START_REF] Morse | Theoretical Acoustics[END_REF]. We would project it onto another polynomial space, satisfying the appropriate boundary condition (see a closely related problem in Vidal et al. 2019b). Then, as a long-term endeavour, the full acoustic spectrum could be synthetically reproduced and used for the flow inversion [START_REF] Moldover | Gas-filled spherical resonators: Theory and experiment[END_REF][START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF].

Accounting for spatially varying background states appears also desirable in several contexts. In acoustics, non-isentropic background states (as driven by thermal diffusion in experiments) modify the acoustic spectrum. A first step has been achieved by [START_REF] Koulakis | Acoustic resonances in gas-filled spherical bulb with parabolic temperature profile[END_REF]. They found the analytical description of the spherical acoustic modes with a parabolic temperature profile. Yet, their method cannot be extended (a priori) to other profiles (without using perturbation theory), and also not in the presence of an ellipsoidal boundary. The polynomial description provides a natural way to investigate this problem, even in spherical geometry. Indeed, spectral decomposition (3.3) and its polynomial description are not limited to spatially uniform background states.

Similarly, we have entirely neglected the centrifugal effects. Because of this omission, we have not completely treated all rotational effects. More generally, centrifugal effects have been poorly studied in experimentally driven studies, except for instance in [START_REF] Horn | Regimes of Coriolis-centrifugal convection[END_REF], 2019) and [START_REF] Menaut | Experimental study of convection in the compressible regime[END_REF]. Firstly, they modify the shape of celestial fluid bodies (e.g. [START_REF] Chandrasekhar | Ellipsoidal Figures of Equilibrium[END_REF][START_REF] Rieutord | An algorithm for computing the 2D structure of fast rotating stars[END_REF]. Secondly, the centrifugal force could modify the frequencies of the acoustic modes. The gravest modes would be largely unaffected by the centrifugal force [START_REF] Ecotiere | Inertial coupling of resonant normal modes in rotating cavities: Acoustic gyrometers for high rotation rates[END_REF]). Yet, contrary to the Coriolis force, the effects of the centrifugal force may increase with the radial-like complexity [START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF]. The spatial structure of the high-frequency acoustic modes may be also altered, with an equatorial concentration of the wave energy (for astrophysical flows, see [START_REF] Lignières | Acoustic oscillations of rapidly rotating polytropic stars-I. Effects of the centrifugal distortion[END_REF]. Moreover, avoided crossing might be also favoured for rapid rotation (e.g. [START_REF] Lignières | Acoustic oscillations of rapidly rotating polytropic stars-I. Effects of the centrifugal distortion[END_REF][START_REF] Reese | Acoustic oscillations of rapidly rotating polytropic stars-II. Effects of the Coriolis and centrifugal accelerations[END_REF]. Thus, we should strive including exactly centrifugal effects within the polynomial approach, by considering a spatially varying background state and an effective gravitational force. For all the aforementioned reasons, we advocate to use the polynomial method for future applications of MAV in full ellipsoids. Finally, we also believe that the polynomial method would hold out promise for accurate planetary-driven reduced models, to go beyond the Boussinesq (e.g. [START_REF] Vidal | Magnetic fields driven by tidal mixing in radiative stars[END_REF] and anelastic (e.g. [START_REF] Clausen | Elliptical instability of compressible flow in ellipsoids[END_REF] descriptions for the flow dynamics. Computations at ReS = 4 × 10 6 and ratio νB/ν = 0.6 in a sphere (isothermal boundary at 20 • C). Solid lines show ∆ω and dashed lines the relative error on the decay rate ∆σ. The reference solutions have been computed by using quintic elements for the both the velocity and the pressure (P5-P5) and #DOF= 880644. Black triangles: elements P2-P1 (quadratic for the velocity, linear for the pressure). Red circles: elements P3-P2 (cubic for the velocity, quadratic for the pressure). Green squares: elements P4-P3 (quartic for the velocity, cubic for the pressure).

and σ. We have checked that the diffusive solutions approach the diffusionless ones when the diffusion is reduced, by gradually increasing the numerical resolution (not shown).

In non-magnetic rotating cases, we have solved equations (2.5) with Lagrange P4-P3 elements (i.e. quartic for the velocity and cubic for the pressure).

A.2. Magnetic computations (M A = 0) Lagrange elements are not suited to solve induction equation (2.5c). We employ Nédélec's elements for the magnetic field, which have been already validated for hydromagnetic computations (e.g. [START_REF] Cébron | Magnetohydrodynamic simulations of the elliptical instability in triaxial ellipsoids[END_REF]. We have supplemented the built-in adiabatic (and viscous) formulation of COMSOL with the induction equation. We have performed only 3-D computations in this case (even in axisymmetric containers). We assume that the rigid boundary is a perfect conductor, in which no magnetic field is trapped. Hence, we consider the magnetic boundary conditions b 1 • 1 n = 0 and (∇ × b 1 ) × 1 n = 0 on ∂V.

(A 2)

With perfectly conducting conditions (A 2), the magnetic boundary layer has weaker effects than for finite values of the electrical conductivity (e.g. [START_REF] Roberts | On the diffusive instability of some simple steady magnetohydrodynamic flows[END_REF]. Consequently, we can expect the diffusive magnetic computations to be closer to the diffusionless theory. We have also replaced induction equation (2.5c) by the equation for the magnetic vector potential b 1 = ∇ × Λ 1 . This formulation enforces the divergenceless condition of the magnetic field in finite-element simulations (e.g [START_REF] Cébron | Magnetohydrodynamic simulations of the elliptical instability in triaxial ellipsoids[END_REF]).

The leading-order terms (i.e. M S = 0) have been considered in the main text. They exactly account for the key physical ingredients, namely compressibility and global rotation. The background flow is introduced at the next asymptotic order (M S ≪ 1). For non-degenerate eigenvalues λ (0) , the first-order correction λ (1) obeys a variational principle given by the solvability condition of the system at the order M S . Since the zeroth-order (infinite-dimensional) operators in (B 1) are Hermitian (e.g. [START_REF] Lynden-Bell | On the stability of differentially rotating bodies[END_REF], we get the variational principle in ellipsoids At the leading asymptotic order (M S = 0), the unperturbed quantities [λ (0) , α (0) ] are solution of the QEP given in the main text. At the next asymptotic order, the solvability condition gives (e.g. formula 4.10 in Seyranian 1993)

λ (1) α † (0) ⊤ • 2λ (0) M + C (0) α (0) = -α † (0) ⊤ • λ (0) C (1) + K (1) α (0) . ( B 4)
Formula (B 4) is the finite-dimensional analogue of variational principle (B 2). However, mathematical complexities in perturbation theory occur for degenerate eigenvalues λ (0) . Indeed, we have shown in the main text that many acoustic modes are degenerate, that is, have the same angular frequency even in spheres or spheroids for different azimuthal wavenumbers |m|. The degeneracy of the differential operator can be often avoided in principle (B 2), by solving separately the eigenvalue problem for each azimuthal wavenumber in axisymmetric containers. This strategy cannot be pursued with the polynomial description, but the difficulties can be circumvented (e.g. [START_REF] Seyranian | Sensitivity analysis of multiple eigenvalues[END_REF]. When the unperturbed eigenvalue λ (0) is degenerate, formula (B 4) ought to be modified. If λ (0) is a r-multiple root, then the solution α (0) can be written as a linear combination of the set of r-degenerate modes {α j,(0) } with j = 1, 2, . . . , r. Then, the solvability conditions for each degenerate mode α j,(0) give a linear system made of the individual equations α † j,(0)

⊤ • λ (0) C (1) + K (1) α (0) + λ (1) α † j,(0) ⊤ • 2λ (0) M + C (0) α (0) = 0.
(B 5)

This linear system admits r non-trivial solutions provided that the determinant vanishes.

The resulting condition is used to determine the r roots for the first-order correction λ (1) . We refer the reader to [START_REF] Seyranian | Sensitivity analysis of multiple eigenvalues[END_REF] for further details.

velocity perturbation, written as u * 1 = ∇Φ * 1 + ζ * 1 where ζ * 1 is a vortical hydrodynamic contribution. This is motivated by the Clebsch representation, but does not depend on it [START_REF] Bergliaffa | Wave equation for sound in fluids with vorticity[END_REF]. Then, the dimensional equations read [START_REF] Bergliaffa | Wave equation for sound in fluids with vorticity[END_REF] correctly obtained the hybridisation of acoustic waves with the Coriolis waves in an unbounded medium, in agreement with the plane-wave analysis of the primitive equations (not shown).

We may naively apply extended Goldstein equations (C 1) to compute the rotational splitting, but this approach is not consistent in rigid ellipsoids. This comes from the description of ζ * 1 . In unbounded fluids, we can expand ζ * 1 by using the Clebsch representation, such that ∇ • ζ * 1 = 0. Yet, the Clebsch representation is not as powerful as the Helmholtz-Hodge decomposition in bounded geometries. This decomposition may not (i) exist globally and (ii) satisfy the non-penetration boundary condition. From general decomposition (3.3) in rigid ellipsoids, u * 1 lies in the union V W, that is the sum of a (divergenceless) vortical space and an irrotational space with a non-zero divergence (see §3.1). Hence, we necessarily get ∇ • ζ * 1 = 0 in rigid ellipsoids. Then, equation (C 1b) becomes uncoupled (when ρ * 0 is homogeneous) and equation (C 1a) reduces to the standard acoustic equation (only valid for potential flows). This shows that equations (C 1) cannot be used in rigid ellipsoids (with homogeneous background density profiles) to compute the acoustic modes in the presence of global rotation.

Finally, extended Goldstein equations (C 1) may not be inconsistent in any bounded geometry. For instance, their relevance remains to be assessed in ellipsoids with freesurface boundary conditions. Indeed, the relevant spectral decomposition is the sum of three vector spaces U V W such that space U is spanned by divergenceless scalar potentials [START_REF] Lebovitz | The stability equations for rotating, inviscid fluids: Galerkin methods and orthogonal bases[END_REF]. This would keep the coupling between the two equations (C 1), such that rotational effects may be described by the Goldstein equations.

Figure 1 .

 1 Figure 1. Sketch of the general problem. The unit Cartesian vectors are [1x, 1y, 1z]. Left: Front view of the fluid-filled ellipsoid with (dimensional) semi-axes [a * , b * , c * ]. Right: Top view (equatorial slice). On the boundary, the velocity satisfies the non-penetration condition.

Figure 3 .Figure 4 .

 34 Figure 3. Schematic structure (non-zero entries) of the sparse matrices [M , C, K] of QEP (4.1) with MS = 0, MA = 0 and 1Ω = 1z. Polynomial degree n = 10 with the dimensions of the finite-dimensional spaces dim V [0, 10] = 495 and dim W [0, 10] = 219. Dashed (red) lines illustrate the block structures associated with the coupling between the trial and test subspaces. Upper left block: test and trial elements in V [0, 10]. Upper right block: test elements in V [0, 10] and trial elements in W [0, 10]. Lower left block: test elements in W [0, 10] and trial elements in V [0, 10]. Lower right block: test and trial elements in W [0, 10].

Figure 5 .

 5 Figure 5. Spectral convergence of the polynomial solutions in W [0, n], as a function of the truncation degree n. Relative error ∆ω between numerics ω and high-precision predictions ω i,l,m in (a) spheres (a = b = c = 1) and (b) an oblate spheroid (a = b = 1, c = 0.5). Black dashed lines show the spectral convergence ∆ω ∝ exp(-2.4n). We have used for ω i,l,m either high-precision analytical predictions[START_REF] Chang | Natural resonant frequencies of an oblate acoustical resonator[END_REF] or high-resolution computations at n = 25.

Figure 6 .

 6 Figure 6. Rotating compressible modes in spheres (a = b = c = 1). Red crosses: polynomial solutions W [0, 20] in (a) and V [0, 20] W [0, 20] in (b). Empty (black) squares: diffusive computations with COMSOL (separate m, ReS = 4 × 10 6 , νB/ν = 0.6), see details in appendix A. Dashed black lines show the upper bound ω = 2 MΩ for inertial modes. The experimental range of parameters is MΩ 2 × 10 -1 (see table 1).

Figure 7 .

 7 Figure 7. Angular frequency ω of non-rotating magneto-acoustic modes as a function of MA. (a) Plane-wave analysis for the wavenumber k = (0.1, 1.0, 0.3) ⊤ . AW: Acoustic Waves. MW: Magnetic (Alfvén) Waves. MAW: Magneto-Acoustic Waves. (b) Polynomial solutions in spheres for V [0, 15] W [0, 15]. Empty (black) circles: three-dimensional COMSOL diffusive solutions for ReS = 10 3 and RmS = 10 (#DOF = 416 456).

Figure 8 .

 8 Figure 8. Ellipsoidal splitting of the first low-frequency acoustic modes in spheroids and triaxial ellipsoids. Polynomial solutions belonging to W [0, 20]. Ellipsoidal geometries a = b = 1, c = 1-ǫ in (a) and a = 1, b = 1 + ǫ, c = 1ǫ in (b). Blue triangles: diffusionless COMSOL computations (using the built-in acoustic scalar equation).

Figure 9 .

 9 Figure 9. Accidental degeneracy of acoustic modes due to flattening in spheroids (a = b = 1 = c). Acoustic angular frequency ω as a function of the polar flattening ǫ. Comparison between polynomial solutions in W [0, 20] (squares) and finite-element computations (solid lines) with COMSOL. Vertical dashed (magenta) line shows the polar flattening ǫ = 0.05 of the ZoRo experiment (Su et al. 2020).

ωFigure 10 .

 10 Figure 10. Meridional sections of the acoustic pressure p1 ∝ -∇ • ξ of the modes near the accidental degeneracies observed in figure 9. Polynomial solutions in W [0, 20]. Oblate spheroids (a = b = 1) with ǫ = 0.05, as in the ZoRo experiment[START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF]. Each section is taken in a meridional plane containing the z-axis, at the longitude φ measured from the long x-axis. Horizontal axis shows x 2 + y 2 at the longitude φ. Colour bar shows p1, with arbitrary amplitudes and phase shifts. Solid grey lines are positive iso-contours, whereas dashed grey lines represent negative iso-contours.

Figure 11 .

 11 Figure 11. Rotational splitting of acoustic modes in spheres. Comparison between polynomial solutions belonging to W [0, 20] (black squares) and V [0, 20] W [0, 20] (red lines). Blue triangles: COMSOL diffusive computations (separate m, ReS = 4 × 10 6 , νB/ν = 0.6). The experimental range of parameters is MΩ 2 × 10 -1 (see table1).

Figure 12 .

 12 Figure 12. Scaling of the rotational shift of acoustic modes in spheres (a = b = c = 1). Ratio |δΩ|/MΩ as a function of MΩ. Dashed black lines (respectively solid red lines) are polynomial solutions belonging to W [0, 20] (respectively V [0, 20] W [0, 20]).

Figure 13 .

 13 Figure 13. Combined effects of rotation and ellipticity on acoustic modes in spheroids (a = b = c). (a) Non-rotating acoustic branches, as a function of the polar flattening ǫ. Solid (red) lines show the polynomial solutions in W [0, 20]. Dashed (black) lines show the second-order perturbation theory, obtained with a quadratic fit of a few polynomial solutions when ǫ ≪ 0.1. (b) Angular frequency ω as a function of MΩ. Red curves: non-perturbation solutions ω(ǫ, MΩ), computed from V [0, 20] W [0, 20]. Black squares: perturbation solutions (5.2), obtained by adding the frequency shift δǫ (due to ellipticity) and the rotating spherical solutions ω(0, MΩ). Frequency shift δǫ for each branch has been computed with formula (5.3), by using the non-rotating polynomial solutions shown in (a). Rotating spherical modes ω(0, MΩ) computed in V [0, 20] W [0, 20].

Figure 14 .

 14 Figure 14. Combined effects of rotation and flattening on acoustic modes in spheroids. Angular frequency ω as a function of MΩ. Red curves: non-perturbation solutions ω(ǫ, MΩ), computed in V [0, 20] W [0, 20]. Black squares: perturbation solutions (5.2), obtained by adding the frequency shift δǫ (due to ellipticity) to the rotating spherical solutions ω(0, MΩ). Frequency shift δǫ has been computed for each acoustic branch with formula (5.3), by using the non-rotating solutions in W [0, 20]. Rotating spherical modes ω(0, MΩ) computed in V [0, 20] W [0, 20].Blue triangles: COMSOL diffusive computations (separate m, ReS = 4 × 10 6 , νB/ν = 0.6). The experimental range of parameters is MΩ 2 × 10 -1 (see table1).

Figure 15 .

 15 Figure 15. Combined effects of rotation and ellipticity on acoustic modes in spheroids (ǫ = 0.05). (a) Angular frequency ω as a function of MΩ. (b) Absolute splitting δmω between the acoustic branches of azimuthal wavenumbers m and -m. Red curves: non-perturbation solutions ω(ǫ, MΩ), computed from V [0, 20] W [0, 20]. Dashed (black) curves: linear fit in Mω, modelling both first-order Coriolis effects δΩ and cross-effects δǫΩ between ellipticity and rotation. The experimental range of parameters is MΩ 2 × 10 -1 (see table1).

Figure 16 .

 16 Figure 16. Rotational splitting in the experimental range MΩ 2 × 10 -2 of low-frequency acoustic modes, for the spheroidal geometry of the ZoRo experiment (Su et al. 2020). Semi-axes a = b = 1 and c = 0.95, yielding ǫ = 0.05. Red curves are the polynomial solutions belonging to V [0, 20] W [0, 20]. Blue squares: accidental degeneracies. Black triangles: COMSOL diffusive computations (separate m, ReS = 4 × 10 6 , νB/ν = 0.6). The experimental range of parameters is MΩ 2 × 10 -1 (see table1).

ωFigure 17 .

 17 Figure 17. Meridional sections of the acoustic pressure p1 ∝ -∇ • ξ of the modes near the accidental degeneracies observed in figure 14. Polynomial solutions in V [0, 20] W [0, 20] for MΩ = 0.05. Each section is taken in a meridional plane containing the z-axis, at the longitude φ measured from the long x-axis. Horizontal axis shows x 2 + y 2 at the longitude φ. Colour bar shows p1, with arbitrary amplitudes and phase shifts. Solid grey lines are positive iso-contours, whereas dashed grey lines represent negative iso-contours.

Figure 18 .

 18 Figure 18. Numerical convergence of the 2-D solutions (separate m) computed with COMSOL, as a function of the number of DOF. Relative error ∆ω as a function of (#DOF) 1/2 , which is a proxy of the number of unknown in a given spatial direction. (a) Solutions of acoustic equation (A 1) in an oblate spheroid with a = b = 1 and c = 0.5. Dashed black lines show the algebraic convergence (here with a -4.5 slope in log-log representation). The considered modes have m = 1 for ω i,l,m ∈ {5.4446, 7.5561, 8.3503} and m = 2 for ω i,l,m = 6.7743. (b) Diffusive non-rotating acoustic modes, by solving the compressible Navier-Stokes equations with COMSOL (separatem). Computations at ReS = 4 × 10 6 and ratio νB/ν = 0.6 in a sphere (isothermal boundary at 20 • C). Solid lines show ∆ω and dashed lines the relative error on the decay rate ∆σ. The reference solutions have been computed by using quintic elements for the both the velocity and the pressure (P5-P5) and #DOF= 880644. Black triangles: elements P2-P1 (quadratic for the velocity, linear for the pressure). Red circles: elements P3-P2 (cubic for the velocity, quadratic for the pressure). Green squares: elements P4-P3 (quartic for the velocity, cubic for the pressure).

  • 2λ (0) + C (0) ξ (0) dV = -V ξ † (0) • λ (0) C (1) + K (1) ξ (0) dV. (B 2)From principle (B 2), we recover formula (3.332) in[START_REF] Aerts | Asteroseismology[END_REF] by considering a spatially uniform density profile and perturbations of non-rotating solutions (i.e. by setting 1 Ω = 0 for M S = 0). Principle (B 2) can also be formulated as a finite-dimensional problem within the framework of the spectral polynomial method. General problem (2.13) is recast as a finite-dimensional QEP with the real-valued matrices [M , C, K]. Then, formal expansion (B 1) yields the finite-dimensional quantities [λ, α, C, K] ≃ λ (0) , α (0) , C (0) , K (0) + M S λ (1) , α (1) , C (1) , K (1) . (B 3)

  dt * = ∂/∂t * + (U * 0 • ∇) the material derivative along the background flow.

Table 1 .

 1 Typical values of the dimensionless numbers for laboratory experiments in spherical-like domains. Maryland experiment: Kelley et al. (2007), Rieutord et al. (2012), Triana et al. (2014) and Zimmerman et al. (2014). Derviche Tourneur Sodium (DTS)

experiment: Figueroa et al. (2013) and Tigrine et al. (2019). ZoRo experiment: Su et al. (2020). Marseille ellipsoidal experiments: Grannan et al. (2016), Lemasquerier et al. (2017) and Le Reun et al. (2019).

  Any truncated expansion (3.16) is not an exact solution of problem (3.1). The best examples are acoustic modes, described by ∇Φ w ∈ W with Φ w solution of the scalar Helmholtz equation. They are expressed in terms of separable transcendental functions, that do not admit exact finite polynomial expressions

Table 2 .

 2 Angular frequency ω of acoustic modes in spheroids (a = b). Comparison between analytical predictions ω i,l,m (from Chang 1972) and polynomial solutions belonging to W [0, n]

† Email address for correspondence: j.n.vidal@leeds.ac.uk

Acknowledgements

JV was partly funded by STFC Grant ST/R00059X/1. SS and DC were supported by ANR-13-BS06-0010 (TuDy). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). JV acknowledges Dr B. Valette for illuminating discussions on the modes' properties. The authors also acknowledge Dr P. Cardin and Dr H.-C. Nataf for fruitful discussions on acoustics. The authors thank the three reviewers. Their comments substantially improved the quality of the paper.

Declaration of interests

The authors report no conflict of interest.

Appendix A. Finite-element computations

For the numerical benchmarks, we have used the finite-element commercial code COMSOL. The ellipsoidal domain is discretised by using an unstructured mesh with tetrahedral elements. The mesh elements are the standard Lagrange elements for the pressure and the velocity. We have solved in practice different versions of the governing equations (2.5), depending on the presence of rotational or magnetic effects.

A.1. Non-magnetic configurations (M A = 0)

In non-rotating and non-magnetic cases (M Ω = M S = M A = 0), the governing equations reduce to the standard (dimensionless) acoustic equation for the pressure

(A 1) Equation (A 1) is directly solved in COMSOL with a built-in eigenvalue solver, by using either three-dimensional (3-D) computations or a Fourier decomposition exp(imφ) in the azimuthal direction to separate the different m (in axisymmetric containers only). In the latter case, we have extended the built-in implementation of (A 1) in COMSOL to account for non-axisymmetric modes m = 0. We have used cubic Lagrange elements for the computations. The numerical convergence is shown in figure 18(a), as a function of the number of DOF (degrees of freedom). We recover the (slow) algebraic convergence of finite-element computations, here with a -4.5 slope in log-log representation. The convergence is much weaker than the exponential convergence of the spectral polynomial method (compared to figure 5). We also obtain a -4.5 slope for the convergence of the 3-D computations, but the relative error is larger than for the 2-D computations (not shown). Acoustic equation (A 1) is only valid when M Ω = M S = M A = 0. In any other case, we must go back to the primitive equations, which necessarily include diffusion (e.g. viscosity) for the convergence of the numerical results. This is due to the piecewise continuous polynomials used in the finite-element computations, which introduce numerical diffusion polluting artificially the results. This is strongly different from our global polynomial spectral method. The latter relies on infinitely differentiable polynomials, probing accurately the diffusionless regime. We have solved the non-magnetic equations (2.5a)-(2.5b) with the built-in adiabatic (and viscous) formulation of COMSOL. We have also extended the COMSOL formulation to solve separately for non-axisymmetric modes m = 0 (with a Fourier decomposition). On the rigid boundary, we prescribe the no-slip condition u 1 = 0. We show in figure 18(b) the typical numerical convergence of the diffusive eigenvalue λ = σ + iω, with σ < 0 the diffusive decay rate and ω ∈ R the angular frequency, as a function of the number of DOF. We have fixed the diffusion at Re S = 4 × 10 6 and ν B /ν = 0.6 (for non-rotating computations). Note that the diffusive frequencies ω are always smaller than the diffusionless ones, as expected theoretically (e.g. [START_REF] Moldover | Gas-filled spherical resonators: Theory and experiment[END_REF]. We also obtain an algebraic convergence. The angular frequency ω converges faster than the decay rate σ towards the expected diffusive solutions (determined by high-resolution computations). We have also varied the polynomial order. High polynomial orders are more accurate for σ but, surprisingly, less accurate for ω. Hence, we must choose the polynomial order to have the desired convergence of both ω Then, the (dimensionless) induction equation reads

To solve equation (A 3), the first magnetic condition in (A 2) directly translates into Λ 1 × 1 n = 0 on the boundary, whereas the second is naturally satisfied by the chosen Nédélec's elements on the boundary. When a magnetic field was considered, we have used quadratic Lagrange P2-P2 elements (quadratic for the velocity and the pressure), combined with cubic Nédélec's elements for the magnetic field. Finally, we have also neglected the bulk viscosity (ν B /ν = 0) in our magnetic computations.

Appendix B. Additional details on MAV

B.1. Magnetic shift

We quantify the frequency shift δ B due to magnetic effects. We have already observed in figure 7 that the magnetic field lifts the azimuthal degeneracy of the spherical acoustic modes. The magnetic shift scales in M 2 A (not shown), as expected from perturbation theory (e.g. [START_REF] Gough | The effect of rotation and a buried magnetic field on stellar oscillations[END_REF] for different magnetic fields). For our magnetic field, the numerical prefactor in the scaling law δ B ∝ M 2

A is of order unity (at least for the first low-frequency acoustic modes). This yields the raw estimate δ B ≃ M 2 A . Typical experimental values are M A 5 × 10 -3 (see table 1), leading to the upper bound δ B 10 -5 (in dimensionless units).

In dimensional units, taking the ZoRo apparatus as a reference [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF], the typical frequency resolution is 0.1 -1 Hz (for modes with a dimensional frequency f * = ω * /(2π) of a few kHz). With ambient air (C * 0 ≃ 343 m/s), the frequency shift would be as small as 10 -3 Hz in dimensional units within the ZoRo geometry (radius a * = 0.2 m). The situation is not very favourable in liquid metals (in which the speed of sound is higher). Therefore, we conclude that the frequency shift δ B due to magnetic fields is entirely negligible for experimental conditions.

B.2. Splitting due to a background velocity

We provide details on the method to calculate the effects of a background flow on the acoustic eigenmodes when M S = 0 (and M A = 0 since magnetic effects are negligible). We could solve directly problem (2.13) with the polynomial description. However, the presence of a background flow U 0 (with an arbitrary spatial complexity) may reduce the numerical convergence of the polynomial solutions. Indeed, the terms involving the background flow are responsible for a cascade on higher-order polynomial bases (unless U 0 is linear in the Cartesian coordinates). Moreover, flow instabilities would be certainly triggered (e.g. [START_REF] Kerswell | The instability of precessing flow[END_REF][START_REF] Kerswell | Elliptical instability[END_REF]Vidal & Cébron 2017). In experimental conditions, the large-scale background velocity is generally sub-sonic (M S ≪ 1), and also often smaller than solid-body rotation (i.e. Ro = M S /M Ω 1). This suggests seeking the solutions in powers of M S , to account for the effects of U 0 on the acoustic spectrum only as small perturbations. This reads in dimensionless form

Appendix C. Extended Goldstein equations in rigid ellipsoids

The standard acoustic equation, which is only valid for potential flows [START_REF] Pierce | Wave equation for sound in fluids with unsteady inhomogeneous flow[END_REF], cannot be used to determine any rotational effect. Acousticians have attempted to extend the acoustic equation, to account for vortical perturbations [START_REF] Goldstein | Unsteady vortical and entropic distortions of potential flows round arbitrary obstacles[END_REF]) and vortical mean flows U * 0 [START_REF] Bergliaffa | Wave equation for sound in fluids with vorticity[END_REF]). The resulting equations bear the name of (extended) Goldstein equations [START_REF] Bensalah | Well-posedness of a generalized time-harmonic transport equation for acoustics in flow[END_REF]. They have for unknown the