David Cébron

Jérémie Vidal

Nathanaël Schaeffer

Antonin Borderies

Alban Sauret

Mean zonal flows induced by weak mechanical forcings in rotating spheroids

Keywords: rotating flows, boundary layer, topographic effects

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Physical context

Global rotation tends to sustain two-dimensional mean flows that are almost invariant along the rotation axis in rapidly rotating systems. These mean flows are indeed often obtained in various models of rotating turbulence (e.g. [START_REF] Guervilly | Large-scale vortices in rapidly rotating Rayleigh-Bénard convection[END_REF][START_REF] Godeferd | Structure and dynamics of rotating turbulence: a review of recent experimental and numerical results[END_REF] and planetary core flows (e.g. [START_REF] Aubert | Steady zonal flows in spherical shell dynamos[END_REF][START_REF] Schaeffer | Turbulent geodynamo simulations: a leap towards Earth's core[END_REF][START_REF] Monville | Rotating double-diffusive convection in stably stratified planetary cores[END_REF]). In the latter context, they are believed to play an important role in the exchange of angular momentum between liquid layers and surrounding solid domains (e.g. [START_REF] Roberts | On the theory of core-mantle coupling[END_REF], which drives the long-term dynamical evolution of planetary bodies. Moreover, mean flows could be unstable in the rapidly rotating regime (e.g. [START_REF] Sauret | Tide-driven shear instability in planetary liquid cores[END_REF][START_REF] Favier | Non-linear evolution of tidally forced inertial waves in rotating fluid bodies[END_REF], which could sustain space-filling turbulence and mixing. Therefore, understanding the formation of mean flows is essential to model the fluid dynamics of many rapidly rotating systems.

A commonly observed feature of geostrophic flows is that they are spontaneously generated by nonlinear effects, for instance involving small-scale eddies (e.g. [START_REF] Christensen | Zonal flow driven by strongly supercritical convection in rotating spherical shells[END_REF][START_REF] Aubert | Observations of zonal flow created by potential vorticity mixing in a rotating fluid[END_REF] or waves. Rapidly rotating fluids are indeed characterised by the ubiquitous presence of inertial waves (e.g. [START_REF] Zhang | Theory and modeling of rotating fluids: convection, inertial waves and precession[END_REF], whose restoring force is the Coriolis force. However, [START_REF] Greenspan | On the non-linear interaction of inertial modes[END_REF] demonstrated that inviscid nonlinear interactions of inertial waves do not produce significant geostrophic flows in the rapidly rotating regime. The combination of some nonlinear interactions and viscous effects is thus essential to generate mean geostrophic flows, and various wave-induced mechanisms have been explored. Local wave interactions in the weakly viscous interior could transfer energy from the inertial waves to the geostrophic flows, either through wave-wave interactions (e.g. [START_REF] Newell | Rossby wave packet interactions[END_REF][START_REF] Smith | Transfer of energy to two-dimensional large scales in forced, rotating three-dimensional turbulence[END_REF] or wave-induced secondary instabilities (e.g. [START_REF] Kerswell | Secondary instabilities in rapidly rotating fluids: inertial wave breakdown[END_REF][START_REF] Brunet | Shortcut to geostrophy in wave-driven rotating turbulence: the quartetic instability[END_REF]. The aforementioned mechanisms have been explored in Cartesian or cylindrical geometries for computational simplicity. In these previous studies, the container depth does not vary in the direction perpendicular to the rotation axis. However, this so-called beta effect is known to be important for planetary configurations [START_REF] Busse | Thermal instabilities in rapidly rotating systems[END_REF], and also strongly modifies the geostrophic flows [START_REF] Greenspan | The theory of rotating fluids[END_REF]). Thus, although these local mechanisms are certainly generic, the geostrophic flows investigated in these studies are not directly relevant for (large-scale) planetary core flows.

Another mechanism, which is relevant for planetary applications, has been proposed by Busse (1968b). Most planetary fluid bodies are subject to mechanical forcings (e.g. librations, precession, or tides) because of the presence of orbital companions. Mechanical forcings have received a renewed interest in fluid mechanics, because of their non-negligible contribution in the internal fluid dynamics of planetary bodies (e.g. Le [START_REF] Bars | Flows driven by libration, precession, and tides[END_REF]. They are indeed responsible for differential motions of the rigid boundary with respect to the fluid. These motions can be transmitted to the bulk by viscous coupling, generating inertial waves (e.g. [START_REF] Aldridge | Axisymmetric inertial oscillations of a fluid in a rotating spherical container[END_REF]Noir et al. 2001a;Sauret et al. 2013), and mean geostrophic flows resulting from nonlinear interactions of the flows within the Ekman boundary layer (as considered in Busse 1968b). The latter mechanism has been then confirmed experimentally and numerically for various mechanical forcings (e.g. Noir et al. 2001b[START_REF] Noir | Experimental study of libration-driven zonal flows in non-axisymmetric containers[END_REF][START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF].

Motivations

Rotating flows are usually characterised by the Ekman number E, which compares viscous to rotational effects. As outlined in Busse (1968b[START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF], mechanical forcings of typical amplitude O( ) can induce a mean zonal flow in the bulk of typical amplitude O( 2 ), which is independent of E in the regime E → 0 (as shown below). Hence, this mechanism gives a non-zero bulk flow driven by viscous effects that survives in the planetary regime of vanishing viscosity E 1. However, most studies about mean zonal flows in spherical-like domains have employed laboratory experiments or direct numerical simulations (DNS) with moderate values E 10 -6 , whereas rapidly rotating planetary flows are characterised by much smaller values (typically E = 10 -15 -10 -12 ). Thus, since viscous effects are overestimated in experimental or numerical works, an analytical study is directly relevant to obtain rigorous results about zonal flows in the planetary regime.

Only a few theoretical studies have hitherto investigated mean zonal flows driven by mechanical forcings. The case of a rotating cylindrical tank subject to longitudinal librations has been recently revisited analytically [START_REF] Sauret | Mean zonal flow generated by azimuthal harmonic forcing in a rotating cylinder[END_REF], and the theory has been convincingly compared with experiments [START_REF] Wang | Cylindrical tank of fluid oscillating about a state of steady rotation[END_REF]) and simulations (Sauret et al. Experiments [START_REF] Sauret | Experimental and numerical study of mean zonal flows generated by librations of a rotating spherical cavity[END_REF] DNS [START_REF] Sauret | Experimental and numerical study of mean zonal flows generated by librations of a rotating spherical cavity[END_REF] Theory [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF] Theory [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF] 0 0.2 0.4 0.6 0.8 1

s -0.2 0 0.2 V φ ( ×10 3
)

Experiments [START_REF] Malkus | Precession of the Earth as the cause of geomagnetism: Experiments lend support to the proposal that precessional torques drive the Earth's dynamo[END_REF] DNS (Noir et al., 2001) Theory (Busse, 1968) (a) (b)

Figure 1. Comparison of the azimuthal component of the time-averaged zonal flow V φ , as a function of cylindrical radius s in the equatorial plane z0 = 0, between theory (dashed curves) and experiments/numerics (solid curves). (a) Libration-driven zonal flows with dimensionless angular frequency ω = 0.1 in a full sphere, extracted from figure 3b and figure 4 in [START_REF] Sauret | Experimental and numerical study of mean zonal flows generated by librations of a rotating spherical cavity[END_REF]. Experiment with E = 1.15 × 10 -5 , and = 0.08. DNS with E = 5 × 10 -5 , and = 0.2. The vertical axis has been normalised by 2 , where is the dimensionless forcing amplitude. (b) Precession-driven zonal flows extracted from figure 8 in Noir et al. (2001b).

2012

). However, a successful validation of mechanically driven zonal flows is generally missing in spherical and ellipsoidal geometries when E is vanishingly small. For instance, considering longitudinal librations in spheres, two different results have been obtained for low libration frequencies [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF][START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF]. As shown in figure 1(a), they can both explain the experimental results of [START_REF] Sauret | Experimental and numerical study of mean zonal flows generated by librations of a rotating spherical cavity[END_REF]. Concerning precession, experiments [START_REF] Malkus | Precession of the Earth as the cause of geomagnetism: Experiments lend support to the proposal that precessional torques drive the Earth's dynamo[END_REF]) or numerical simulations (Noir et al. 2001b) have never properly validated the theory of Busse (1968b), as illustrated in figure 1(b).

Similarly, the theoretical zonal flows driven by tides [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF] do not agree with experimental findings (as we will show below). Consequently, theoretical predictions remain to be thoroughly validated before they can be extrapolated for planets. Finally, singularities have been found in the boundary-layer calculation due to the presence of the critical latitudes, where the flows should be smoothed out by additional viscous effects (e.g. [START_REF] Kerswell | On the internal shear layers spawned by the critical regions in oscillatory ekman boundary layers[END_REF][START_REF] Kida | Steady flow in a rapidly rotating sphere with weak precession[END_REF]) not taken into account in the theory (as in Busse 1968b;Sauret et al. 2013). Around these locations, the mean flows are known to take the form of narrow geostrophic shear layers aligned with the axis of rotation (e.g. [START_REF] Calkins | Axisymmetric simulations of libration-driven fluid dynamics in a spherical shell geometry[END_REF]. The variations of the geostrophic shear amplitude with the Ekman number are however still disputed (Noir et al. 2001b;[START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF], such that planetary extrapolations remain speculative. Thus, targeted DNS in the regime E 1 are also worth performing to explore the behaviour of the geostrophic shear layers.

Solving the full mathematical problem of mechanically driven flows is complex, but analytical progress can be made for planetary parameters (as undertaken in Busse 1968b[START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF]. Since planetary interiors are characterised by small forcing amplitudes 1 and small viscous effects E 1, we will employ asymptotic theory in and E. Moreover, our zonal flow calculation will also assume that the spin-up time scale of the fluid [START_REF] Greenspan | The theory of rotating fluids[END_REF]) is much longer than the characteristic time scale of the mechanical forcing (in the fluid rotating frame), such that no global spin-up of the fluid will occur during the dynamics. Finally, following previous works on mean zonal flows (e.g. [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF][START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF], we neglect in the theory viscous effects at the critical latitudes (associated with internal shear layers) and our theoretical bulk basic flow is taken as a solid-body rotation (for its spatial dependency). Thus, we assume that no inertial mode is excited by the forcing on top of this basic flow (which is exact if the forcing frequency is larger than twice the mean fluid rotation rate, e.g. [START_REF] Greenspan | The theory of rotating fluids[END_REF]). Using DNS, where these effects are fully taken into account, we will revisit the proposed associated scaling laws for planetary extrapolations (e.g. Noir et al. 2001b;[START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF]. The paper is organised as follows. We introduce the problem and the methods in §2. We describe the asymptotic weakly nonlinear analysis in §3, and present the theoretical and numerical results in §4. We discuss the results in §5, and we finally conclude the paper in §6.

Description of the problem and methods

Mathematical description

We consider an incompressible and homogeneous Newtonian fluid of kinematic viscosity ν and density ρ, enclosed in a spheroidal container of semi-axis length r pol along the revolution axis, while the other one is noted r eq (the spheroid is oblate when r eq > r pol , and prolate when r eq < r pol ). We introduce the Cartesian basis vector ( x I , y I , z I ) of the inertial frame, whose origin O is the centre of the spheroidal container. In the following, we work in a frame of reference where the spheroidal shape of the container boundary is stationary. We use a Cartesian basis vectors ( x R , y R , z R ) where z R is aligned with the spheroid revolution axis, as illustrated in figure 2. The rotation vector of this reference frame, denoted Ω * c (t) in the following, is along z I in the absence of perturbation. In this reference frame, the velocity V * satisfies the no-slip boundary conditions (BC) on Σ

V * • n| Σ = 0, V * × n| Σ = V * Σ , (2.1a,b)
where n is the unit vector normal to the boundary, and V * Σ is mainly a solid-body rotation at Ω * 0 around z I , possibly perturbed by a small flow v * Σ . We thus consider

V * Σ = Ω * 0 z R × r * + v * Σ , (2.2)
where r * is the position vector, and v * Σ is an imposed tangential velocity related to the considered mechanical forcing. One can choose a frame of reference where the spheroidal boundary is steady, that is with V * Σ = 0. This frame is referred as the wall frame, or the mantle frame in planetary sciences. We denote the associated Cartesian basis vectors ( x M , y M , z M ) centred on O, where z M is along the spheroid revolution axis.

We work below using dimensionless units, denoting the dimensionless variables without the superscript * for the sake of clarity. We choose r eq as the length scale, and

|Ω * s | -1 = (|Ω * 0 | + |Ω * c |) -1
as the time scale, where Ω * c is the time average of Ω * c . In the reference frame rotating at Ω c , the dimensionless fluid velocity V is governed by

∂ t V + (V • ∇) V + 2 Ω c × V + Ωc × r = -∇Π + E ∇ 2 V , ∇ • V = 0, (2.3a,b)
where Π is the reduced pressure (taking into account the centrifugal effects), Ωc = d t Ω c is the time derivative of Ω c , ∂ t V is the partial time derivative of V , and E = ν/(Ω * s r 2 eq ) is the (dimensionless) Ekman number.

Mechanical forcings

We note U the inviscid bulk flow driven by the forcing. In the absence of any mechanical forcings, that is with V Σ = 0 and Ωc = 0, U reduces to a solid-body rotation with the angular velocity Ω 0 z R . However, the latter flow is perturbed by weak harmonic perturbations generated by mechanical forcings. The general framework described in this study allows us to consider various mechanical forcings that are described below.

(i) Multipolar tidal-like forcing corresponds to Ω c = Ω c z R (with here z R = z I ), and

v Σ = s q cos(mφ) cos(ωt) z R × r, (2.4)
with the azimuthal angle φ with respect to z R (see figure 2), the forcing amplitude , and the azimuthal wavenumber m of the spatial deformation (this boundary velocity has also been considered by [START_REF] Greenspan | The theory of rotating fluids[END_REF], see e.g. his equation 2.14.2). In our calculation, q is taken as an independent parameter but, for regularity along the rotation axis [START_REF] Lewis | Physical constraints on the coefficients of Fourier expansions in cylindrical coordinates[END_REF], and to consider multipolar flows (e.g. [START_REF] Cébron | Libration-driven multipolar instabilities[END_REF][START_REF] Sauret | Mean zonal flow generated by azimuthal harmonic forcing in a rotating cylinder[END_REF], we must consider q = |m -1|. In expression (2.4), the case m = 2, Ω c = ω = 0 has been considered in [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF], which is extended here to account for both multipolar deformations and oscillations at the frequency ω (e.g. [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF].

(ii) Longitudinal librations are investigated with Ω s = 1 and z R = z I . Introducing the forcing amplitude , rotating spheroids can be studied, in an equivalent way, either (i) in the mantle frame of reference (e.g. [START_REF] Favier | Generation and maintenance of bulk turbulence by libration-driven elliptical instability[END_REF] with Ω 0 = 0, Ω c = [1+ cos(ωt)] z R and V Σ = 0, (ii) in the mean rotating frame of reference (e.g. [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF] with Ω 0 = 0, Ω c = z R and V Σ = cos(ωt) z R × r, or (iii) in the inertial frame of reference with Ω 0 = 1, Ω c = 0 and V Σ = [1 + cos(ωt)] z R × r . Note that the case (ii) can actually be recovered with the particular case m = Ω 0 = 0 of the multipolar tidal-like forcing.

(iii) Latitudinal librations are modelled with Ω 0 = 0, and the general case of a rigid spheroidal container can only be studied in the mantle frame, where the Σ is stationary and V Σ = 0. The corresponding forcing in this frame is (see in [START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF])

Ω c = Θ x M + sin(Θ) y M + cos(Θ) z M , (2.5a) Ωc = Θ x M + Θ cos(Θ) y M -Θ sin(Θ) z M , (2.5b)
where Θ = ( /ω) sin(ωt) is the instantaneous libration angle, and where is the forcing amplitude. In the limit 1 considered for the analytical calculations performed in this work, these expressions read (at first order in )

Ω c = [cos(ωt) x M + ω -1 sin(ωt) y M ] + z M , (2.6a) Ωc = [-ω sin(ωt) x M + cos(ωt) y M ],
(2.6b)

Note that the particular case of the sphere can also be studied analytically and numerically in the mean rotating frame with

Ω c = z R = z I and V Σ = cos(ωt) x R × r.
(iv) Precession can be considered in the precession frame (e.g.

V Σ = Ω 0 z R × r with Ω 0 = (1 + P o ) -1 ,
and

Ω c = Ω 0 P o [sin(α) x R + cos(α) z R ], (2.7)
where α is the precession angle, and P o is the Poincaré number (ratio of the precession and the boundary rotation rates). The associated bulk flow U is then mainly a tilted (stationary) solid-body rotation U = ω f × r. For weak precession forcing, ω f ≈ (1 + P o ) -1 z R at the order , where characterises the small misalignment of ω f and z R (i.e. the forcing amplitude, see further details in the seminal work of Busse 1968b).

Numerical modelling

We integrate equations (2.3) using two open-source codes. Equations in spherical geometries are solved using the parallel pseudo-spectral code xshells (e.g. [START_REF] Schaeffer | Turbulent geodynamo simulations: a leap towards Earth's core[END_REF], based on a poloidal-toroidal decomposition of the velocity field onto spherical harmonics of degree l l max and azimuthal wavenumber m m max using the shtns library [START_REF] Schaeffer | Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations[END_REF], and second-order finite differences with N r points are used in the radial direction. The code has been validated for full-sphere computations, including flows crossing the origin [START_REF] Marti | Full sphere hydrodynamic and dynamo benchmarks[END_REF], and details about the implementation at the centre are given in appendix A. To solve the dynamical equations, the code can use several semi-implicit time-stepping schemes, which treat the diffusive terms implicitly and the other ones explicitly. Most of the simulations have been performed using the accurate semi-implicit backward difference formula of order 3 (SBDF3, see [START_REF] Ascher | Implicit-explicit methods for timedependent partial differential equations[END_REF]. The typical spatial resolution at E = 10 -7 is N r = 576, l max = 159, m max = 5.

In spheroidal geometries, we solve the nonlinear equations in their weak variational form using the spectral-element code Nek5000 (e.g. [START_REF] Fischer | Simulation of high-Reynolds number vascular flows[END_REF], which combines the geometrical flexibility of finite element methods with the accuracy of spectral methods. The computational domain is made of E = 3584 non-overlapping hexahedral elements in coreless geometries (or E = 3840 in spheroidal shells, see below). Moreover, the velocity (and pressure) is represented within each element as Lagrange polynomials of order N = 13 (respectively, N -2) on the Gauss-Lobatto-Legendre (Gauss-Legendre) points. Temporal discretisation is accomplished by a third-order method, based on an adaptive and semi-implicit scheme in which the nonlinear and Coriolis terms are treated explicitly, and the remaining linear terms are treated implicitly. We have checked the numerical accuracy in targeted simulations by varying the polynomial order from N = 13 to N = 15, and found that the resolution of the Ekman boundary layers is appropriate with at least ten grid points within the layer.

Extracting the mean zonal flows from DNS

In the planetary limit 1 considered in this work, the forced flow is mainly a solid-body rotation in the bulk for all the aforementioned forcings. Consequently, when E → 0, the mean zonal flows tend to geostrophic flows, which are invariant along the fluid rotation axis and are established on the dimensionless spin-up time scale E -1/2 (Greenspan 1968). Thus, for every DNS, we have simulated the dynamics over several spin-up time scale, ensuring that the mean zonal flows are well established. We have also used typical time steps dt = 10 -3 -10 -2 , which were sufficient to integrate the dynamics.

The mean zonal flow is computed from the three-dimensional flow V by considering the cylindrical radial variation of V φ (s, z = z 0 ) in the horizontal plane z = z 0 , where X and X are the time and azimuthal averages of the quantity X (respectively). With xshells, the mean zonal flow is computed from the time-averaged m = 0 component of the toroidal scalar in the plane z = z 0 . In spheroids, the azimuthal component V φ of the mean zonal flow is estimated as [START_REF] Favier | Generation and maintenance of bulk turbulence by libration-driven elliptical instability[END_REF])

V φ (s, z = z 0 ) 1 N |z-z0| zmax 0 φ 2π s-ds<s<s+ds V • φ, (2.8)
where N is the total number of grid points used to evaluate expression (2.8). The Nek5000 DNS have been performed at E 10 -6 , contrary to the xshells DNS performed at E 10 -7 . Hence, the Nek5000 DNS are more influenced by Ekman pumping when approaching the boundary. To properly estimate the geostrophic components, we z-average the flows over the vertical positions |zz 0 | z max . In the xshells DNS, the mean zonal flow is defined as the value of the m = 0 azimuthal velocity in the plane z = z 0 . To be consistent, we have here considered z max = 0.1. We have also checked that the mean flow computations are unchanged when using z max /r pol 0.4 -0.5. Moreover, the approximate number of points in each direction in the Nek5000 DNS is here E 1/3 N ≈ 200. We have thus averaged the azimuthal component over one hundred different shells along the cylindrical radius s, and over the vertical positions |z -z 0 | 0.1. We show in figure 3 the mean zonal flows computed from DNS in spheres with E = 10 -4 and = 10 -2 , with Nek5000 in the mantle frame of reference and with xshells in the frame rotating at z R . We find a very good agreement between the two codes, which validates our procedure to compute the mean zonal flows.

Perturbation approach

Current DNS cannot be performed at the extremely small values of E reached in planetary liquid cores. Thus, we solve analytically equations (2.3) to gain physical insights into the asymptotic regime E 1. We assume that there is no significant shear in the interior and use viscous boundary-layer theory (BLT, e.g. [START_REF] Greenspan | The theory of rotating fluids[END_REF]) to write [V , Π] = [U , P ] + [u, p], where U describes the interior flow for which viscous effects can be neglected, and a boundary-layer flow u. The latter contribution takes into account the viscous effects near the outer boundary, and decays exponentially towards the interior of the container. The governing equations, obtained from (2.3), are in the limit E 1

∂ t U + (U • ∇) U + 2 Ω c × U + Ωc × r = -∇P, (2.9a) ∂ t u + (u • ∇) U + (u + U ) • ∇u + 2 Ω c × u = -∇p + E ∇ 2 u, (2.9b)
Within the boundary layer, we introduce the stretched coordinate ζ = (r |Σr) • n/E 1/2 , where r |Σ is the position vector on the boundary Σ. We also assume that the field gradients along the boundary are negligible compared to the gradients normal to the boundary, that is n • ∇ -E 1/2 ∂ ζ (e.g. [START_REF] Greenspan | The theory of rotating fluids[END_REF]). Then, the mass conservation equation reduces to its usual boundary layer approximation (see p. 25 in [START_REF] Greenspan | The theory of rotating fluids[END_REF])

∇ • u -E 1/2 ∂ ζ (u • n) + n • ∇ × ( n × u) = 0.
(2.10)

To solve the BLT equations, we use asymptotic theory with the small forcing amplitude 1. (2.11a,b) where the time average of Ω c is Ω c = Ω 0 c z R (since the perturbations are harmonic). The BC (2.1) imposes then V k Σ • n = 0 at every order k . To perform the boundary-layer and perturbation calculations, we also expand [V , Π] in double power series involving the asymptotic parameters 1 and E 1. Note that we formally neglect the possible critical latitudes (e.g. [START_REF] Kerswell | On the internal shear layers spawned by the critical regions in oscillatory ekman boundary layers[END_REF], although they can modify the mean zonal flows (as previously found in cylinders by [START_REF] Sauret | Fluid flows in a librating cylinder[END_REF], see also below). Since the Ekman layer scales as E 1/2 outside the critical latitudes, we use the double power series expansions for all our unknowns

Noting V 0 Σ = Ω 0 z R × r, we write V Σ = V 0 Σ + V 1 Σ + 2 V 2 Σ + ..., Ω c = Ω 0 c z R + Ω 1 c + 2 Ω 2 c + ...,
[U , u] = +∞ i,j=0 E 1 2 i j U j i , u j i , [P, p] = +∞ i,j=0
E 1 2 i j P j i , p j i .

(2.12a,b)

The substitution of equations (2.12) in the governing equations (2.9)-(2.10) leads to a sequence of equations for the interior and boundary-layer flows. We anticipate that our flows may also a priori vary slowly on the time scale τ = O(E -1/2 ) and, thus, we also expand below the time in powers of E -1/2 (as e.g. done when calculating the Ekman layer damping of inertial modes, see [START_REF] Greenspan | The theory of rotating fluids[END_REF]).

Governing equations

At the leading order 0 E 0 , a natural solution of the interior zeroth-order equation is U 0 0 = Ω 0 z R × r, which satisfies the BC since Ω 0 c is constant and v Σ = 0 at this order. Since U 0 0 verifies the BC, we obtain that u 0 0 = 0. Noting that, for an arbitrary velocity field v, we have (2.13) the first-order interior flow equations at the order E 0 are then

(U 0 0 • ∇) v + (v • ∇) U 0 0 = 2 Ω 0 z R × v + Ω 0 ∂ φ v,
(∂ t + Ω 0 ∂ φ ) U 1 0 + 2 Ω 0 s z R × U 1 0 + 2 Ω 1 c × U 0 0 + Ω1 c × r = -∇P 1 0 , (2.14)
together with the divergenceless condition ∇ • U 1 0 = 0 and Ω 0 s = Ω 0 c + Ω 0 , where the BC is U 1 0 • n Σ = 0. Considering now the boundary-layer flows, we first integrate equation (2.10) using the BC u k 0 (ζ) → 0 when ζ → +∞ and, since the first term is of order E -1/2 , we obtain the zeroth order condition for the boundary-layer flow u k 0 • n = 0 at every order k inside the boundary layer (since

V k Σ • n = v k Σ • n = 0).
Then, the boundary-layer equation at the order E 0 is (2.15a,b) with the linear operator Lu

Lu 1 0 = -∂ ζ p 1 1 n, u 1 0 + U 1 0 = V 1 Σ on Σ,
1 0 = -(∂ t + Ω 0 ∂ φ ) u 1 0 -2 Ω 0 s z R × u 1 0 + ∂ 2 ζζ u 1 0 . At the order E 1/2 , the so-called Ekman circulation U 1 1 is governed by (∂ t + Ω 0 ∂ φ )U 1 1 + ∂ τ U 1 0 + 2 Ω 0 s z R × U 1 1 = -∇P 1 1 , u 1 1 + U 1 1 = V 1 Σ on Σ, (2.16)
and the divergenceless condition ∇ • U 1 1 = 0, where we have anticipated that U 0 1 , and thus u 0 1 , can be set to zero without loss of generality (the flows are forced by the forcing of amplitude , and we will see that all our equations can be verified with the solution U 0 1 = u 0 1 = 0). At this order, the mass conservation imposes (2.17) which allows us to obtain easily the Ekman pumping u 1 1 = u 1 1 n from u 1 0 . At next order 2 E 0 , the bulk flow is divergenceless ∇ • U 2 0 = 0 and given by

∂ ζ (u 1 1 • n) -n • ∇ × ( n × u 1 0 ) = 0,
(∂ t + Ω 0 ∂ φ ) U 2 0 + 2 Ω 0 s z R × U 2 0 = -∇P 2 0 -(U 1 0 • ∇) U 1 0 -2Ω 1 c × U 1 0 -2Ω 2 c × U 0 0 - Ω2 c × r, (2.18)
with the BC U 2 0 • n = 0 on Σ. The boundary-layer equations are

Lu 2 0 = -n ∂ ζ p 2 1 + (u 1 0 • ∇)u 1 0 -(u 1 1 + U 1 1 ) • n ∂ ζ u 1 0 + 2Ω 1 c × u 1 0 +(U 1 0 • ∇)u 1 0 + (u 1 0 • ∇)U 1 0 , (2.19a)
u 2 0 + U 2 0 = V 2 Σ on Σ. (2.19b)
In the equations (2.18)-(2.19) governing the order 2 E 0 , U 1 1 and u 1 1 only appear via their normal components in the boundary layer. Then, equation (2.17) gives directly u 1 1 • n, and thus

U 1 1 • n = -u 1 1 • n| ζ=0 via the BC (u 1 1 + U 1 1 )
• n = 0 on Σ. Finally, the interior flow U k 0 at every order k will be decomposed as U k 0 = P U k 0 + H U k 0 , where P U k 0 is the particular solution forced by non-homogeneous terms, and H U k 0 is the solution of the homogeneous part of the equation that is required to satisfy the BC for the total flow. We will actually see that H U 1 0 = 0 in certain cases (e.g. in the fast libration limit ω E 1/2 ). Moreover, in all the cases considered here, we will show that the theory gives U 2 0 = P U 2 0 + H U 2 0 as an azimuthal flow, which provides the leading-order azimuthal component of V φ . We will thus compare the values of V φ obtained from DNS with the theoretical values of U 2 0 .

Asymptotic analysis

In this section, we aim at calculating the steady axisymmetric component U 2 0 of the interior flow U 2 0 , which requires the full mathematical expressions of U 1 0 , u 1 0 and u 2 0 . To solve the corresponding equations, we employ the spheroidal orthogonal coordinates (q 1 , q 2 , φ), associated with the orthogonal normal unit basis ( q 1 , q 2 , φ) where φ is the usual azimuthal unit vector. We introduce the change of variables x R = aT (q1) sin q 2 cos φ, y R = aT (q1) sin q 2 sin φ, z R = aT (q1) cos q 2 , (3.1a-c) where (q 1 , q 2 , φ) are spheroidal coordinates, a = |1 -(r pol /r eq ) 2 | 1/2 = (T (q1) 2 + T (q1) 2 ) 1/2 is the distance between the centre and the foci of the ellipse. For later use, we also define the cylindrical radius s = (x 2 R + y 2 R ) 1/2 = aT (q1) sin q 2 and the scale factors (h 1 , h 2 , h φ ) for the coordinates (q 1 , q 2 , φ) as (3.2a,b) which gives h = sinh 2 q 1 + cos 2 q 2 and h = cosh 2 q 1cos 2 q 2 for oblate and prolate spheroidal coordinates, respectively. Note that we recover the usual spherical case with q 1 → ∞, giving for instance a h → a exp(q 1 )/2 ≈ r or aT (q1) → r, with the spherical radius r (see e.g. [START_REF] Schmitt | Numerical study of a rotating fluid in a spheroidal container[END_REF]. This definition of the spheroidal coordinates allows us to encompass both oblate and prolate spheroidal coordinates in a single framework, by using respectively T (q1) = cosh q 1 when r eq > r pol , and T (q1) = sinh q 1 otherwise.

h 1 = h 2 = a T (q1) 2 + (T (q1) 2 -T (q1) 2 ) cos 2 q 2 = a h, h φ = aT (q1) sin q 2 ,
Here, we note T the derivative of T with respect to q 1 . In spheroidal coordinates, the semi-axes r eq and r pol are given by r eq = aT (Q1) and r pol = aT (Q1) , where Q 1 is the value of the radial-like coordinate q 1 at the boundary. In appendix B.1, we give various useful expressions related to the spheroidal coordinates used in this work.

First-order flows

Considering first the interior flow U 1 0 = P U 1 0 + H U 1 0 , the particular solutions P U 1 0 are usually sought as uniform-vorticity flows because of the spatial dependency of the Poincaré term Ω1 c × r. For instance, such solutions for P U 1 0 in ellipsoids have been successfully obtained for latitudinal libration [START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF], precession (e.g. [START_REF] Noir | Precession-driven flows in non-axisymmetric ellipsoids[END_REF], and if we consider longitudinal librations in the mantle frame of reference, a natural solution is

P U 1 0 = -cos(ωt) z R × r.
For the sake of our asymptotic analysis, we consider below a generic uniform-vorticity flow P U 1 0 (see equations B 2-B 3 in appendix B.2), which encompasses all the various cases. Naturally, P U 1 0 = 0 in absence of non-homogeneous forcing terms, as this is for instance the case in the mean rotating frame for longitudinal librations in the spheroid or latitudinal librations in the sphere.

Considering now H U 1 0 , the governing equations are then (3.3a,b) which has to be integrated together with the BC

(∂ t + Ω 0 ∂ φ ) H U 1 0 + 2 Ω 0 s z R × H U 1 0 = -∇P 1 0 , ∇ • ( H U 1 0 ) = 0,
u 1 0 + U 1 0 = V 1 Σ on Σ.
We assume in the mean flow computation below that the spin-up time scale E -1/2 of the fluid is much longer than the characteristic time scale of the mechanical forcing (in the fluid rotating frame), which implies H U 1 0 → 0 (as in Busse 2010; Sauret & Le Dizès 2013). In appendix B.3, we investigate the validity of this assumption, that is how this limit is approached when ω/E 1/2 is increased (for the particular case of longitudinal librations). Note also that assuming H U 1 0 = 0, as in the following, is not valid when bulk flows are generated by the forcing at this order. Considering for instance longitudinal librations (as in Aldridge & Toomre 1969), we detail in appendix B.4 how the excitation of an inertial mode flow H U 1 0 = 0 can indeed modify the mean zonal flow. Then, since U 1 0 = P U 1 0 is known, we can solve equations (2.15) to obtain u 1 0 . The computations of the first-order boundary-layer flow are detailed in appendix B.5, but here we only outline the essential steps. The pressure term in equation (2.15a) is usually removed by multiplying the equation by n × (. . . ) and -i n × ( n × . . . ), which gives L( n × u 1 0 + iu 1 0 ) = 0, with the imaginary number i. While the no-penetration condition u 1 0 • n = 0 gives directly the first component of u 1 0 as u 1 0 • q 1 = 0, the two other components Y = (u 1 0 • q 2 , u 1 0 • φ) can then be obtained by integrating this equation. Given the spatio-temporal periodicity of the perturbation, Y is sought as the linear combination

Y = k Y k , where the individual terms Y k are in the form exp[i(m k φ + ω k t)] to
encompass the various sign possibilities. Equation (2.15a) then reduces to

∂ 2 ζζ Y k = M Y k , (3.4)
where the anti-symmetric matrix M reads (noting

γ k = (ω k + m k Ω 0 )/2) M = 2iγ k -2γ 1 2γ 1 2iγ k (3.5) γ 1 = Ω 0 s z R • q 1 = (Ω 0 + Ω c ) • n = h-1 Ω 0 s T (q1) cos q 2 . (3.6)
where the spherical geometry is recovered for q 1 → ∞, with T (q1) / h → 1.

As detailed in appendix B.5, the linear system (3.4) can be solved together with the no-slip BC to obtain u 1 0 . The resulting mathematical expression shows that the boundarylayer thickness is singular when

γ 1 ± γ ± = 0, that is z R • q 1 = ± γ ± Ω 0 s = ± mΩ 0 ± ω 2 Ω 0 s , (3.7)
with γ ± = (mΩ 0 ± ω)/2. The presence of these singularities shows that boundary-layer theory is not valid at this order of approximation, and their description requires the introduction of new scalings near these so-called critical latitudes (e.g. Kida 2020). The calculations performed in this work are thus strictly valid when |ω ± m| > 2|Ω 0 s |, to prevent the generation of internal shear layers (and the excitation of inertial waves or modes, e.g. [START_REF] Aldridge | Axisymmetric inertial oscillations of a fluid in a rotating spherical container[END_REF]Sauret et al. 2013).

Weakly nonlinear analysis: second-order bulk flows

Using the decomposition

U 2 0 = P U 2 0 + H U 2 0 , the average of equations (2.18) gives 2 Ω 0 s z R × H U 2 0 = -∇ H P 2 0 , ∇ • H U 2 0 = 0, (3.8a,b)
together with BC (2.19b). Note that H U 2 0 is related to viscous effects, and thus, contrary to P U 2 0 , it vanishes when E = 0 (but is non-zero for E 1). Equation (3.8a) admits a solution of the form (e.g. Busse 1968b)

H U 2 0 = sf (s) φ, (3.9)
where the rotation rate f (s) of the mean zonal flow has to be determined. Considering now P U 2 0 , the inhomogeneous forcing term in (2.18) is linear in the Cartesian coordinates [x, y, z], such that we can seek P U 2 0 as a uniform-vorticity flow. For all the forcings considered in this work, the time average of P U 2 0 can be written as sg φ, such that the mean zonal flow reduces to

U 2 0 = s [f (s) + g] φ, (3.10)
where g is a constant, found to be g = 0 in all cases studied here. At the order 2 E 1/2 , the mean zonal component of the bulk flow U 2 1 is governed by

( z R • ∇) U 2 1 = 0, (3.11)
which is actually the Taylor-Proudman theorem [START_REF] Greenspan | The theory of rotating fluids[END_REF]). It implies that the flux ejected out of the boundary layer through the interior (which is symmetric with respect to the axis and anti-symmetric with respect to the equatorial plane) vanishes at every distance from the axis, such that n • u 2 1 | ζ=0 = 0 at the order 2 E 1/2 . Moreover, the continuity equation at order 2 E 0 reads

∂ ζ (u 2 1 • n) = n • ∇ × ( n × u 2 0 ). (3.12)
Finally, integrating equation (3.12) between ζ = 0 and ζ → +∞ yields (e.g. Busse 2010)

n • u 2 1 | ζ=0 = -n • ∇ × +∞ 0 n × u 2 0 dζ = 0, (3.13)
which is used to determine the unknown function f (s).

Weakly nonlinear analysis: second-order boundary-layer flows

To obtain u 2 0 , we separate equation (2.19a) in three distinct problems by considering three distinct velocity fields u A , u B and u C such that u 2 0 = u A + u B + u C . In the first problem, we seek a velocity field u A satisfying the homogeneous equations and the inhomogeneous BC, that is

Hu A + n ∂ ζ Φ 2 1 = 0, u A + U 2 0 = 0 on Σ, (3.14a,b)
where we have defined the linear operator Hu 2 0 = Lu 2 0 , using the operator L defined below equation (2.15). Then, we seek the velocity fields u B and u C that satisfy the homogeneous BC u B = u C = 0 at ζ = 0 and the inhomogenous equations given by

Hu B + n ∂ ζ Φ 2 1 = (u 1 0 • ∇)u 1 0 + 2Ω 1 c × u 1 0 + (U 1 0 • ∇)u 1 0 + (u 1 0 • ∇)U 1 0 , (3.15) Hu C + n ∂ ζ Φ 2 1 = -(u 1 1 + U 1 1 ) • n ∂ ζ u 1 0 . (3.16)
For tidal forcing, [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF] claimed erroneously that the term (u

1 1 + U 1 1 ) • n ∂ ζ u 1 0 vanishes in equation (3.16
), such that the contribution of u C could be discarded. This would be correct if the normal velocity were zero at all orders, but this term only vanishes at the boundary and not everywhere in the boundary layer. We will instead demonstrate that a non-zero u C is required to balance the singularity of u B on the rotation axis. The equations governing u A are formally similar to boundary-layer equations. Similarly, we obtain in the spheroidal coordinates (noting

λ = [1 + i sgn(γ 1 )] |γ 1 |) u A = - s 2 f (s)   0 i (e -λ ζ -e -λ * ζ ) e -λ ζ + e -λ * ζ   .
(3.17)

The calculation of u B is more laborious. After some algebra, equation (3.15) reduces to the following scalar equation

(∂ 2 ζζ -λ 2 )F 0 = E, (3.18)
where

E = (u 1 0 • ∇) u 1 0 + 2Ω 1 c × u 1 0 + (U 1 0 • ∇) u 1 0 + (u 1 0 • ∇) U 1 0 • ( q 2 + i φ), (3.19) with F 0 = u B • ( q 2 + i φ) and λ 2 = 2 i γ 1 .
To calculate E, we first consider each term separately, that is (u

1 0 • ∇) u 1 0 , 2Ω 1 c × u 1 0 , (U 1 0 • ∇) u 1 0 , and (u 1 0 • ∇) U 1 0 .
Then, we decompose each term as a constant term, which contributes to the mean zonal average, and terms proportional to exp(±2iωt), exp(±2imφ) and exp(±2i(mφ ± ωt)) that only contribute to the average if m = 0 or ω = 0.

Considering each term of E separately, the problem is made simpler by making the ζ dependency explicit, that is by rewriting the equation as (3.20) where the complex coefficients (κ k , ϑ k ) and ς k (which is a linear combination of λ ± , κ ± and their complex conjugates, see equation B 29) are independent of ζ. We can then integrate equation (3.20) by considering each term of the sum separately, which gives

(∂ 2 ζζ -λ 2 )F 0 = k (κ k + ϑ k ζ) e -ς k ζ ,
F 0 = k e -ς k ζ [(ς 2 k -λ 2 ) (κ k + ϑ k ζ) + 2ϑ k ς k ] -e -λ ζ [(ς 2 k -λ 2 )κ k + 2ϑ k ς k ] (ς 2 k -λ 2 ) 2 , (3.21)
and then u B can be obtained using u B • q 2 = e (F 0 ) and

u B • φ = m (F 0 ).
Similarly, the calculation of u C can be reduced to the integration of

(∂ 2 ζζ -λ 2 )H 0 = F (3.22) with H 0 = u C • ( q 2 + i φ)
and where the right-hand side is given by

F = -(u 1 1 + U 1 1 ) • n ∂ ζ [u 1 0 • ( q 2 + i φ)] . (3.23)
The calculation of F requires the expression of (u

1 1 + U 1 1 ) • n. The Ekman pumping u 1 1 • n is obtained using the continuity equation -∂ ζ (u 1 1 • n) + n • ∇ × ( n × u 1 0 ) = 0. (3.24)
Using the expression of u 1 0 , we then obtain u 1 1 • n by integration. Together with the Ekman pumping u 1 1 , an Ekman (bulk) circulation U 1 1 is generated via the no-penetration of the fluid at the boundary, such that (u

1 1 + U 1 1 ) • n = 0 at ζ = 0. We thus obtain U 1 1 • n = -u 1 1 • n| ζ=0 .
Using a similar procedure for u B , we can now calculate the analytical expression of F by considering the terms contributing to the average, in particular when ω = 0 or m = 0. We obtain similarly H 0 , and thus u C , by summing all the solutions.

Mean axisymmetric zonal flow

Having explicitly obtained u B and u C in section 3.3, one now use equation (3.13),

n • ∇ × +∞ 0 n × (u A + u B + u C ) dζ = 0, (3.25)
to obtain the unknown rotation rate f (s) present in u A . Therefore, using the expression of u A given by equation (3.17), we obtain

∂ ∂q 2 - sgn(γ 1 ) a T (q1) sin 2 q 2 2 |γ 1 | f (s) + +∞ 0 sin q 2 [ e (F 0 ) + e (H 0 )] dζ = 0, (3.26) which gives f (s) = 2 |γ 1 | sgn(γ 1 ) a T (q1) sin q 2 e ∞ 0 [F 0 + H 0 ] dζ (3.27)
where the integration constant has to be taken equal to 0 to avoid the divergence of the zonal flow when s → 0. From a practical point of view, one can notice that the primitive function G of F 0 + H 0 tends to 0 for ζ → ∞ in order to ensure a zero flux at ζ = ∞, such that equation (3.27) simplifies into

f (s) = - 2 |γ 1 | sgn(γ 1 ) s e G (ζ=0) , (3.28)
which gives the axisymmetric mean zonal flow through equation (3.10).

0 0.2 0.4 0.6 0.8 1 s -0.4 -0.2 0 0.2 f ω = 0.1 ω = 3/2 ω = 2 (a) 0 0.2 0.4 0.6 0.8 1 s -0.3 -0.15 0 0.15 f ω = 2 ω = 3 ω = 4 (b)

Results

In the DNS, we find that the geostrophic flows are produced in an O(E -1/2 ) interval of time, where E -1/2 is the spin-up time scale [START_REF] Greenspan | The theory of rotating fluids[END_REF]). Therefore, for the sake of numerical convergence, we have first integrated the nonlinear equations during a few spin-up times, and then time-averaged the flows over a few tens of forcing periods 2π/ω to extract the mean geostrophic component from the three-dimensional velocity field.

Longitudinal librations

We consider weak longitudinal librations in the mean rotating frame with Ω c = z R , U = 0, and V Σ = cos(ωt) z R × r. In this reference frame, the zonal flow has first been studied theoretically in the sphere by [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF] in the limit of vanishing libration frequency ω → 0. Using a mathematical description in terms of a stream function, [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF] extended the spherical theory to spherical shells and with an arbitrary libration frequency (but still neglecting the shear layers). To avoid the presence of critical shear layers, we only present here results for libration frequencies |ω| 2, The shear layers indeed excite inertial waves occurring when |ω| < 2, as obtained from (3.7), and modify the zonal flow [START_REF] Sauret | Fluid flows in a librating cylinder[END_REF], in the cylindrical geometry). We obtain an excellent quantitative agreement with the results of [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF] in a full sphere, as shown in figure 4(a), which validates our analytical theory. We naturally obtain the same results when calculating the theory in the mantle frame, in which the three last terms of equation 2.19a are now non zero (their contributions to f balance each other).

Since our theoretical approach closely follows the one of [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF], we aim at comparing our results with Busse's theoretical zonal flow, which surprisingly differs from the one obtained by [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF] for the full sphere librating at ω → 0. Indeed, [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF] [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF], the two expressions differ significantly when s > 0.7. The latter authors attributed this difference to their assumption ω 1, supposedly different from the assumption ω of [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF]. Actually, our asymptotic theory follows closely [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF], but our results are in exact agreement with the zonal flow profile of [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF] as shown for ω = 0.1 in figure 4(a). We have thus investigated the origin of this intriguing discrepancy by replicating step by step the calculations of [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF]. We found that his equations are correct, contrary to his integration of the weakly nonlinear inhomogeneous equations, i.e. equations (A5)-(A7) are erroneous. Performing the calculations of [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF] with a computer algebra system gives indeed (4.1b) in the relevant limit ω → 0. To further assess the validity of the asymptotic theory, we have also performed DNS with xshells in the same frame of reference, rotating at Ω c = z R . Considering extremely small viscosity and forcing amplitude (i.e. E = 10 -7 , = 10 -4 ), the numerical flows agree very well with the theoretical predictions (figure 4b). This clearly confirms the agreement already obtained by [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF] at more moderate parameters. For ω = 2, note the small discrepancy at s = 0, due to the presence of the critical latitude.

Finally, we investigate how the zonal flows are modified in spheroids. Rapidly rotating planetary bodies are indeed deformed into ellipsoids due to centrifugal deformations, and several laboratory experiments have been designed such as the ZoRo experiment [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF][START_REF] Vidal | Compressible fluid modes in rigid ellipsoids: towards modal acoustic velocimetry[END_REF]) with r pol /r eq = 0.95. We perform DNS with Nek5000 in the mantle frame of reference (where the boundary velocity is zero), and present the mean zonal flows obtained from spheroidal DNS at E = 10 -5 for various values of the ratio r pol /r eq in figure 5(a). We also compare the results to the theoretical profiles that have been obtained in spheroidal coordinates. Overall, we find a good quantitative agreement, even if the DNS have not been performed in the regime E 1. The numerical results convincingly validate our asymptotic theory of libration-driven zonal flows in spheroids. It is worth noting that significant departure from the spherical profile is found, even for moderate spheroidal deformations r pol /r eq 1 as often considered experimentally (e.g. 0.7 in [START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF]. The theory also allows us to explore more extreme spheroidal configurations that cannot be simulated numerically, as illustrated in figure 5(b). Two points are worthy of comments. We find that the mean zonal flow reaches a constant value at s = 0. The latter value actually corresponds to the constant profile obtained in the cylinder [START_REF] Wang | Cylindrical tank of fluid oscillating about a state of steady rotation[END_REF][START_REF] Sauret | Fluid flows in a librating cylinder[END_REF], which gives a lower bound for f . In the interior 0 < s < 1, f tends again to the cylindrical value in the disc limit, that is r pol /r eq → 0, whereas it vanishes in the infinite cylinder limit r pol /r eq → ∞. Moreover, our results illustrate that the cylindrical geometry cannot be faithfully used as a reduced model of the spheroid. Therefore, results obtained in a cylindrical geometry should be interpreted with caution for planetary applications.

Latitudinal librations

We now consider the flows driven by latitudinal librations, which have only received scant attention so far [START_REF] Chan | Simulations of fluid motion in spheroidal planetary cores driven by latitudinal libration[END_REF][START_REF] Zhang | Asymptotic theory of resonant flow in a spheroidal cavity driven by latitudinal libration[END_REF][START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF]. In particular, the mean zonal flows have only been computed numerically at moderate values of E [START_REF] Chan | Simulations of fluid motion in spheroidal planetary cores driven by latitudinal libration[END_REF], and never compared to theoretical predictions. Temporal and spatial perturbations must be indeed considered simultaneously, respectively at the frequency ω and at the azimuthal wavenumber m = 1. This approach contrasts with previous theoretical studies of zonal flows, where only one kind of perturbations was considered, and fully justifies the generic theoretical framework presented in §3.

We consider for simplicity the spherical geometry, and we perform our analytical and numerical calculations in the mean rotating frame with

Ω c = z R = z I , V Σ = cos(ωt)
x R × r, and U = 0. Using our asymptotic approach, we uncover the theoretical zonal flow associated with this forcing in the relevant limit of vanishing viscosity. We compare the associated theoretical profiles with DNS in figure 6(a). We obtain an excellent agreement for the three different libration frequencies and for very small perturbation and viscosity (E = 10 -7 , = 10 -4 ) in the regime ω 2 (i.e. without critical latitudes). Note that the rotation rate is always regular at s = 0 in the DNS (as mathematically expected from Lewis & Bellan 1990, see also appendix A), but our theoretical profile diverges at s = 0 for ω = 2. Indeed, the mathematical singularity associated with the critical latitude is located on the rotation axis for ω = 2 (see equation 3.7). This mathematical singularity is smoothed out by viscosity in the DNS but, to regularise our asymptotic theory and obtain a regular rotation rate profile everywhere in space, additional viscous effects (e.g. [START_REF] Kida | Steady flow in a rapidly rotating sphere with weak precession[END_REF]) should be taken into account at the critical latitudes. Finally, we can explore with the theory how the rotation rate evolves with the libration frequency. Similarly to figure 4 for longitudinal librations, we illustrate in figure 6(b) the theoretical zonal flows for various libration frequencies in the particular regime ω < 2 (where the theory may not be valid, which will be further discussed in §5). Even if higher-order viscous effects are expected to smooth out the singularity, one can already notice that, at this order, the width of the divergence seems to increase as ω → 2. The influence of the critical latitudes on the zonal flows is further discussed below.

Precession-driven zonal flows

In his seminal work, Busse (1968b) considered a precessing sphere and found that the first-order bulk flow is a solid-body rotation ω f × r, tilted from the boundary rotation vector Ω 0 z R . Having shown that 2 = (Ω 0 z Rω f ) 2 = Ω 2 0ω 2 f , he showed that the component of Ω c normal to ω f is of the order E 1/2 , which can thus be neglected at the order of the mean zonal flow calculation. To calculate the mean zonal flow, he then neglected Ω c • ω f with respect to ω 2 f = |ω f | 2 for simplicity, yielding finally Ω c = 0 (his mean zonal flow results are thus obtained in the inertial frame of reference). However, only a crude agreement has been found experimentally [START_REF] Malkus | Precession of the Earth as the cause of geomagnetism: Experiments lend support to the proposal that precessional torques drive the Earth's dynamo[END_REF]) and numerically (Noir et al. 2001b) with his theoretical zonal flow (see figure 1). To carefully compare theory and numerics, we have performed DNS in spheres to explore smaller values of E than in ellipsoids, and also to consider very small precession angles. Moreover, we have taken Ω c • ω f into account in our theoretical calculations (contrary to Busse 1968b).

While DNS are performed in the precessing frame described in section 2.2, it is more convenient for theoretical calculations to consider the reference frame where the z-axis is along ω f . Then, the boundary velocity can be written as

V Σ = Ω 0 z R × r + x R × r, with Ω 0 = (1+P o ) -1 .
Since only the component of Ω c parallel to ω f is important at this order, equation (2.7) can be simplified to obtain Ω c Ω 0 P o cos(α) z R . without loss of generality. We show in figure 7(a) that the results of Busse (1968b) are recovered in the limit P o = 0 (where the frame of reference reduces to the inertial frame). We also show the mean flow for small P o = 0 (i.e. when Ω c • ω f is taken into account). Considering P o = 0 could be relevant for DNS performed at moderate values of P o and α, but the mean flow is expected to be only marginally modified in the planetary regime P o 1. We have next performed DNS at low values of E to validate the asymptotic theory of Busse (1968b), since previous studies have failed to recover the theoretical mean zonal flows (see figure 1 above). The simulations have been computed in the precession frame, where the rotation vector ω f of the basic flow is accurately given by [START_REF] Cébron | Precessing spherical shells: flows, dissipation, dynamo and the lunar core[END_REF])

ω f • x R = 1 1 + P o [λ i + γ cos(α)] γ sin(α) γ[γ + 2λ i cos(α)] + |λ| 2 , (4.2a
)

ω f • y R = - 1 1 + P o γλ r sin(α) γ[γ + 2λ i cos(α)] + |λ| 2 , (4.2b) ω f • z R = 1 1 + P o γ[γ cos 2 (α) + 2λ i cos(α)] + |λ| 2 γ[γ + 2λ i cos(α)] + |λ| 2 , ( 4.2c)
with γ = P o /E 1/2 and where λ = λ r + iλ i is the spin-over damping factor given by [START_REF] Hollerbach | Oscillatory internal shear layers in rotating and precessing flows[END_REF]Noir et al. 2001b)

λ r ≈ -2.62 -1.36 E 0.27 , λ i ≈ 0.258 + 1.25 E 0.21 . (4.3a,b)
To compare the simulations with the theory, which assumes that ω f is along z I at leading order (Busse 1968b), we post-process the data as follows. We introduce a new reference frame, called the fluid frame, where the new z-axis is along the axis of rotation given by equation (4.2). Then, we rotate the velocity field into that frame and, to isolate the mean zonal component of order 2 from the leading-order steady uniform-vorticity flow given in equation (4.2), we compute the rotation rate of the mean zonal flow as from the axis (when 0.3 s 0.7). We now investigate how these axial values f (0) vary with the control parameters, and the oscillations of f are discussed further in section 5.

f (s) = V φ (z = 0) -ω f s 2 s , ( 4
Our DNS show that the axial value f (0) depends linearly on γ, as illustrated in figure 8(a). However, note that the azimuthal velocity sf (s) always vanishes at s = 0. Thus, the value f (0) only governs the weak slope of the velocity at s = 0, which can vary in the DNS. In the following, for every E, we have varied P o to obtain the value of γ which cancels out f (0). The corresponding DNS are shown in figure 8(b). Similarly to figure 7(b), the numerical profiles are in good agreement with the theory for E 10 -6 for 0.3 s 0.7. Near the critical latitude located at s = √ 3/2 0.866, the mathematical singularity is smoothed out by viscous effects. However, the lower the viscosity, the better the agreement with the theory on both sides of the singularity, with the mean zonal flow converging towards the theory when E is reduced (see figure 8b).

4.4.

Revisiting the tidal-like forcing of [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF] We finally consider the tidal-like forcing, first considered analytically by [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF], assuming ω = 0 and m = 2 in boundary flow (2.4) together with Ω c = 0. Suess (1971) investigated theoretically and experimentally the mean zonal flow generated by this forcing. [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF] predicted theoretically the generation of a strong retrograde vortex along the rotation axis, and his prediction was surprisingly in broad agreement with his experimental results as reproduced in figure 9(a). However, the experimental profile cannot be singular at s = 0 (see appendix A), and the reported broad agreement with theory has thus to be erroneous. Actually, a thorough analysis reveals that even his theory is incorrect because (i) he erroneously discarded the contribution from u C and, (ii) because he made some errors in his calculation (e.g. his equation 42 is incorrect). A possible singularity at s = 0 may actually be expected from the divergence of the firstorder boundary-layer solution u 1 0 at s = 0, and figure 9(b) shows that the contributions of u B and u C are indeed singular at s = 0. Nevertheless, figure 9(b) clearly shows that the total zonal flow, given by the sum of the contributions, is smooth because they exactly balance each other on the axis at s = 0. This confirms the crucial role of the flow u C , which was discarded by [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF]. (1971). Dashed blue curve: erroneous theory of [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF]. Inset shows the azimuthal velocity sf . (b) Various contributions to f (solid red curve), made of the sum of two contributions: the one due to uB (dotted curve), which was mistakenly considered to be f by [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF], and the one due to uC (dashed curve).

We have drawn the correct theoretical solution in figure 9(a), clearly showing that it does not agree with the experimental measurements (especially at s = 0). One may wonder whether the observed strong retrograde flow near the axis of rotation is reminiscent of the non-zero values of f (0) reported above for precession and libration at ω 2. However, since the axial value sf (0) of the azimuthal velocity does not seem to vanish (contrary to the corresponding value for precessing flows), this flow may have a different origin. Indeed, the experimental results might instead exhibit the saturation of an elliptical instability in the bulk of the fluid, which is not taken into account in the theory (as it was not known at this time, see the review in [START_REF] Kerswell | Elliptical instability[END_REF]. The elliptical instability can be excited when the streamline ellipticity β is large compared to the viscous dissipation (with β = 2 in Suess 1971). More quantitatively, the growth rate σ of the elliptical instability can be estimated as σ = 9β/16 -KE 1/2 , with the typical values 1 K 10 related to viscous dissipation in the Ekman boundary layer (e.g. [START_REF] Cébron | Elliptical instability in terrestrial planets and moons[END_REF]. Since the experimental results of [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF] shown in figure 9(a) have been obtained for E 1/2 = 5.3 × 10 -3 and β = 2.5 × 10 -2 , the expected growth rate is -0.04 σ 10 -2 . Thus, an elliptical instability might have been excited in the experiment, which would certainly modify the observed zonal flow (see e.g. figure 5 in Grannan et al. 2017, obtained for Ω c = -1).

Discussion

Mean zonal flows when ω 2

We have so far successfully validated the theoretical mean zonal flows driven by librations when ω > 2, that is in the absence of inertial waves and critical latitudes. A successful validation was, however, less straightforward to obtain for the precession forcing. Precession-driven flows are indeed strongly affected by the presence of inertial waves and conical shear layers spawned from the critical latitudes (Noir et al. 2001a). However, we have still found an overall good quantitative agreement with the analytical model, even if the latter has been obtained by neglecting these two additional effects. . The curve at E = 10 -7 is indistinguishable from the one obtained with = 10 -5 (not shown). In (b), DNS have been performed at E = 10 -7 , and we have checked that the numerical profile f at ω = 1 is unchanged for = 10 -3 and = 10 -6 (not shown). Insets show that f (0) remains finite in DNS.

By analogy, more complicated flow structures are also expected for libration-driven flows when ω 2, as previously reported in cylinders [START_REF] Sauret | Fluid flows in a librating cylinder[END_REF]) and spheres [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF]. Critical latitudes indeed also exist for libration-driven flows when ω 2 but, based on our findings for precession, one may still expect a rather good agreement between the analytical libration-driven mean zonal flows and DNS.

We illustrate in figure 10 the mean zonal flows in rotating spheres subject to latitudinal librations in panel (a) and longitudinal librations in panel (b). Outside the region affected by the critical latitudes, the rotation rate f exhibits oscillations in the bulk, which are very similar to the ones reported for precession in §4.3, that are superimposed onto the theoretical profile. In particular, the fluid rotation rate at the axis f (0) can be quite different than the theoretical profile. Because the finite difference method loses some accuracy at r = 0 (see appendix A), we have carefully checked numerical convergence at E = 10 -7 by refining the grid (especially at the centre) and also by increasing the maximum spherical harmonic degree up to l max = 339 (and up to l max = 425 at E = 10 -8 ). The variations are within the thickness of the curve.

In the case of the latitudinal librations shown in figure 10(a), the same stationary oscillations are also present when only |m| 1 are considered in the DNS. Careful inspection of the instantaneous flow reveals that forced inertial modes (equatorially antisymmetric m = ±1 for latitudinal librations, and equatorially symmetric m = 0 for longitudinal librations) are present, which have a number of radial zeros. We conjecture that these modes have almost the same frequency as the excitation (here ω 0.1) and produce the multiple jets by nonlinear interaction, which carry the signature of inertial modes. In addition, we have computed the same case at a lower viscosity (E = 10 -8 ) also shown in figure 10a. It highlights that once viscosity is low enough for the oscillations to appear, their amplitude is nearly independent of the Ekman number. Having changed given by the theory (black dotted curve for the sphere, black dashed curve for r pol /req = 0.8) and DNS (coloured solid curves) at E = 10 -5 and = 10 -4 . (b) Comparison of the theoretical rotation rate between the sphere, an oblate (r pol /req = 0.2) or prolate (r pol /req = 4) spheroid.

in figure 10, we obtain that the amplitude of such oscillations does not depend on when 1. The zonal jet velocities scale thus as 2 E 0 , exactly as the theoretical mean zonal flow. Such scaling law is consistent with a zonal flow associated with a librationexcited mode, which has an amplitude of the order E 0 (see appendix B.4). The precise mechanism is beyond the scope of this paper, but it is worth noting that Lin & Noir (2020) also observed an imprint of the excited inertial mode on the mean zonal flow in shells. Despite the presence of these jets, our asymptotic theory fairly reproduces the mean zonal flows found in the DNS (see appendix B.4 for further details).

The same conclusion is drawn in spheroids, as reported in figure 11(a), even if the DNS in spheroidal geometries have been performed for moderate values of E 10 -5 (compared to DNS in spheres with E 10 -7 ). Another striking point in the figure is that the location of the critical latitude varies with the ellipsoidal deformation. Indeed, the unit vector e 1 in the left-hand side of equation (3.7) depends on the spheroidal geometry, such that the spatial position of the critical latitude is modified. This phenomenon is further illustrated with more deformed spheroidal geometries in figure 11(b). This effect has direct consequences for the numerical profiles obtained. Only a small part of the volume is affected by the shear layer when the critical latitude is close to the boundary (see in figure 11a for r pol /r eq = 0.8), such that the oscillations are rather localised around the position of the critical latitude. However, the oscillations can spread in the volume when the critical latitude is far from the boundary (see r pol /r eq = 1 in figure 11a).

Geostrophic shear spawned from the critical latitudes when E 1

The analytical solutions show that the mean zonal flows are singular at the cylindrical radius s c when ω 2 (see e.g. figure 11b), Nevertheless, this singularity is regularised by viscosity such that the mean zonal flow takes the form of a shear layer near s c , as first noticed by Busse (1968b). Moreover, it is known that the amplitude of the geostrophic shear increases as E is reduced, by contrast with the typical horizontal length scale of the shear (e.g. Noir et al. 2001b, for precession). However, the corresponding scaling laws with the Ekman number have been disputed, and no conclusive answer has been obtained yet. We aim at revisiting here that problem numerically with DNS in spheres, to hopefully capture the correct asymptotic behaviour in the relevant regime E 1. [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF] recently explored the geostrophic shear generated by longitudinal librations with ω = 1 in spherical shells, but with a small inner core (s 0.1). We reproduce their results in figure 12(a), and find that their mean zonal flows are in very good agreement with the theoretical solution at ω = 1 in the full sphere. Considering now the width δs and the peak-to-peak amplitude δu g of the geostrophic shear (as defined in figure 12a), [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF] proposed the scaling laws (5.1a,b) using DNS of longitudinal librations at ω = 1 and order-of-magnitude arguments. We confirm numerically scaling laws (5.1) for various libration-driven flows at different ω, but only in a certain parameter range (see figure 12b). One can indeed anticipate a possible change of regime when the geostrophic shear, of typical thickness E 1/5 [START_REF] Kerswell | On the internal shear layers spawned by the critical regions in oscillatory ekman boundary layers[END_REF]Noir et al. 2001b) and centred on the cylindrical radius s c (see inset in figure 12b), interacts with the equatorial boundary layer at s 1. Because the typical Ekman layer thickness E 1/2 is negligible with respect to E 1/5 when E 1, we expect a different behaviour for a certain value ω = ω c given by

δu g / 2 ∝ E -1/10 , δs ∝ E 1/5 ,
s c + O(E 1/5 ) 1, (5.2)
with the cylindrical radius of the critical latitude s c = sin cos -1 (ω c /2) ≈ 1ω 2 c /8 for small values of ω c in the sphere. Equation (5.2) gives ω c = O(E 1/10 ) 1 in the regime E 1 (see figure 12b). When ω E 1/10 , we thus expect the Ekman layer scaling laws (5.3a,b) In the opposite regime ω E 1/10 , the relevant scaling laws should be (5.1) as proposed by [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF]. In order to validate these theoretical considerations, we have

δu g / 2 ∝ E -1/2 , δs ∝ E 1/2 .
10 -7 10 -6 10 -5 E 10 -4 10 -3 10 -2 10 -1 10 0 δs ∝ E 1/5 ∝ E 1/2 ω = 0.1 ω = 1 ω = 0.1 ω = 1.5 10 -7 10 -6 10 -5 E 10 -2 10 -1 10 0 10 1 δu g / 2 ∝ E -1/10 ∝ E -1/2 ∝ E -0.4 ω = 0.1 ω = 1 ω = 0.1 ω = 1.5 (a) (b)
Figure 13. Geostrophic shear associated with the critical latitudes of longitudinal librations (empty points) and latitudinal librations (filled points) in a The symbols ω = 1 performed for longitudinal librations have been extracted from figure 10 in [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF]. The black cross indicates our DNS in a full sphere at ω = 1 and = 10 -4 , which agrees with [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF]. Distance δs between peaks in (a), and peak-to-peak amplitude δug/ 2 in (b), as a function of E.

performed DNS in both regimes at different libration frequencies, especially near the transition between these two configurations (figure 12b). The numerical results, obtained for various libration forcings and frequencies, are summarised in figure 13. The various scaling laws are numerically recovered, as well as the change of regimes. Note that the typical frequency of most planetary bodies subject to longitudinal librations is ω O(1) (e.g. [START_REF] Noir | An experimental and numerical study of librationally driven flow in planetary cores and subsurface oceans[END_REF][START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF], such that scaling laws (5.1) are expected to be relevant for most planetary applications. Finally, the aforementioned scaling laws differ from the ones that have been proposed for precession-driven flows (Noir et al. 2001b), that is

δu g / 2 ∝ E -3/10 , δs ∝ E 1/5 .
(5.4a,b)

We have checked that we also recover scaling laws (5.4) for our DNS of precession-driven flows. Note that the scaling law (5.4a) has also been experimentally observed in a rotating sphere subject to a tidal deformation [START_REF] Morize | Experimental determination of zonal winds driven by tides[END_REF]. As outlined in [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF], the origin of the different scaling laws between libration and precession remains puzzling.

Influence of a solid inner core

Planetary fluid layers are often bounded by two solid layers (e.g. the Earth's liquid core is surrounded by the upper mantle and a solid inner core). One can thus wonder how our results, obtained in coreless geometries, could be modified by the presence of an inner boundary. We focus here on libration-driven zonal flows, which have already received attention (e.g. [START_REF] Calkins | Axisymmetric simulations of libration-driven fluid dynamics in a spherical shell geometry[END_REF][START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF][START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF].

We first consider the case where the critical latitudes and inertial waves are absent. This regime has been theoretically investigated for longitudinal librations in [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF], showing that the mean zonal flows in a spherical shell can be entirely deduced from the solutions in the full sphere (if we exclude the Stewartson layers associated with the presence of the inner boundary, see [START_REF] Stewartson | On almost rigid rotations[END_REF]. Here, we consider a (possibly nonhomoeoidal) spheroidal shell, where the inner boundary is spheroidal, with respectively Figure 14. Normalised mean zonal velocity V φ / 2 at various dimensionless heights z = √ 1s 2 in homoeoidal shell geometries (with aspect ratios η pol = ηeq = 0.5) for libration forcings. The DNS profiles are computed at z = z0 for s ηeq, and at z0 = 0 for s > ηeq, giving a single profile for each DNS. Grey area shows the tangent cylinder s ηeq. Forcings with ωin = ω and in = on the no-slip inner boundary, which is subject to the same forcing as the outer boundary (except for the stress-free case, labelled SF). Black dashed curves indicate the theoretical profiles in the shell geometry. (a) Spherical shell with a no-slip (cyan and orange solid curves) or a stress-free (red and pink solid curves) inner boundary. Blued dotted dashed curves illustrate the full-sphere analytical profiles. DNS performed at E = 10 -6 , ω = 3, and = 10 -6 for the two forcings. (b) Homoeoidal (i.e. η pol = ηeq) spheroidal shells subject to longitudinal librations (solid coloured curves). DNS performed at E = 2.5 × 10 -6 , ω = π, and = 5 × 10 -4 . the (dimensional) inner equatorial r in eq = η eq r eq and polar r in pol = η pol r pol axes, where [η eq , η pol ] are the equatorial and polar shell aspect ratios (η eq = η pol in homoeoidal shells). We impose on the inner boundary a harmonic tangential velocity of magnitude in and angular frequency ω in , which may differ from the forcing at the outer boundary (with the amplitude and angular frequency ω as above). Following [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF], we find that the rotation rate of the mean zonal flow is given in dimensionless form by

V φ s =      2 in (1 -s 2 ) 1/4 f icb sp (s/r in eq ) + 2 [1 -(s/r in eq ) 2 ] 1/4 f cmb sp (s) (1 -s 2 ) 1/4 + [1 -(s/r in eq ) 2 ] 1/4
for s < r in eq , 2 f cmb sp for s > r in eq ,

(5.5)

where f icb sp (respectively f cmb sp ) is the rotation rate profile of the mean zonal flow in a coreless geometry when the forcing at the inner (respectively outer) boundary is considered. According to equation (5.5), the presence of an inner boundary at s = r in eq is not expected to modify the mean zonal flow for s > r in eq (i.e. outside the tangent cylinder). Note that [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF] only considered the particular situation ω i = ω for inner and outer boundaries subject to longitudinal librations, our expression (5.5) naturally agrees with their formula (4.24) in this case. DNS in spherical shells are in excellent agreement with formula (5.5) as shown in figure 14(a), even for latitudinal librations not considered in [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF]. Note that the observed discontinuity at s = r in eq is related to the presence of the Stewartson layers due to the velocity mismatch between the zonal bulk flow and the inner boundary (these layers are absent for a stress-free inner boundary, as found in figure 14a). Introducing nested viscous layers would be required to smooth out the singularity at the Stewartson layer [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF]. Another striking point in figure 14(a) is that considering a stress-free inner boundary does not Grey area shows the tangent cylinder s ηeq. (a) Non-homoeoidal shells subject to longitudinal librations (solid coloured curves), with a spherical inner boundary (i.e. r in pol = r in eq ) and a spheroidal outer boundary (i.e. r pol = req). DNS performed at E = 2.5 × 10 -6 , ω = ωin = π, = in = 5 × 10 -4 . (b) Homoeoidal shells with η pol = ηeq = 0.5. Case 1: DNS at E = 10 -6 in a spherical shell subject to latitudinal librations at outer boundary (with ω = 3 and = 10 -6 ), and to longitudinal librations at inner boundary (with ωin = 2, in/ = 2). Case 2: Spheroidal shell with r pol /req = 0.8 subject to longitudinal librations. DNS at E = 5 × 10 -6 with = 5 × 10 -4 and ω = π at outer boundary, and with in/ = ωin/ω = 2 at inner boundary. modify the mean zonal flow when s r in eq (red curve). Indeed, if the flow obeys stress-free conditions on the inner boundary, the corresponding mean flow in the tangent cylinder is only generated by nonlinear interactions within the Ekman layer at the outer boundary in formula (5.5), and so we recover the coreless solution when s r in eq . The agreement with DNS is also very good in homoeoidal shells (i.e. η pol = η eq ), see in figure 14(b).

Figure 15 shows that formula (5.5) is also valid for other configurations. An excellent agreement is found in non-homeoidal spheroidal shells (i.e. η pol = η eq with r eq = r pol ) as shown in figure 15(a) or, as illustrated in figure 15(b), in homoeoidal shells with distinct angular frequencies ω in = ω and magnitudes in = (purple curve), or in the presence of different kind of mechanical forcings at inner and outer boundaries (orange curve, when the inner boundary undergoes longitudinal librations and the outer one latitudinal librations). This confirms that formula (5.5) is valid even for such complicated cases.

We have obtained and validated so far formula (5.5) for libration angular frequencies larger than 2 (see also [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF]. One can thus wonder how this formula compares with DNS when critical latitudes are present. Such a situation has been recently considered for longitudinal librations in [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF], for the particular libration frequency √ 2 which is associated with conical shear layers spawned form the critical latitudes leading to a simple closed trajectory for the forced inertial wave in a spherical shell with η pol = η eq = 0.35 (e.g. [START_REF] Rieutord | Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum[END_REF]. We revisit their results in figure 16, reproducing their figure 13(a) in our panel (a), for which only the inner boundary is subject to librations, and their figure 17(a) in our panel (b) that corresponds to the opposite situation. While it is very challenging to analytically tackle properly this problem, we find that formula (5.5) provides a reasonably good agreement with DNS, by capturing the essential features of the mean zonal flow profile. [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF], where V φ is computed as in formula (2.8), with zmax = 1, but only considering only the m = 0 component. Theory is given by formula (5.5). Grey area indicate the tangent cylinder s 0.35.

Therefore, even in the presence of additional complicated flow structures and waves in shell geometries, nonlinear interactions within the Ekman boundary layers still make a significant contribution to the mean zonal flows. This agrees with previous findings in librating cylinders (see figure 17 in [START_REF] Sauret | Fluid flows in a librating cylinder[END_REF], which showed that analytical theory obtained in the regime ω > 2 provides the general trend for V φ at ω < 2, on which additional mean flow contributions can be superimposed.

Planetary applications

We have shown that our theory fairly predicts the mean zonal flows in rotations ellipsoids and shells, and for various mechanical forcings. The relevance of these mean zonal flows ought now to be addressed for planetary applications. First, the Ekman boundary layers must be laminar for our theory to be valid. Various mechanisms are known to destabilise laminar Ekman boundary layers, such as Taylor-Görtler instability [START_REF] Noir | An experimental and numerical study of librationally driven flow in planetary cores and subsurface oceans[END_REF][START_REF] Calkins | Axisymmetric simulations of libration-driven fluid dynamics in a spherical shell geometry[END_REF] or local shear instabilities (e.g. [START_REF] Lorenzani | Fluid instabilities in precessing spheroidal cavities[END_REF]. The transition between laminar and turbulent Ekman boundary layers would here occur when ∼ K E 1/2 [START_REF] Cébron | Precessing spherical shells: flows, dissipation, dynamo and the lunar core[END_REF], where K is a numerical prefactor. The first boundary-layer instabilities would occur when K = 20 -55 (e.g. [START_REF] Lorenzani | Fluid instabilities in precessing spheroidal cavities[END_REF][START_REF] Noir | An experimental and numerical study of librationally driven flow in planetary cores and subsurface oceans[END_REF][START_REF] Calkins | Axisymmetric simulations of libration-driven fluid dynamics in a spherical shell geometry[END_REF], and fully turbulent boundary layers are expected when K 150 (e.g. [START_REF] Caldwell | Characteristics of Ekman boundary layer instabilities[END_REF][START_REF] Sous | Friction law and turbulent properties in a laboratory Ekman boundary layer[END_REF]. For precession and latitudinal librations, note that the conical shear layers spawned from the critical latitudes (either at inner or outer boundaries) can also be prone to shear instabilities (e.g. the conical shear instability, see in [START_REF] Lin | Shear-driven parametric instability in a precessing sphere[END_REF][START_REF] Horimoto | Conical shear-driven parametric instability of steady flow in precessing spheroids[END_REF].

However, the Ekman boundary layers are expected to become turbulent before the onset of such instabilities when E 1 (see figure 6 in [START_REF] Cébron | Precessing spherical shells: flows, dissipation, dynamo and the lunar core[END_REF]. Next, in addition to laminar boundary layers, our theory also assumes laminar bulk flows. Bulk turbulence may indeed alter the mean zonal flows, as previously reported for strong tidal or libration forcings (e.g. [START_REF] Favier | Generation and maintenance of bulk turbulence by libration-driven elliptical instability[END_REF][START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF]. To characterise the forcings, we introduce the dimensionless Rossby number Ro = U/(Ω s r eq ), where U is the typical amplitude of the nonlinear flows driven by the mechanical forcings (based on control parameters). Laminar bulk flows are known to be destabilised by several instabilities when Ro KE 1/2 (e.g. the elliptical instability), where Ro is also here the typical inviscid growth rate of the instability and K = 1 -10 is a numerical pre-factor due to Ekman pumping [START_REF] Lemasquerier | Libration-driven flows in ellipsoidal shells[END_REF]. In ellipsoids, we have at leading order Ro ∼ β for tides (e.g. [START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF]Vidal & Cébron 2017), and Ro ∼ β for topographic precession [START_REF] Kerswell | The instability of precessing flow[END_REF] or libration forcings [START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF][START_REF] Vidal | Fossil field decay due to nonlinear tides in massive binaries[END_REF], where β is here a typical measure of the boundary (equatorial or polar) ellipticity. Several secondary instability mechanisms could then occur to sustain bulk-driven zonal flows. Although different in nature, these various scenarios are due to nonlinear bulk interactions and apparently all operate on the dimensionless time scale of order (kRo) -2 in the rapidly rotating planetary regime Ro 1 (e.g. [START_REF] Kerswell | Secondary instabilities in rapidly rotating fluids: inertial wave breakdown[END_REF][START_REF] Brunet | Shortcut to geostrophy in wave-driven rotating turbulence: the quartetic instability[END_REF][START_REF] Le Reun | Near-resonant instability of geostrophic modes: beyond Greenspan's theorem[END_REF], where k is a typical wavenumber of the flow. On the contrary, boundary-layer interactions establish geostrophic flows on the spin-up time scale E -1/2 . The bulk mechanisms should thus operate faster than our boundary-driven mechanism when k 2 Ro 2 E 1/2 , giving Ro E 1/4 /k (which is also the threshold onset for the secondary instabilities, see in [START_REF] Kerswell | Secondary instabilities in rapidly rotating fluids: inertial wave breakdown[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 typical wavenumber for the aforementioned bulk mechanisms is poorly constrained from previous studies (as these mechanisms have only been explored without taking the beta effect into account), such that a rigorous scaling law is still unknown.

For a direct comparison with the boundary-driven zonal flows, we assume k = 1 -10 (i.e. to focus on the large-scale components of the geostrophic flows) and illustrate the corresponding stability diagrams for the boundary-driven and bulk-driven mechanisms in figure 17. Typical planetary values are E = 10 -15 -10 -12 and = 10 -7 -10 -3 , depending on the considered forcing. We thus expect laminar Ekman layers in several planetary bodies, as observed in figure 17(a). Figure 17(b) then clearly indicates that the bulk-driven mechanisms are likely irrelevant to explain the occurrence of mean zonal flows in planetary bodies. On the contrary, several planetary bodies may have simultaneously laminar boundary layers and no bulk-driven turbulence, such that nonlinear interactions within the laminar Ekman layers could be important to generate mean zonal flows.

Conclusion

Summary

In this work, we have investigated the mean zonal geostrophic flows in rapidly rotating spheres and spheroids subject to weak mechanical forcings (librations, precession and tides). Geostrophic flows are indeed often encountered in geophysical or astrophysical systems, which are usually attributed to nonlinear interactions occurring at a small scales (e.g. [START_REF] Christensen | Zonal flow driven by strongly supercritical convection in rotating spherical shells[END_REF][START_REF] Aubert | Observations of zonal flow created by potential vorticity mixing in a rotating fluid[END_REF]. However, the external mechanical forcings can generate large-scale geostrophic flows in the bulk by nonlinear viscous effects, as considered here. We have presented a generic asymptotic theory accounting for the various forcings, in the double limit of small Ekman numbers E and small (dimensionless) forcing amplitude , and we have analytically considered simultaneously azimuthal and temporal variations of the forcings. We have also assessed the range of validity of the analytical profiles as a function of the forcing frequency ω, using targeted DNS.

For all the forcings, we have shown that the leading-order mean zonal flows in the bulk scale as 2 , and are independent of the Ekman number when ω is greater than twice the rotation rate in the absence of inertial waves (i.e. ω 2 in dimensionless spin units), as previously found in spherical geometries for precession (Busse 1968b) and longitudinal librations [START_REF] Busse | Mean zonal flows generated by librations of a rotating spherical cavity[END_REF][START_REF] Calkins | Axisymmetric simulations of libration-driven fluid dynamics in a spherical shell geometry[END_REF][START_REF] Sauret | Experimental and numerical study of mean zonal flows generated by librations of a rotating spherical cavity[END_REF][START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF][START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF]. Moreover, we have shown that these flows can be significantly modified in spheroids subject to longitudinal librations. Our asymptotic theory provides thus a reliable point of comparison for forthcoming experimental measurements, for instance in the ZoRo experiment (Zonal jets in Rotating fluids, see in [START_REF] Su | Acoustic spectra of a gas-filled rotating spheroid[END_REF][START_REF] Vidal | Compressible fluid modes in rigid ellipsoids: towards modal acoustic velocimetry[END_REF]) that is currently used to investigate libration-driven zonal flows.

Then, the existence of critical latitudes and inertial waves when ω < 2 is known to lead to more complicated mean zonal flows in terms of amplitude and structure. Indeed, the critical shear layers spawned from the critical latitudes (e.g. [START_REF] Kerswell | On the internal shear layers spawned by the critical regions in oscillatory ekman boundary layers[END_REF] are responsible for zonal geostrophic shears at the singular points of the theoretical profiles. We have numerically confirmed that the geostrophic shear driven by longitudinal and latitudinal librations has a typical width δs ∝ E 1/5 and a characteristic amplitude δu g / 2 ∝ E -1/10 when ω E 1/10 , which contrasts with the scaling law δu g / 2 ∝ E -3/10 for the geostrophic shear driven by precession.

Finally, we have investigated how the mean zonal flows are modified in the presence of a solid inner core. We have focused on libration-driven flows to revisit previous numerical findings at low E in shells [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF]. Interestingly, we have shown that the mean zonal flows in homoeoidal shells can be fairly estimated from the coreless solutions, in agreement with previous analytical works [START_REF] Sauret | Fluid flows in a librating cylinder[END_REF][START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF].

Perspectives

Further work remains to be done to get a more complete description of the generation of mean zonal flows in rapidly rotating bodies. In particular, the competition between bulkdriven and boundary-driven zonal flows should be quantitatively investigated in future studies, to go beyond the qualitative picture discussed above. The forcing amplitude was indeed set here to be small enough, to filter out any fluid instabilities that can grow in the bulk (e.g. [START_REF] Kerswell | Elliptical instability[END_REF][START_REF] Lin | Shear-driven parametric instability in a precessing sphere[END_REF][START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF]Vidal & Cébron 2017;[START_REF] Nobili | Hysteresis and instabilities in a spheroid in precession near the resonance with the tilt-over mode[END_REF]. Only a few experimental or numerical works have hitherto studied mean zonal flows in the presence of bulk turbulence (e.g. [START_REF] Favier | Generation and maintenance of bulk turbulence by libration-driven elliptical instability[END_REF][START_REF] Grannan | Tidally forced turbulence in planetary interiors[END_REF][START_REF] Le Reun | Experimental study of the nonlinear saturation of the elliptical instability: inertial wave turbulence versus geostrophic turbulence[END_REF]). Yet, it is difficult to draw robust planetary conclusions from these studies, which only explored the dynamics for values of and E that were not representative of the planetary regime. Therefore, the competition between bulk and boundary mechanisms remains to be explored in the geophysically relevant regime of small forcing amplitudes 1 and small Ekman numbers E 1. To do so, note that the curvature of the boundaries (i.e. the beta effect) should be included to obtain realistic large-scale zonal flows for planetary systems. However, this effect cannot be consistently taken into account in any local Cartesian models that are commonly used in turbulence (e.g. [START_REF] Godeferd | Structure and dynamics of rotating turbulence: a review of recent experimental and numerical results[END_REF]. Consequently, we should strive considering global geometries to develop more realistic models of planetary bodies.

Apart from the interplay with bulk turbulence, the boundary layers could also become turbulent in the presence of strong enough forcings (e.g. [START_REF] Noir | An experimental and numerical study of librationally driven flow in planetary cores and subsurface oceans[END_REF][START_REF] Calkins | Axisymmetric simulations of libration-driven fluid dynamics in a spherical shell geometry[END_REF][START_REF] Sous | Friction law and turbulent properties in a laboratory Ekman boundary layer[END_REF], which could modify the boundary-driven geostrophic flows. The geostrophic shear attached to the critical latitudes should also be further characterised. For instance, a naive extrapolation of the aforementioned scaling laws would predict an amplitude for the geostrophic shear velocity of ∼ 10 -1 m.s -1 for the lunar precession (figure 17), which is an order of magnitude larger than the expected differential velocity between the lunar core and mantle [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF]. Therefore, the observed scaling laws cannot be valid in the asymptotic regime of very low Ekman numbers, as previously reported in preliminary experiments [START_REF] Morize | Experimental determination of zonal winds driven by tides[END_REF]. The intense geostrophic differential rotation at the critical latitudes could also become unstable [START_REF] Sauret | Tide-driven shear instability in planetary liquid cores[END_REF], for instance due to shear instabilities (e.g. Busse 1968a;[START_REF] Schaeffer | Quasigeostrophic model of the instabilities of the Stewartson layer in flat and depth-varying containers[END_REF], which may lead to space-filling turbulence and mixing. Moreover, the shell geometry should be further explored for more accurate planetary applications. Additional viscous effects are indeed expected due to the presence of an inner core (e.g. the reflection of inertial waves, see in [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF], such that exploring shell geometries should be further continued. Moreover, we have only validated formula (5.5) for a few librationdriven zonal flows, but it could apply to other forcings (e.g. precession) in shell geometries or possibly other geometries. For instance, ignoring the need for joining corner regions (as in [START_REF] Wedemeyer | Viscous corrections to stewartson's stability criterion(viscous corrections to stewartson theory on stability of spinning top containing liquid[END_REF], mean zonal flow can be calculated in no-slip half-spheroids (as e.g. in [START_REF] Noir | Experimental study of libration-driven zonal flows in non-axisymmetric containers[END_REF]) by summing the contribution of the plane boundary layer (i.e. half the mean zonal flow in the cylinder, given by [START_REF] Wang | Cylindrical tank of fluid oscillating about a state of steady rotation[END_REF]) and the contribution of the curved boundary (i.e. half the one of the full spheroid).

Finally, the core-mantle boundary of most planets exhibits roughness [START_REF] Narteau | On a small-scale roughness of the core-mantle boundary[END_REF][START_REF] Le Mouël | Dissipation at the core-mantle boundary on a small-scale topography[END_REF]), but scant attention has been given to the flow dynamics in the presence of small-scale topography (e.g. [START_REF] Burmann | Effects of bottom topography on the spin-up in a cylinder[END_REF]. However, our asymptotic theory could be used to get further physical insights into topographic effects for planetary applications. A small-scale azimuthal roughness could be mimicked here using the multipolar tidal-like forcing (2.4) with ω = 0 and Ω c = 0 (such that U = z R × r). The mathematical problem is tractable in the short azimuthal wavelength approximation (i.e. m 1), and we obtain the mean zonal flow f (s) = s 2(m-2) /4 → 0 when m → ∞ (the mean zonal flow driven by weak librations of a rotating sphere also vanishes in the limit ω 1, see in [START_REF] Sauret | Libration-induced mean flow in a spherical shell[END_REF]. Therefore, it appears that a small-scale azimuthal roughness is unlikely to sustain significant mean zonal flows in planetary interiors via this mechanism. Investigating this problem deserves further numerical work, as well as exploring the flows driven by other small-scale topographies.

1 h ∂ ∂q 2 sin q 2 ∞ 0 v 2 dζ + 1 T (q1) ∂ ∂φ ∞ 0 v φ dζ , ( B 1)
where the last term vanishes when calculating the mean zonal component of v.

B.2. Generic solution for the forced interior flow P U 1 0

Using the Cartesian coordinates of our frame of reference, the uniform-vorticity flow satisfying the no-penetration condition is given by (e.g. [START_REF] Noir | Precession-driven flows in non-axisymmetric ellipsoids[END_REF])

U 1 0 =   f 2 0 Q y z R -Q z y R -f 2 0 Q x z R + Q z x R f 2 0 (r pol /r eq ) 2 [Q x y R -Q y x R ]   (B 2)
with f 2 0 = 2/(1 + r 2 pol /r 2 eq ), and where the (constant

) rotation vector Q = (∇ × U 1 0 )/2 has the Cartesian components Q = Q x x R + Q y y R + Q z z R given by Q x = Q x cos(ωt) + Q x sin(ωt), ( B 3a
)

Q y = ( Q y sin ωt + Q y cos ωt), ( B 3b
)

Q z = -Q z cos(ωt). (B 3c)
This generic form encompasses all the particular cases considered in this article, and is compliant with the ansatz used below to integrate equations (2.15). For instance, the flow U 1 0 driven by longitudinal librations in the mantle frame of reference is obtained with Q =cos(ωt) z R . For latitudinal librations, the viscous flow in the mantle frame of reference is obtained from [START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF]. We have corrected a few typos in their expressions (3.29)-(3.31), which leads to

Q x = (ω 2 -f 2 )(f 2 0 -ω 2 ) -K 2 (f 2 0 + ω 2 ) (ω 2 -f 2 ) 2 + 2K 2 (f 2 + ω 2 ) + K 4 , (B 4a) Q x = Kω ω 2 + f 2 -2f 2 0 + K 2 (ω 2 -f 2 ) 2 + 2K 2 (f 2 + ω 2 ) + K 4 , ( B 4b
)

Q y = ω β bc (ω 2 -f 2 ) -K 2 (1 + f 2 /f 2 0 ) (ω 2 -f 2 ) 2 + 2K 2 (f 2 + ω 2 ) + K 4 , ( B 4c
)

Q y = K ω 2 + f 2 -2ω 2 f 2 /f 2 0 + K 2 (ω 2 -f 2 ) 2 + 2K 2 (f 2 + ω 2 ) + K 4 , ( B 4d )
with

β bc = (f 2 -f 2 0 )/f 2 0 , K = ΛE 1/2
where Λ 0 is the viscous damping factor, and with the eigenfrequency f of the spin-over mode, that is f = f 2 0 for the spheroid as in [START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF]). Yet, note the erroneous presence in their formula (3.30) of β bc in Q y , as well as that of z M instead of z I in their expression (A5). In the inviscid limit E = 0, we have Q x = Q y = 0 and, for a spheroid,

Q x = (f -ω 2 )/(ω 2 -f 2 ), Q y = ωβ/(ω 2 -f 2 ), (B 5a,b)
where β = (r 2 eqr 2 pol )/(r 2 eq + r 2 pol ), giving a diverging flow at ω = f for any spheroids (except the sphere β = 0, where the solution reduces to Q x = -1 and Q y = 0). Finally, the precession-driven flow U 1 0 can be written under this form in the mean rotating frame (with ω = 0), also called precessing frame, or in the mantle frame [START_REF] Noir | Precession-driven flows in non-axisymmetric ellipsoids[END_REF]. B.3. Calculation of H U 1 0 for longitudinal librations In their analytical study, [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] considered longitudinal librations at ω 1 in axisymmetric arbitrary containers of revolution around the z R , and obtained H U 1 0 in the form of a quasi-geostrophic flow. Their assumption ω 1 implies a steady boundary layer at leading order, whereas the boundary layer can be time dependent when ω E 1/2 (as considered in the main text), in particular when ω O(1). Since the two regimes overlap for E 1/2 ω 1, we aim to understand the transition between these two situations. To do so, we determine below the quasigeostrophic component g U 1 0 of H U 1 0 for arbitrary values of ω, extending the study of [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] to time-dependent boundary-layer flows. The slight differences between H U 1 0 and g U 1 0 will then be briefly studied by performing the exact calculation in the particular case ω > 2 (i.e. without any critical latitude).

Using cylindrical coordinates (s, φ, z) in the mean rotating frame, we consider a fluid within an arbitrary axisymmetric container g1 (s) z g2 (s), with g1 (s) 0 and g2 (s) 0. The zonal (i.e. m = 0) component zo U 1 0 of H U 1 0 can be written as

zo U 1 0 = -∇ × (Ψ 1 0 (s, z, t) φ) + V 1 0 (s, z, t) φ = [∂ z Ψ 1 0 , V 1 0 , -∂ s (s Ψ 1 0 )/s] , ( B 6)
which satisfies divergenceless condition (3.3b). Then, the curl of equation (3.3a), that is

(∂ t + Ω 0 ∂ φ ) ∇ × ( H U 1 0 ) -2 Ω 0 s ( z R • ∇) H U 1 0 = 0, ( B 7)
can be written in this case

∂ t V 1 0 = -2 ∂ z Ψ 1 0 , ∂ t ( LΨ 1 0 ) = 2 ∂ z V 1 0 , (B 8a,b) with LΨ 1 0 = (∇ 2 -s -2 )Ψ 1 0 = ∂ s [s -1 ∂ s (sΨ 1 0 )] + ∂ 2 zz Ψ 1 0 .
To obtain the quasi-geostrophic component g U 1 0 of H U 1 0 , the velocity components perpendicular to the rotation axis are assumed to be z-invariant [START_REF] Labbé | On magnetostrophic inertia-less waves in quasigeostrophic models of planetary cores[END_REF]. Then, equation (B 8) naturally reduces to the Taylor-Proudman constraint ∂ z ( g U 1 0 ) = 0. Using the axisymmetric decompositions

g U 1 0 + E 1/2 U 1 1 = -E 1/2 ∇ × (χ I φ) + V I φ, (B 9a) g u 1 0 + E 1/2 u 1 1 = -E 1/2 ∇ × (χ B φ) + V B φ, ( B 9b)
where we have anticipated that the meridional stream functions are of the order E 1/2 (as in [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF], we obtain ∂ τ V I + 2 ∂ z χ I = 0 and then

χ I = -(z/2) ∂ τ V I + χ 0 I , ( B 10)
as previously obtained by [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF]. As we follow closely the approach and the notations of [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF], we do not remind below all the intermediate steps. In the boundary layer, equations (5.7)-(5.10) of [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] are modified into (see also equation

(3.2) in Sauret & Le Dizès 2013) ∂ 2 ζζ -∂ t V B -2 n • z R ∂ ζ χ B = 0, ∂ 2 ζζ -∂ t ∂ 2 ζζ χ B + 2 n • z R ∂ ζ V B = 0, (B 11a,b
) which can be integrated by considering the ansatz e iωt for [χ I , χ B , V I , V B ], leading to

V B = s (1 -f 0 ) 4 e -λ+ζ + e -λ * -ζ e i ωt + c.c., ( B 12)
with V I = m f 0 s e i ωt , and c.c. the complex conjugate. Note that equation (B 12) can naturally be exactly recovered with equation (B 27) by considering the BC u

1 0 +U 1 0 = V 1 Σ with a non-zero g U 1 0 . From equation (B 11) we get χ B (ζ = 0)| z=gi = sgn(g i ) s(1 -f 0 ) Λi 4 e i ωt + c.c., (B 13a) Λi = 1 2 | n i • z R | λ + + λ * --i ω 1 λ + + 1 λ * - , (B 13b) with | n i • z R | = [1 + (d s gi ) 2 ] -1/2
, and where the terms λ + and λ * -involved in Λi have to be calculated with the associated n i • z R (noting d s gi the derivative of the one-variable functions g1 (s) z g2 (s) describing the container geometry). Then, the BC χ B +χ I = 0 can be written using equations (B 10) and (B 13) as

∂ τ V I = 2 [χ B (ζ = 0)| z=g2 -χ B (ζ = 0)| z=g1 ][g 2 -g1 ] -1 , ( B 14a)
χ 0 I = [g 1 χ B (ζ = 0)| z=g2 -g2 χ B (ζ = 0)| z=g1 ][g 2 -g1 ] -1 , ( B 14b)
which provides χ 0 I = m A 0 e i ωt and

f 0 = (1 + i ω S -1 E -1/2 ) -1 , S = ( Λ1 + Λ2 )/(g 2 -g1 ). (B 15a,b)
Equation (B 15) gives V I , but also χ I using equation (B 10), providing both the geostrophic part of U 1 0 and the associated component of U 1 1 . The asymptotic regime studied by [START_REF] Greenspan | On a time-dependent motion of a rotating fluid[END_REF] is recovered by using ω = 0 in S (i.e. ω 1), that is when the boundary layer can be assumed to be steady (giving e.g. S ≈ (1-s 2 ) -3/4 within the sphere in this regime). In figure 18(a), we confirm that this approximation is in excellent agreement with DNS for ω 1. By contrast, figure 18(b) shows that only the more complete theory developed here is in excellent agreement with DNS for ω O(1), as expected. We now come back to our initial question of the decrease of the bulk flow V I in the limit ω E 1/2 , for the opposite regimes ω 1 and ω 1. Since S is independent of ω in the particular regime ω 1, equation (B 15) shows that V I

V φ /( s) (×10 2 ) (a) ω = 10 -2 (b) ω = 3
Figure 18. Normalised (instantaneous) rotation rate V φ /( s) as a function of time for longitudinal librations of a sphere (E = 10 -6 and = 10 -4 ), at position (s = 0.9, q2 = π/2) in the equatorial plane. DNS are given by the solid red lines and the theory V φ /( s) = m f0 e i ωt , where f0 is given by equation (B 15a), is shown as dotted blue and dashed green lines, when calculating S with equation (B 15b) or with its approximation for ω = 0, respectively.

decreases toward 0 as V ∼ E 1/2 /ω. By contrast, S ∼ √ ω for ω 1, showing that V I decreases toward 0 as V I ∼ E/ω in this regime.

To confirm further that g U 1 0 is a good approximation of H U 1 0 , we now aim at obtaining H U 1 0 directly from equation (B 8) without the quasi-geostrophic assumption. To this end, we have adapted the calculation of [START_REF] Wang | Cylindrical tank of fluid oscillating about a state of steady rotation[END_REF]

performed in cylinders. With the ansatz (V 1 0 , Ψ 1 0 ) = ( V 1 0 , Ψ 1 0 )e iωt + c.c. and noting ω = 1 -4/ω 2 , equation (B 8) gives ∂ s [∂ s (s Ψ 1 0 )/s] + ω ∂ 2 zz Ψ 1 0 = 0, ( B 16)
whose general solution is (imposing regularity at s = 0)

Ψ 1 0 = k J 1 ( Ck s) Ãk e Ck z/ω 1/2 + Bk e -Ck z/ω 1/2 , ( B 17)
with J 1 the Bessel function of the first kind and with the constants [ Ãk , Bk , Ck ]. In the quasi-static regime ω → ∞, the leading order naturally recovers the linear z-dependency of quasi-geostrophic solution (B 10). In symmetric containers with g1 = -g 2 , we have Ψ 1 0 (z = 0) = 0 by symmetry, such that Ãk = -Bk . Equation (B 17) then reduces to

Ψ 1 0 = k Ãk J 1 ( Ck s) sinh Ck z ω1/2 . ( B 18)
The (complex-valued) constants [ Ãk , Ck ] are then fixed by the BC Ψ 1 0 + E 1/2 χ B = 0 at ζ = 0, where χ B is given by equation (B 13). Contrary to the cylinder, the properties of Fourier-Bessel series cannot be used to obtain the constants [START_REF] Wang | Cylindrical tank of fluid oscillating about a state of steady rotation[END_REF]

. Using f 0 = v φ /s and v φ = -2∂ z Ψ 1 0 /(iω) in χ B (ζ = 0), the BC reads Λ2 4 s + k ãk J 1 ( Ck s) sinh Ck z ω1/2 - i Λ2 Ck 2ω ω1/2 cosh Ck z ω1/2 = 0, ( B 19)
with ãk = Ãk /E 1/2 . When ω < 2, many terms are required in equation (B 19) to accurately fulfil this BC (due to the divergence at the critical latitude). In this case, bulk inertial modes can also be excited (e.g. [START_REF] Aldridge | Axisymmetric inertial oscillations of a fluid in a rotating spherical container[END_REF][START_REF] Zhang | Theory and modeling of rotating fluids: convection, inertial waves and precession[END_REF], and they can then constitute a better basis to describe Ψ 1 0 . By contrast, we find that only few terms are necessary for ω > 2, and there is no excitation of inertial mode. Considering for instance ω = 3 in the sphere (s = sin θ, z = cos θ) and the first three terms of Ψ 1 0 , we impose the BC at six equally spaced values of the colatitude θ ∈]0, π/2[. The numerical integration of this nonlinear system of six equations provides then the 6 constants, with ã1 ≈ 0.26 + 0.65i and C1 ≈ 0.63 -0.21i for the leading-order term (ã 2 and ã3 are respectively ∼ 10 and ∼ 10 3 times smaller, showing convergence of the series). One can then check a posteriori that BC (B 19) is verified on the whole range θ ∈ [0, π/2] with a maximum error < 2 × 10 -6 . Using these constants, we can then calculate the three components of velocity and compare with DNS. This solution naturally recovers the excellent agreement shown in figure 18(b), and even agrees slightly better with DNS as shown in figure 19(a) for s = 0.1. As can be expected, this difference is due to the weak departure of the flow from quasi-geostrophy, as shown in figure 19(b) that is in perfect agreement with our DNS. To conclude this section, note that the mean zonal component of ( H U 1 0 • ∇) H U 1 0 is non-zero, and H U 1 0 should thus a priori modify the bulk mean zonal flow when not considering the no spin-up regime ω E 1/2 .

B.4. Inertial modes excited by longitudinal librations and mean zonal flows

Longitudinal librations can excite inertial modes through the Ekman pumping u 1 1 generated by the oscillating Ekman layer (see sections 2.12 and 2.14 in [START_REF] Greenspan | The theory of rotating fluids[END_REF]). The amplitude of such forced inertial modes is of the order O( E 0 ) when the forcing frequency ω E 1/2 matches the frequency of an inertial mode [START_REF] Greenspan | The theory of rotating fluids[END_REF][START_REF] Aldridge | Axisymmetric inertial oscillations of a fluid in a rotating spherical container[END_REF][START_REF] Zhang | Theory and modeling of rotating fluids: convection, inertial waves and precession[END_REF]. This inertial mode excitation is not an inviscid resonance (as encountered e.g. with latitudinal librations, see in [START_REF] Vantieghem | Latitudinal libration driven flows in triaxial ellipsoids[END_REF], and thus its amplitude remains finite in the limit E → 0 (but vanishes when E = 0). For a sphere in longitudinal libration at the inertial mode frequency ω = ω 12 = 12/7 ≈ 1.309, [START_REF] Zhang | Theory and modeling of rotating fluids: convection, inertial waves and precession[END_REF] obtained the following flow in spherical with Ã12 = 0.034156 -i0.13641, and noting u 1 0 = P u 1 0 + H u 1 0 where P u 1 0 (resp. H u 1 0 ) is the boundary-layer flow associated with P U 1 0 (resp. H U 1 0 ). In the container frame considered in [START_REF] Zhang | Theory and modeling of rotating fluids: convection, inertial waves and precession[END_REF], the flow is mainly an apparent one, that is the oscillating solid-body rotation P U 1 0 directly related to the frame motion (there is no spin up since ω E 1/2 ). By contrast with the findings of [START_REF] Zhang | Theory and modeling of rotating fluids: convection, inertial waves and precession[END_REF], the bulk flow reduces to the inertial mode in the mean rotating frame (where P U 1 0 = 0), and is then strongly dependent on ω (similar to the findings of Aldridge & Toomre 1969, who found a way to measure this effect). Note however that P u 1 0 = 0 in the mean rotating frame (due to the oscillating boundary velocity), and this flow generates the mean zonal flow in the absence of any other bulk flows (e.g. when no inertial mode is excited with ω > 2).

Beyond the mean zonal flow obtained by assuming H U 1 0 = 0 (i.e. generated by P u 1 0 only), as mainly considered in this work, one can use equation (B 20) to calculate how this mean flow is modified by the bulk inertial mode (i.e. by the non-zero H U 1 0 and H u 1 0 ). To do so, we first note that the mean zonal component of ( H U 1 0 • ∇) H U 1 0 + 2Ω 1 c × H U 1 0 is zero. One can thus proceed exactly as in sections 3.2-3.4 to integrate equations (2.18)-(2.19). The result is shown in figure 20(a), where it is compared with DNS and with the mean zonal flow generated by P u 1 0 only. Note that the presence of the critical latitude is difficult to spot on the DNS mean zonal flow (by contrast with ω = 0.1 or ω = 1, see in figure 10b). Moreover, the disagreement suggests that considering a single forced inertial mode H U 1 0 in equations (2.18)-(2.19) is not sufficient to explain the generation of the observed jets. Contributions neglected in the theory are thus expected to be significant, such as the presence of several inertial modes in the bulk, or nonlinear interactions within the internal shear layers (as in [START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF].

To perform a more systematic comparison, we show in figure 20(b) the mean zonal flow rotation rate (at s = 0.6) given by the theory (assuming H U 1 0 = 0) and by DNS. One can observe the expected divergence at ω = 2 cos sin -1 0.6 = 1.6, due to the presence of the critical latitudes. Note also the overall rather good agreement, except near inertial mode frequencies excited by the forcing (the apparent good agreement at the eigenfrequency ω = 12/7 is coincidental and related to the choice s = 0.6, as shown in figure 20a). where A k and B k are directly obtained using BC (2.15b). The velocity V 1 Σ imposed at the spheroidal boundary depends on the problem at hand. Here, we consider that V 1 Σ can generically be written as

V 1 Σ = V m + V st + V un , V m = Qz cos(ωt) s q-1 cos(mφ) z R × r, (B 23a,b)
where V m is the multipolar and oscillating extension of the BC used by [START_REF] Suess | Viscous flow in a deformable rotating container[END_REF], who considered the particular case ω = 0 and m = 2. In equation (B 23), V un is a uniform-vorticity flow, of the form (B 2) but with the rotation vector q = (q x , qy , qz ), and the velocity V st = V sbn • V sb is given by the tangential components on Σ of the solid-body rotation V sb = Qx cos(ωt) x R × r. Since U 1 0 -V un is a uniform-vorticity flow, we assume that q is of the form of Q, that is

∇ × (U 1 0 -V un ) 2 = Q -q =   Q x cos ωt + Q x sin ωt Q y sin ωt + Q y cos ωt -( Q z + qz ) cos ωt   . ( B 24)
The last component in (B 24) allows us to reproduce the particular case m = 0 of expression (B 23b). Since non-zero Q z + qz will only be considered when m = 0, we simplify the equations by replacing Qz by Q z in expression (B 23b), with Q z = Qz + qz -Q z , and by putting Q z + qz = 0 in equation (B 24). Then, replacing the m = 1 dependency of V un + V st by the formal m dependency imposed by the ansatz of Y , equation (2.15b) gives

u 1 0 = V m + V st + V un -U 1 0 = 1 4     0 A (1) -C 1 1 -A (-1) - C -1 -1 + A (1) + C -1 1 -A (-1) + C 1 -1

B

(1)

-C 1 1 + B (-1) - C -1 -1 + B (1) + C -1 1 + B (-1) + C 1 -1     (B 25)
at the boundary q 1 = Q 1 , with C l k = exp[i(kmφ + lωt)] and

A (j) ± = i Qx a sinh 2Q 1 2 h + a h tanh(2Q 1 )[∓j(Q x ∓ Q y ) -i (Q x ± Q y )], (B 26a) B (j) ± = -Qx aT (Q1) cos q 2 + Q z s q +a tanh(2Q 1 )T (Q1) cos q 2 [(Q x ± Q y ) ∓ i j (Q x ∓ Q y )]. (B 26b)
Using BC (B 25), equation (B 22) gives Y , and u 1 0 reads u 1 0 • q 1 = 0, (B 27a) -± iB

u 1 0 • q 2 = 1 8 (A (
(j) - = a tanh(2Q 1 )[-i(Q x -Q y ) -j(Q x + Q y )]( h ∓ T (Q1) cos q 2 )
+i Qx a[sinh(2Q 1 )/(2 h) ∓ T (Q1) cos q 2 ] ± i Q z s q (B 28a)

B (j) ± = A (j) + ± iB (j) + = a tanh(2Q 1 )[-i(Q x + Q y ) + j(Q x -Q y )]( h ∓ T (Q1) cos q 2 )
+i Qx a[sinh(2Q 1 )/(2 h) ∓ T (Q1) cos q 2 ] ± i Q z s q , (B 28b) and

λ ± = [1 + i sgn(γ 1 ± γ + )] |γ 1 ± γ + |, κ ± = [1 + i sgn(γ 1 ± γ -)] |γ 1 ± γ -|, (B 29a,b)
where we have introduced γ ± = (mΩ 0 ± ω)/2. Note also the following identities λ ± = i sgn(γ 1 ± γ + )λ * ± , λ ± /λ * ± = -λ * ± /λ ± , and λ ± /λ * ± 3 = λ * ± /λ 3 ± (see also in Busse 2010), with similar identities for κ ± . In order to ease the cumbersome calculation of the mean zonal flow, it is also useful to already note that ∂λ ± ∂q 2 = -sgn(γ 1 ±γ + ) γ 1 T (q1) 2 tan q 2 h2 λ * ± , ∂κ ± ∂q 2 = -sgn(γ 1 ±γ -) γ 1 T (q1) 2 tan q 2 h2 κ * ± , (B 30a,b) which also gives the derivative of the complex conjugates by swapping λ ± and λ * ± (as well as κ ± and κ * ± ). The boundary-layer flow u 1 0 is actually generated by the differential velocity u 1 0 =

Figure 2 .

 2 Figure 2. Left: Spheroidal geometry of the forced problem. Right: equatorial (circular) section.

Figure 3 .

 3 Figure 3. Azimuthal component of the time-averaged zonal flow V φ / 2 as a function of cylindrical radius s in the equatorial plane z0 = 0. DNS in a sphere with E = 10 -4 , = 10 -2 and forcing frequency ω = π, using Nek5000 (red solid curves) and xshells (blue dashed curves). (a) Longitudinal librations. (b) Latitudinal librations.

Figure 4 .

 4 Figure 4. Rotation rate of the mean zonal flow, as a function of the cylindrical radius s, in rotating spheres subject to weak longitudinal librations (with Ω0 = 0, m = 0). (a) Our analytical results (solid coloured curves) are compared with the theoretical results (black dashed curves) of Sauret & Le Dizès (2013). (b). Comparison between the theory (black dashed curves) and DNS (solid coloured curves) at E = 10 -7 and = 10 -4 when ω 2.

Figure 5 .

 5 Figure5. Rotation rate of the mean zonal flow, as a function of the cylindrical radius s, in rotating spheroids subject to longitudinal librations with ω = π. (a) Comparison between DNS with E = 10 -5 and = 5 × 10 -4 (solid coloured curves) and theory (dashed black curves). Numerical profiles have been truncated for s 0.05, because there are not enough grid points near the centre to get reliable numerical averages with formula (2.8) for the rotation rate. (b) Comparison of the theoretical rotation rate between the sphere, oblate spheroids (r pol /req = 0.2), and prolate spheroids (r pol /req = 10). Horizontal dashed blue curve indicates the theoretical profile for the cylinder[START_REF] Wang | Cylindrical tank of fluid oscillating about a state of steady rotation[END_REF][START_REF] Sauret | Fluid flows in a librating cylinder[END_REF].

Figure 6 .

 6 Figure 6. Rotation rate of the mean zonal flow, as a function of the cylindrical radius s, for a rotating sphere subject to weak latitudinal librations. (a) Comparison between theory (black dashed curves) and DNS (coloured solid curves) at E = 10 -7 and = 10 -4 in the regime ω 2. (a) Theoretical zonal flows when ω < 2 with divergent flows at the critical latitudes.

Figure 7 .

 7 Figure 7. Rotation rate of the mean zonal flow, as a function of the cylindrical radius s, for a precessing sphere. (a) Comparison of theoretical profiles between our theory and Busse's predictions (black squares, extracted from figure 1 in Busse 1968b). (b) Comparison between our theory (black dashed curve) and DNS (solid curves) for different values of Po with γ = Po/E 1/2 .

Figure 8 .

 8 Figure 8. Mean zonal flows driven by precession in spheres. (a) Value at s = 0 of f , as a function of γ = Po/E 1/2 in DNS (symbols). Dashed blue line shows a slope -0.15. (b) Rotation rate of the mean zonal flow, as a function of the cylindrical radius s. Comparison between theory (black dashed curve) and DNS (solid coloured curves) choosing values of γ such that f (0) = 0.

Figure 9 .

 9 Figure9. Rotation rate of the mean zonal flows driven by the tidal-like forcing of[START_REF] Suess | Viscous flow in a deformable rotating container[END_REF], obtained with Ωc = ω = 0 and m = 2 in expression (2.4). (a) Thick solid red curve: present theory. Solid orange curve: experimental results of[START_REF] Suess | Viscous flow in a deformable rotating container[END_REF]. Dashed blue curve: erroneous theory of[START_REF] Suess | Viscous flow in a deformable rotating container[END_REF]. Inset shows the azimuthal velocity sf . (b) Various contributions to f (solid red curve), made of the sum of two contributions: the one due to uB (dotted curve), which was mistakenly considered to be f by[START_REF] Suess | Viscous flow in a deformable rotating container[END_REF], and the one due to uC (dashed curve).

Figure 10 .

 10 Figure10. Geostrophic velocity sf (or angular velocity f in the insets) of the mean zonal flow, as a function of the cylindrical radius s, for a sphere subject to (a) latitudinal librations at ω = 0.1 and (b) longitudinal librations for two ω. Comparison between theory (black dashed curves) and DNS (solid coloured curves) at = 10 -4 . The curve at E = 10 -7 is indistinguishable from the one obtained with = 10 -5 (not shown). In (b), DNS have been performed at E = 10 -7 , and we have checked that the numerical profile f at ω = 1 is unchanged for = 10 -3 and = 10 -6 (not shown). Insets show that f (0) remains finite in DNS.

Figure 11 .

 11 Figure 11. Mean zonal flows in rotating spheroids subject to longitudinal librations in the presence of critical latitudes (ω = π/2 2). (a) Comparison of the mean zonal velocity sf (s)given by the theory (black dotted curve for the sphere, black dashed curve for r pol /req = 0.8) and DNS (coloured solid curves) at E = 10 -5 and = 10 -4 . (b) Comparison of the theoretical rotation rate between the sphere, an oblate (r pol /req = 0.2) or prolate (r pol /req = 4) spheroid.

Figure 12 .

 12 Figure 12. Geostrophic shear associated with the critical latitudes of longitudinal librations with ω 2. (a) Comparison between theory and DNS (Lin & Noir 2020) in a spherical shell for libration-driven zonal flows computed with = 10 -2 . The grey area indicates the tangential cylinder s 0.1 associated with the inner core in Lin & Noir (2020). (b) Schematic regime diagram for the evolution of δug and δs, as a function of ω and E. Pink area indicates the regime dominated by the scalings of the Ekman boundary layer. In the blue area the dominant scalings are given by (5.1). Inset illustrates the geostrophic shear centred on sc in a meridional section, where the blue area represents the Ekman boundary layer of thickness E 1/2 .

r

  pol /r eq = 0.8 r pol /r eq = 1.2 (a) Spherical shell (r pol /req = 1) (b) Spheroidal shells (r pol /req = 1)

Figure 15 .

 15 Figure15. Normalised geostrophic mean velocity V φ / 2 at various heights z = √ 1s 2 in shells with ηeq = 0.5. Black dashed curves indicate the theory in the shell. The DNS profiles computed at z = z0 for s ηeq and at z0 = 0 for s > ηeq, giving a single profile for each DNS. Grey area shows the tangent cylinder s ηeq. (a) Non-homoeoidal shells subject to longitudinal librations (solid coloured curves), with a spherical inner boundary (i.e. r in pol = r in eq ) and a spheroidal outer boundary (i.e. r pol = req). DNS performed at E = 2.5 × 10 -6 , ω = ωin = π, = in = 5 × 10 -4 . (b) Homoeoidal shells with η pol = ηeq = 0.5. Case 1: DNS at E = 10 -6 in a spherical shell subject to latitudinal librations at outer boundary (with ω = 3 and = 10 -6 ), and to longitudinal librations at inner boundary (with ωin = 2, in/ = 2). Case 2: Spheroidal shell with r pol /req = 0.8 subject to longitudinal librations. DNS at E = 5 × 10 -6 with = 5 × 10 -4 and ω = π at outer boundary, and with in/ = ωin/ω = 2 at inner boundary.

Figure 16 .

 16 Figure 16. Normalised geostrophic mean velocity in a spherical shell (η pol = ηeq = 0.35) subject to longitudinal librations at ω = √ 2. In panel (a), in = √ 2/100 and = 0, whereas, in panel (b), in = 0 and = √ 2/100. Legend in panel (b) is identical to the one in panel (a). Coloured curves have been reproduced from figures 13(a) and 17(a) in[START_REF] Lin | Libration-driven inertial waves and mean zonal flows in spherical shells[END_REF], where V φ is computed as in formula (2.8), with zmax = 1, but only considering only the m = 0 component. Theory is given by formula (5.5). Grey area indicate the tangent cylinder s 0.35.

Figure 17 .

 17 Figure 17. Schematic regime diagrams for the existence of mean zonal flows driven by mechanical forcings (planetary values estimated from Noir et al. 2009; Cébron et al. 2012; Lin et al. 2015; Vantieghem et al. 2015). Red circles: precession. Blue squares: longitudinal librations. Black triangles: latitudinal librations. Magenta stars: tides. (a) Transition between laminar and turbulent boundary layers (BL). Transition occurs when ∼ K E 1/2 (with the typical values K = 20 -150 shown by the grey area, K = 55 by the dashed dotted line). (b) Competition between bulk and boundary driven generation of mean zonal flows. Input Rossby number Ro = O(β) for tidal forcing or Ro = O( β) for precession and libration forcings, where β is the typical boundary (equatorial or polar) ellipticity. Instabilities and bulk turbulence (hatched area) onsets when Ro K E 1/2 (with the values K = 1 -10 shown by the grey area). Thick dashed line shows Ro = E 1/4 /k, with k = 10.

Figure 19 .

 19 Figure 19. (a) Same as in Figure 18(b) but at s = 0.1 and showing the solution obtained from equation (B 18) in black dashed line. (b) Snapshot v φ = -2∂z Ψ 1 0 /(iω) in a meridional section of the theoretical azimuthal flow obtained from equation (B 18). We obtain an excellent agreement with the analogous snapshot from DNS (not shown).

Figure 20 .

 20 Figure20. Mean zonal flow driven by longitudinal librations of a sphere, from theory and DNS (E = 10 -7 , = 10 -4 ). (a) Geostrophic velocity sf , as a function of the cylindrical radius s, at the inertial mode eigenfrequency ω = 12/7[START_REF] Zhang | Theory and modeling of rotating fluids: convection, inertial waves and precession[END_REF]. DNS (solid line) does not agree well with the theory, whether H U 1 0 is taken into account (dashed line) using equation (B 20) or not (dotted line). (b) Rotation rate f at s = 0.6, as a function of the libration frequency ω, given by the theory (solid line) assuming H U 1 0 = 0. Dashed lines show the four least damped inertial modes in figure3of[START_REF] Aldridge | Axisymmetric inertial oscillations of a fluid in a rotating spherical container[END_REF], including ω = 12/7.

B. 5 .

 5 First-order boundary-layer flowIntegrating the coupled ordinary differential equations (3.4), we obtainY k = A 1 e α k+ ζ + A 2 e -α k+ ζ + A 3 e α k-ζ + A 4 e -α k-ζ -i A 1 e α k+ ζ + A 2 e -α k+ ζ + i A 3 e α k-ζ + A 4 e -α k-ζ (B 21a) with α k± = (1 + is ± ) |γ 1 ± γ k |, and s ± = sgn(γ 1 ± γ k ). The boundary layer velocity vanishes when ζ → ∞, imposing A 1 = A 3 = 0. We thus finally obtainY = k A k e -α k+ ζ + B k e -α k-ζ e i(ω k t+m k φ) -i A k e -α k+ ζ -B k e -α k-ζ e i(ω k t+m k φ) , (B 22)

† Email address for correspondence: david.cebron@univ-grenoble-alpes.fr

Acknowledgements. J.V. and D.C. acknowledge Dr Benjamin Favier for sharing his mapping to model non-homoeoidal shells in Nek5000. D.C. acknowledges Dr Loïc Huder for his expert support in Python.

Funding. This work received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme via the theia project (grant agreement no. 847433). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). The xshells code received funding from the European Union's Horizon 2020 research and innovation programme under the ChEESE project, grant agreement no. 823844. This work was granted access to the HPC resources of TGCC and CINES under allocation A0080407382 attributed by GENCI (Grand Equipement National de Calcul Intensif).

Declaration of interests. The authors report no conflict of interest.

Data availability statement.

The open-source codes xshells and Nek5000 are available at https://nschaeff.bitbucket.io/xshells and https://nek5000.mcs.anl.gov/, respectively. The analytical calculations have been checked using Maple software.

Author ORCID. D.C., 0000-0002-3579-8281; J.V., 0000-0002-3654-6633; N.S., 0000-0001-5206-3394; A.S., 0000-0001-7874-5983

Author contributions. The theory was initiated by A.S. during his PhD. D.C then corrected his calculations during the first COVID-19 lockdown, and extended them to account for other forcings and geometries. J.V. conducted the DNS using Nek5000. D.C. and N.S. performed the DNS of libration-driven flows with xshells, whereas the DNS of precession were conducted by N.S. and A.B. The discussion was led by J.V. and D.C., who both drafted the manuscript. N.S and A.S proof-checked the article, and the authors gave final approval for submission.

Appendix A. Central regularity conditions with finite differences

We detail here how the central conditions is implemented with finite differences in xshells. We expand the velocity field V onto the set of spherical harmonics Y m l using the poloidal-toroidal decomposition in spherical coordinates (r, θ, φ)

where [P m l (r), T m l (r)] are the poloidal-toroidal radial scalars. These scalars must satisfy regularity conditions at the centre for the vector field to be regular and infinitely differentiable. To do this, the two scalars and the Cartesian components of the velocity fields must behave like monomials in the Cartesian coordinates (x, y, z). This is ensured by expanding [P m l (r), T m l (r)] in the regular form (e.g. Dudley & James 1989)

where [A j , B j ] are unknown coefficients. At the centre r = 0, the (l, m) components for the velocity (A 1) then reduce to

Within xshells, the poloidal-decomposition is implemented using vector spherical harmonics that depend on the radial scalar l(l + 1)P m l /r, the spheroidal scalar S m l = (1/r) ∂ r (rP m l ) and the toroidal scalar T m l . Therefore, we obtain from (A 3) the following regularity conditions for a finite difference scheme

This allows a non-zero velocity at the centre, corresponding to a flow going through it. These regularity conditions apply to the velocity, but also to the vorticity W m l , which is related to the velocity

Applying the same reasoning as above leads to condition (A 4) with P m l = 0 replaced by T m l = 0 and T m l = 0 replaced by ∆P m l = 0. For l = 1, this also leads to ∂ rr P m 1 | r=0 = 0. This vanishing second-order derivative of P m 1 ensures that the error for the 2-point finite difference approximation of S m 1 (r = 0) is of order 2

The above conditions are actually simpler to implement with finite differences than with some spectral descriptions in radius (e.g. see the discussion in [START_REF] Livermore | A comparison of numerical schemes to solve the magnetic induction eigenvalue problem in a spherical geometry[END_REF].

In addition, to avoid a stringent restriction on the time-step size, the spherical harmonic expansion is truncated near the centre at lower degrees

where l max is the maximum spherical harmonic degree in the DNS, and α = 0.05 is found to be an appropriate parameter to avoid spurious numerical errors near r = 0 while allowing large enough time steps. In practice, our resolution for all the DNS ensured that the truncation quickly jumps to l 6 at the second radial point (not shown), ensuring a sufficient numerical resolution. Finally, the above numerical implementation is accurate enough to determine the values of the mean flow rotation rate f near the centre. Indeed, defining

with the cylindrical radius s = r sin θ, we obtain the value of the rotation rate in the equatorial plane θ = π/2 (i.e. z 0 = 0) from the toroidal component T 0 1 (r)

Expression (A 10) shows that, as expected, the rotation rate is regular on the axis of rotation (and more generally in coreless geometries, e.g. see [START_REF] Lewis | Physical constraints on the coefficients of Fourier expansions in cylindrical coordinates[END_REF], with a possible non-zero value ∝ ∂ r T 0 1 | r=0 due to the l = 1 spherical harmonic in full spheres. The error for the 2-point finite difference approximation of f (r = 0) ∝ ∂ r T 0 1 (r = 0) is guaranteed to be of order 1

We took special care to refine the grid near r = 0 to ensure the reported values for f (r = 0) are meaningful, although they are less accurate than the values away from 0.

Appendix B. Details on the theoretical calculations

In this appendix, we first provide some useful formulas in the spheroidal coordinates used in this work (section B.1). We then provide details on the calculation of the firstorder flows in the bulk (sections B.2-B.4) and in the boundary layer (section B.5).

B.1. Miscellaneous formula in our spheroidal coordinates

The ratio of the spheroid axes is r pol /r eq = T (Q1) /T (Q1) . For oblate and prolate spheroidal coordinates, it gives respectively r pol /r eq = tanh Q 1 and r pol /r eq = coth Q 1 . Note that we only have n = q 1 at the boundary q 1 = Q 1 . Using these coordinates, the leading-order bulk flow is then Ω 0 z R × r = Ω 0 s φ, where s = (x 2 + y 2 ) 1/2 = aT (q1) sin q 2 is the cylindrical radius. Finally, the iso-surface for q 1 becomes more and more spherical when q 1 becomes large, with a spherical radius r given by a exp(q 1 ) ≈ 2r.

Note also that, for every boundary layer-flow v of components v = (0, v 2 , v φ ) in our spheroidal coordinates, we have

Σ -U 1 0 at the boundary, and can thus be the same in various frames of reference. Considering for instance longitudinal librations (m = 0), we have Qz = 1 and Q z = 0 in the mean rotating frame (i.e. an oscillating boundary with a zero basic flow), which gives Q z = 1. By comparison, we have Q z = 1 and Qz = 0 in the wall frame (oscillating basic flow with a zero boundary velocity), which also gives Q z = 1 and thus u 1 0 is formally the same in both frames of reference.

As illustrating examples, we give below the expression of u 1 0 in few particular simpler cases. Considering for instance the latitudinal librations of a sphere in the mean rotating frame of reference, V 1 Σ = V un = cos(ωt) x R × r can be imposed using Q x = -1 and m = 1, discarding the other possible contributions to

Then, in the spherical geometry limit q 1 → ∞, equation (B 27) reduces to u 1 0r = 0, (B 31a)

+ e -λ+ζ + A

(1)

-e -λ-ζ )e -i(φ+ωt)

+ e -κ+ζ + A

(1)

-e -κ-ζ )e -i(φ-ωt) (B 31b)

+ e -λ+ζ -A

(1)

+ e -κ+ζ -A

(1)

where the coordinates (q 1 , q 2 , φ) are naturally mapped on the usual spherical coordinates (r, θ, φ) used here, with A

± = i r (1∓ cos θ), and where equation (B 29) can be simplified using γ ± = ±ω/2 and γ 1 = cos θ.

Considering now the multipolar tidal-like forcing of a spheroid, V 1 Σ = cos(ωt) z R × r can be imposed using Q z = 1, discarding the other possible contributions to Equation (B 32) can even be further simplified in the two following particular cases of interest, that are (i) longitudinal librations, where m = 0 and thus λ ± = κ ∓ , and (ii) the steady tidal-like forcing, where ω = 0 and thus λ ± = κ ± .