key: cord-0336353-obnvavl3
authors: C'ebron, David; Vidal, J'er'emie; Schaeffer, Nathanael; Borderies, Antonin; Sauret, Alban
title: Mean zonal flows induced by weak mechanical forcings in rotating spheroids
date: 2021-03-18
journal: nan
DOI: 10.1017/jfm.2021.220
sha: 95dfccff749cc1c7d32a580631fbe88a9d6988e8
doc_id: 336353
cord_uid: obnvavl3

The generation of mean flows is a long-standing issue in rotating fluids. Motivated by planetary objects, we consider here a rapidly rotating fluid-filled spheroid, which is subject to weak perturbations of either the boundary (e.g. tides) or the rotation vector (e.g. in direction by precession, or in magnitude by longitudinal librations). Using boundary-layer theory, we determine the mean zonal flows generated by nonlinear interactions within the viscous Ekman layer. These flows are of interest because they survive in the relevant planetary regime of both vanishing forcings and viscous effects. We extend the theory to take into account (i) the combination of spatial and temporal perturbations, providing new mechanically driven zonal flows (e.g. driven by latitudinal librations), and (ii) the spheroidal geometry relevant for planetary bodies. Wherever possible, our analytical predictions are validated with direct numerical simulations. The theoretical solutions are in good quantitative agreement with the simulations, with expected discrepancies (zonal jets) in the presence of inertial waves generated at the critical latitudes (as for precession). Moreover, we find that the mean zonal flows can be strongly affected in spheroids. Guided by planetary applications, we also revisit the scaling laws for the geostrophic shear layers at the critical latitudes, and the influence of a solid inner core.

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. Guervilly et al. 2014; Godeferd & Moisy 2015) and planetary core flows (e.g. Aubert 2005; Schaeffer et al. 2017; Monville et al. 2019) . 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. Roberts & Aurnou 2012) , which drives the long-term dynamical evolution of planetary bodies. Moreover, mean flows could be unstable in the rapidly rotating regime (e.g. Sauret et al. 2014; Favier et al. 2014) , which could sustain space-filling turbulence and † Email address for correspondence: david.cebron@univ-grenoble-alpes.fr arXiv:2103.10260v3 [physics.flu-dyn] 24 Sep 2021 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. Christensen 2002; Aubert et al. 2002) or waves. Rapidly rotating fluids are indeed characterised by the ubiquitous presence of inertial waves (e.g. Zhang & Liao 2017) , whose restoring force is the Coriolis force. However, Greenspan (1969) 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. Newell 1969; Smith & Waleffe 1999) or wave-induced secondary instabilities (e.g. Kerswell 1999; Brunet et al. 2020) . 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 (Busse 1970) , and also strongly modifies the geostrophic flows (Greenspan 1968 ). 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 . 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. Aldridge & Toomre 1969; Noir et al. 2001a; , 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 Noir et al. , 2012 Lin & Noir 2020 ).

Rotating flows are usually characterised by the Ekman number E, which compares viscous to rotational effects. As outlined in Busse (1968b Busse ( , 2010 , 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 (Sauret 2015) , and the theory has been convincingly compared with experiments (Wang 1970 ) and simulations (Sauret et al. Experiments (Sauret et al., 2010) DNS (Sauret et al., 2010) Theory (Busse, 2010) Theory (Sauret & Le Dizès, 2013) 0 0.2 0.4 0.6 0.8 1

Experiments (Malkus, 1968) DNS (Noir et al., 2001) Theory (Busse, 1968) (a) (b) 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 (Busse 2010; Sauret & Le Dizès 2013) . As shown in figure 1(a), they can both explain the experimental results of Sauret et al. (2010) . Concerning precession, experiments (Malkus 1968) 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 (Suess 1971) 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. Kerswell 1995; Kida 2011 ) not taken into account in the theory (as in Busse 1968b; . 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. Calkins et al. 2010) . The variations of the geostrophic shear amplitude with the Ekman number are however still disputed (Noir et al. 2001b; Lin & Noir 2020) , 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 Busse , 2010 . 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 (Greenspan 1968 ) 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. Busse 2010; Sauret & Le Dizès 2013) , 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. Greenspan 1968 ). 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; Lin & Noir 2020) . 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.

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 Σ

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

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

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.

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 Greenspan 1968 , 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 (Lewis & Bellan 1990) , and to consider multipolar flows (e.g. Cébron et al. 2014; Sauret 2015) , we must consider q = |m − 1|. In expression (2.4), the case m = 2, Ω c = ω = 0 has been considered in Suess (1971) , which is extended here to account for both multipolar deformations and oscillations at the frequency ω (e.g. Sauret & Le Dizès 2013).

(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. Favier et al. 2015) with Ω 0 = 0, Ω c = [1+ cos(ωt)] z R and V Σ = 0, (ii) in the mean rotating frame of reference (e.g. Busse 2010) 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 Vantieghem et al. 2015 )

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 )

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. Cébron et al. 2019) by

and

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

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. Schaeffer et al. 2017) , 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 (Schaeffer 2013) , 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 (Marti et al. 2014) , 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 Ascher et al. 1995) . 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. Fischer et al. 2007) , 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.

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 (Favier et al. 2015 )

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 |z − z 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.

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

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. Greenspan 1968 ). Then, the mass conservation equation reduces to its usual boundary layer approximation (see p. 25 in Greenspan 1968 )

(2.10)

To solve the BLT equations, we use asymptotic theory with the small forcing amplitude 1.

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. Kerswell 1995) , although they can modify the mean zonal flows (as previously found in cylinders by Sauret et al. 2012 , 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

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

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

the first-order interior flow equations at the order E 0 are then

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

with the linear operator Lu 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 (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

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

with the BC U 2 0 · n = 0 on Σ. The boundary-layer equations are

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 .

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

which givesh = sinh 2 q 1 + cos 2 q 2 andh = cosh 2 q 1 − cos 2 q 2 for oblate and prolate spheroidal coordinates, respectively. Note that we recover the usual spherical case with q 1 → ∞, giving for instance ah → a exp(q 1 )/2 ≈ r or aT (q1) → r, with the spherical radius r (see e.g. Schmitt & Jault 2004) . 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.

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.

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 (Vantieghem et al. 2015) , precession (e.g. Noir & Cébron 2013) , 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

which has to be integrated together with the BC 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

encompass the various sign possibilities. Equation (2.15a) then reduces to

where the anti-symmetric matrix M reads (noting γ k = (ω k + m k Ω 0 )/2)

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

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. Aldridge & Toomre 1969; ).

Using the decomposition U 2 0 = P U 2 0 + H U 2 0 , the average of equations (2.18) gives

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)

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

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

which is actually the Taylor-Proudman theorem (Greenspan 1968) . 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

Finally, integrating equation (3.12) between ζ = 0 and ζ → +∞ yields (e.g. Busse 2010)

which is used to determine the unknown function f (s).

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

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

For tidal forcing, Suess (1971) 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 |)

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

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

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

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

and where the right-hand side is given by

(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

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.

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

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

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

which gives the axisymmetric mean zonal flow through equation (3.10). 

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 (Greenspan 1968 ). 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.

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 Busse (2010) in the limit of vanishing libration frequency ω → 0. Using a mathematical description in terms of a stream function, Sauret & Le Dizès (2013) 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 , in the cylindrical geometry). We obtain an excellent quantitative agreement with the results of Sauret & Le Dizès (2013) 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 Busse (2010), we aim at comparing our results with Busse's theoretical zonal flow, which surprisingly differs from the one obtained by Sauret & Le Dizès (2013) for the full sphere librating at ω → 0. Indeed, Busse (2010) which were illustrated in figure 1(a). The two profiles are indistinguishable near the rotation axis, and are actually in overall good agreement with the experimental and numerical results of Sauret et al. (2010) for small but finite values of ω 1 Nevertheless, as already noticed by Sauret & Le Dizès (2013) , 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 Busse (2010) . Actually, our asymptotic theory follows closely Busse (2010), but our results are in exact agreement with the zonal flow profile of Sauret & Le Dizès (2013) 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 Busse (2010) . 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 Busse (2010) 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 Sauret & Le Dizès (2013) 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 Vidal et al. 2020) 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 Grannan et al. 2017) . 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 (Wang 1970; Sauret et al. 2012) , 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.

We now consider the flows driven by latitudinal librations, which have only received scant attention so far (Chan et al. 2011; Zhang et al. 2012; Vantieghem et al. 2015) . In particular, the mean zonal flows have only been computed numerically at moderate values of E (Chan et al. 2011) , 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. Kida 2011) 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.

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 (Malkus 1968 ) 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

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 

with γ = P o /E 1/2 and where λ = λ r + iλ i is the spin-over damping factor given by (Hollerbach & Kerswell 1995; Noir et al. 2001b )

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

where V φ (z = 0) is the azimuthal velocity in the equatorial plane z = 0 of the fluid frame of reference considered in this section, ω f = ||ω f || is the norm of the fluid rotation vector ω f given by equation (4.2), 2 = Ω 2 0 − ω 2 f , and s is the cylindrical radius measured from the fluid rotation axis. A typical DNS is illustrated in figure 7(b) . We observe that the axial value f (0) is generally non zero in the DNS, which agrees with previous findings of Noir et al. (2001b) , but the numerical profile is in good agreement with the theory far 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. 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 ).

Revisiting the tidal-like forcing of Suess (1971) We finally consider the tidal-like forcing, first considered analytically by Suess (1971) , 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. Suess (1971) 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 Suess (1971) . Suess (1971) . Dashed blue curve: erroneous theory of Suess (1971) . 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 Suess (1971) , 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 Kerswell 2002) . 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. Cébron et al. 2012) . Since the experimental results of Suess (1971) 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).

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. By analogy, more complicated flow structures are also expected for libration-driven flows when ω 2, as previously reported in cylinders ) and spheres (Lin & Noir 2020) . 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 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).

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. Lin & Noir (2020) 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), Lin & Noir (2020) proposed the scaling laws δu g / 2 ∝ E −1/10 , δs ∝ E 1/5 , (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 (Kerswell 1995; 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

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

In the opposite regime ω E 1/10 , the relevant scaling laws should be (5.1) as proposed by Lin & Noir (2020) . In order to validate these theoretical considerations, we have 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 . As outlined in Lin & Noir (2020) , the origin of the different scaling laws between libration and precession remains puzzling.

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. Calkins et al. 2010; Sauret & Le Dizès 2013; Lin & Noir 2020) .

We first consider the case where the critical latitudes and inertial waves are absent. This regime has been theoretically investigated for longitudinal librations in Sauret & Le Dizès (2013) , 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 Stewartson 1966) . Here, we consider a (possibly nonhomoeoidal) spheroidal shell, where the inner boundary is spheroidal, with respectively . Normalised mean zonal velocity V φ / 2 at various dimensionless heights z = √ 1 − s 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 Sauret & Le Dizès (2013) , we find that the rotation rate of the mean zonal flow is given in dimensionless form by

for s < r in eq , 2 f cmb sp for s > r in eq ,

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 Sauret & Le Dizès (2013) 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 Sauret & Le Dizès (2013) . 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 (Sauret & Le Dizès 2013) . 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 in Sauret & Le Dizès 2013) . 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 Lin & Noir (2020) , 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. Rieutord et al. 2001) . 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. 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 Sauret et al. 2012) , 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.

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 (Noir et al. 2009; Calkins et al. 2010) or local shear instabilities (e.g. Lorenzani & Tilgner 2001) . The transition between laminar and turbulent Ekman boundary layers would here occur when ∼ K E 1/2 , where K is a numerical prefactor. The first boundary-layer instabilities would occur when K = 20 − 55 (e.g. Lorenzani & Tilgner 2001; Noir et al. 2009; Calkins et al. 2010) , and fully turbulent boundary layers are expected when K 150 (e.g. Caldwell & Van Atta 1970; Sous et al. 2013) . 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 Lin et al. 2015; Horimoto et al. 2020) .

However, the Ekman boundary layers are expected to become turbulent before the onset of such instabilities when E 1 (see figure 6 in Cébron et al. 2019 ). 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. Favier et al. 2015; Grannan et al. 2017) . 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 (Lemasquerier et al. 2017) . In ellipsoids, we have at leading order Ro ∼ β for tides (e.g. Grannan et al. 2017; , and Ro ∼ β for topographic precession (Kerswell 1993) or libration forcings (Vantieghem et al. 2015; Vidal et al. 2019) , 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. Kerswell 1999; Brunet et al. 2020; Le Reun et al. 2020) , 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 Kerswell 1999; Le Reun et al. 2019) . 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.

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. Christensen 2002; Aubert et al. 2002) . 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 (Busse 2010; Calkins et al. 2010; Sauret et al. 2010; Sauret & Le Dizès 2013; Lin & Noir 2020 ). 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 Su et al. 2020; Vidal et al. 2020 ) 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. Kerswell 1995) 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 (Lin & Noir 2020 ). 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 Sauret & Le Dizès 2013) .

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. Kerswell 2002; Lin et al. 2015; Vantieghem et al. 2015; Nobili et al. 2021) . Only a few experimental or numerical works have hitherto studied mean zonal flows in the presence of bulk turbulence (e.g. Favier et al. 2015; Grannan et al. 2017; Le Reun et al. 2019 ). 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. Godeferd & Moisy 2015) . 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. Noir et al. 2009; Calkins et al. 2010; Sous et al. 2013) , 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 (Williams et al. 2001) . Therefore, the observed scaling laws cannot be valid in the asymptotic regime of very low Ekman numbers, as previously reported in preliminary experiments . The intense geostrophic differential rotation at the critical latitudes could also become unstable (Sauret et al. 2014 ), for instance due to shear instabilities (e.g. Busse 1968a; Schaeffer & Cardin 2005) , 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 Lin & Noir 2020) , 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 Wedemeyer 1966) , mean zonal flow can be calculated in no-slip half-spheroids (as e.g. in Noir et al. 2012) by summing the contribution of the plane boundary layer (i.e. half the mean zonal flow in the cylinder, given by Wang 1970) 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 (Narteau et al. 2001; Le Mouël et al. 2006 ), but scant attention has been given to the flow dynamics in the presence of small-scale topography (e.g. Burmann & Noir 2018) . 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 Sauret & Le Dizès 2013) . 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.

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. 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 Livermore & Jackson 2005) .

In addition, to avoid a stringent restriction on the time-step size, the spherical harmonic expansion is truncated near the centre at lower degrees l tr (r) = l max r α max(r) + 1, (A 8) 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 Lewis & Bellan 1990) , 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.

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

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

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. Noir & Cébron 2013 )

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

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 Vantieghem et al. (2015) . We have corrected a few typos in their expressions (3.29)-(3.31), which leads to

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 (Vantieghem et al. 2015 ). 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,

where β = (r 2 eq − r 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 (Noir & Cébron 2013 ).

B.3. Calculation of H U 1 0 for longitudinal librations In their analytical study, Greenspan & Howard (1963) 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 Greenspan & Howard (1963) 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 containerg 1 (s) z g 2 (s), withg 1 (s) 0 and g 2 (s) 0. The zonal (i.e. m = 0) component zo U 1 0 of H U 1 0 can be written as

which satisfies divergenceless condition (3.3b). Then, the curl of equation (3.3a), that is

can be written in this case

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 (Labbé et al. 2015) . Then, equation (B 8) naturally reduces to the Taylor-Proudman constraint ∂ z ( g U 1 0 ) = 0. Using the axisymmetric decompositions

where we have anticipated that the meridional stream functions are of the order E 1/2 (as in Greenspan & Howard 1963) , we obtain ∂ τ V I + 2 ∂ z χ I = 0 and then

as previously obtained by Greenspan & Howard (1963) . As we follow closely the approach and the notations of Greenspan & Howard (1963) , we do not remind below all the intermediate steps.

In the boundary layer, equations (5.7)-(5.10) of Greenspan & Howard (1963) are modified into (see also equation (3.2) in Sauret & Le Dizès 2013)

) which can be integrated by considering the ansatz e iωt for [χ I , χ B , V I , V B ], leading to 

with | n i · z R | = [1 + (d sgi ) 2 ] −1/2 , and where the terms λ + and λ * − involved inΛ i have to be calculated with the associated n i · z R (noting d sgi the derivative of the one-variable functionsg 1 (s) z g 2 (s) describing the container geometry). Then, the BC χ B +χ I = 0 can be written using equations (B 10) and (B 13) as

which provides χ 0 I = m A 0 e i ωt and

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 Greenspan & Howard (1963) 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 ( decreases toward 0 as V I ∼ 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 Wang (1970) 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)

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

The (complex-valued) constants [Ã k ,C k ] 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 (Wang 1970 ). Using f 0 = v φ /s and v φ = −2∂ z Ψ 1 0 /(iω) in χ B (ζ = 0), the BC reads Figure 19 . (a) Same as in Figure 18 bulk inertial modes can also be excited (e.g. Aldridge & Toomre 1969; Zhang & Liao 2017) , 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 andC 1 ≈ 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 .

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 Greenspan 1968 ). 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 (Greenspan 1968; Aldridge & Toomre 1969; Zhang & Liao 2017) . This inertial mode excitation is not an inviscid resonance (as encountered e.g. with latitudinal librations, see in Vantieghem et al. 2015) , 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, Zhang & Liao (2017) (1 + cos 2θ − ω 12 cos θ) sin θ(i φ − θ) e λ+ζ + i 2 (1 + cos 2θ + ω 12 cos θ) sin θ (i φ + θ) e λ * − ζ e iω12t , (B 20d ) 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 Zhang & Liao (2017) , 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 Zhang & Liao (2017) , 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 Lin & Noir 2020) .

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

Integrating the coupled ordinary differential equations (3.4), we obtain Y 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 obtain Y = 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) 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

where V m is the multipolar and oscillating extension of the BC used by Suess (1971) , 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 vectorq = (q x ,q y ,q z ), and the velocity V st = V sb − n · V sb is given by the tangential components on Σ of the solid-body rotation V sb =Q x cos(ωt) x R × r. Since U 1 0 − V un is a uniform-vorticity flow, we assume thatq 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 + q z ) cos ωt   .

The last component in (B 24) allows us to reproduce the particular case m = 0 of expression (B 23b). Since non-zero Q z + q z will only be considered when m = 0, we simplify the equations by replacingQ z by Q z in expression (B 23b), with Q z = Q z + q z − Q z , and by putting Q z + q z = 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

at the boundary q 1 = Q 1 , with C l k = exp[i(kmφ + lωt)] and A (j)

Using BC (B 25), equation (B 22) gives Y , and u 1 0 reads u 1 0 · q 1 = 0, (B 27a) 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 h 2 λ * ± , ∂κ ± ∂q 2 = −sgn(γ 1 ±γ − ) γ 1 T (q1) 2 tan q 2 h 2 κ * ± , (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 = V 1 Σ − 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 haveQ z = 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 andQ z = 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 V 1 Σ (i.e. Q y = Q y = Q x = Q z = Q x = 0). Then, in the spherical geometry limit q 1 → ∞, equation (B 27) 

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

(1) ± = 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 V 1 Σ (i.e. Q x = Q y = Q y = Q x =Q x = 0). Then, equation (B 27) reduces to 

Axisymmetric inertial oscillations of a fluid in a rotating spherical container

Implicit-explicit methods for timedependent partial differential equations

Steady zonal flows in spherical shell dynamos

Observations of zonal flow created by potential vorticity mixing in a rotating fluid

Shortcut to geostrophy in wave-driven rotating turbulence: the quartetic instability

Effects of bottom topography on the spin-up in a cylinder

Shear flow instabilities in rotating systems

Steady fluid flow in a precessing spheroidal shell

Thermal instabilities in rapidly rotating systems

Mean zonal flows generated by librations of a rotating spherical cavity

Characteristics of Ekman boundary layer instabilities

Axisymmetric simulations of libration-driven fluid dynamics in a spherical shell geometry

Precessing spherical shells: flows, dissipation, dynamo and the lunar core

Elliptical instability in terrestrial planets and moons

Libration-driven multipolar instabilities

Simulations of fluid motion in spheroidal planetary cores driven by latitudinal libration

Zonal flow driven by strongly supercritical convection in rotating spherical shells

Time-dependent kinematic dynamos with stationary flows

Non-linear evolution of tidally forced inertial waves in rotating fluid bodies

Generation and maintenance of bulk turbulence by libration-driven elliptical instability

Simulation of high-Reynolds number vascular flows

Structure and dynamics of rotating turbulence: a review of recent experimental and numerical results

Tidally forced turbulence in planetary interiors

The theory of rotating fluids

On the non-linear interaction of inertial modes

On a time-dependent motion of a rotating fluid

Large-scale vortices in rapidly rotating Rayleigh-Bénard convection

Oscillatory internal shear layers in rotating and precessing flows

Conical shear-driven parametric instability of steady flow in precessing spheroids

The instability of precessing flow

On the internal shear layers spawned by the critical regions in oscillatory ekman boundary layers

Secondary instabilities in rapidly rotating fluids: inertial wave breakdown

Elliptical instability

Steady flow in a rapidly rotating sphere with weak precession

Steady flow in a rapidly rotating spheroid with weak precession: I

On magnetostrophic inertia-less waves in quasigeostrophic models of planetary cores

Flows driven by libration, precession, and tides

Dissipation at the core-mantle boundary on a small-scale topography

Experimental study of the nonlinear saturation of the elliptical instability: inertial wave turbulence versus geostrophic turbulence

Near-resonant instability of geostrophic modes: beyond Greenspan's theorem

Libration-driven flows in ellipsoidal shells

Physical constraints on the coefficients of Fourier expansions in cylindrical coordinates

Shear-driven parametric instability in a precessing sphere

Libration-driven inertial waves and mean zonal flows in spherical shells

A comparison of numerical schemes to solve the magnetic induction eigenvalue problem in a spherical geometry

Fluid instabilities in precessing spheroidal cavities

Precession of the Earth as the cause of geomagnetism: Experiments lend support to the proposal that precessional torques drive the Earth's dynamo

Full sphere hydrodynamic and dynamo benchmarks

Rotating double-diffusive convection in stably stratified planetary cores

Experimental determination of zonal winds driven by tides

On a small-scale roughness of the core-mantle boundary

Rossby wave packet interactions

Hysteresis and instabilities in a spheroid in precession near the resonance with the tilt-over mode

Experimental evidence of inertial waves in a precessing spheroidal cavity

Precession-driven flows in non-axisymmetric ellipsoids

Experimental study of libration-driven zonal flows in non-axisymmetric containers

An experimental and numerical study of librationally driven flow in planetary cores and subsurface oceans

Numerical study of the motions within a slowly precessing sphere at low Ekman number

Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum

On the theory of core-mantle coupling

Mean zonal flow generated by azimuthal harmonic forcing in a rotating cylinder

Spontaneous generation of inertial waves from boundary turbulence in a librating sphere

Fluid flows in a librating cylinder

Experimental and numerical study of mean zonal flows generated by librations of a rotating spherical cavity

Tide-driven shear instability in planetary liquid cores

Libration-induced mean flow in a spherical shell

Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations

Quasigeostrophic model of the instabilities of the Stewartson layer in flat and depth-varying containers

Turbulent geodynamo simulations: a leap towards Earth's core

Numerical study of a rotating fluid in a spheroidal container

Transfer of energy to two-dimensional large scales in forced, rotating three-dimensional turbulence

Friction law and turbulent properties in a laboratory Ekman boundary layer

On almost rigid rotations. Part 2

Acoustic spectra of a gas-filled rotating spheroid

Viscous flow in a deformable rotating container

Latitudinal libration driven flows in triaxial ellipsoids

Inviscid instabilities in rotating ellipsoids on eccentric Kepler orbits

Fossil field decay due to nonlinear tides in massive binaries

Compressible fluid modes in rigid ellipsoids: towards modal acoustic velocimetry

Cylindrical tank of fluid oscillating about a state of steady rotation

Viscous corrections to stewartson's stability criterion(viscous corrections to stewartson theory on stability of spinning top containing liquid)

Lunar rotational dissipation in solid body and molten core

Asymptotic theory of resonant flow in a spheroidal cavity driven by latitudinal libration

Theory and modeling of rotating fluids: convection, inertial waves and precession

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 towith A 0 = ∂ r P m 1 | r=0 . 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)