Precession-driven flows in stress-free ellipsoids

Abstract

Motivated by modelling rotating turbulence in planetary fluid layers, we investigate precession-driven flows in ellipsoids subject to stress-free boundary conditions (SF-BC). The SF-BC could indeed unlock numerical constraints associated with the no-slip boundary conditions (NS-BC), but are also relevant for some astrophysical applications. Although SF-BC have been employed in the pioneering work of Lorenzani & Tilgner (J. Fluid Mech., vol. 492, 2003, pp. 363–379), they have scarcely been used due to the discovery of some specific mathematical issues associated with angular momentum conservation. We revisit the problem using asymptotic analysis in the low-viscosity regime, which is validated with numerical simulations. First, we extend the reduced model of uniform-vorticity flows in ellipsoids to account for SF-BC. We show that the long-term evolution of angular momentum is affected by viscosity in triaxial geometries, but also in axisymmetric ellipsoids when the mean rotation axis of the fluid is not the symmetry axis. In a regime relevant to planets, we analytically obtain the primary forced flow in triaxial geometries, which exhibits a second inviscid resonance. Then, we investigate the bulk instabilities existing in precessing ellipsoids. We show that using SF-BC would be useful to explore the non-viscous instabilities (e.g. Kerswell, Geophys. Astrophys. Fluid Dyn., vol. 72, 1993, pp. 107–144), which are presumably relevant for planetary applications but are often hampered in experiments or simulations with NS-BC.

1. Introduction

Motivated by numerous natural applications (e.g. Le Bars, Cébron & Le Gal 2015), we aim to explore the long-term dynamics of rapidly rotating fluids enclosed in ellipsoids subject to (harmonic) mechanical forcings. Global rotation is indeed ubiquitous in many planetary fluid layers or stars, which are usually ellipsoidal at the leading order (e.g. due to the combined action of centrifugal effects and gravitational interactions with nearby orbital partners, see Chandrasekhar 1969). In particular, mechanically driven flows in ellipsoids (e.g. flows driven by precession or tides) have received much attention in the fluid community. Mechanical forcings can indeed sustain bulk instabilities (e.g. Kerswell 1993,  2002), turbulence (e.g. Grannan et al. 2017; Le Reun, Favier & Le Bars 2019) and possibly dynamo magnetic fields (e.g. Reddy, Favier & Le Bars 2018; Vidal et al. 2018). These works have also renewed interest in a key fundamental question in the theory of rotating fluids, which is the generation of two-dimensional geostrophic motions (Greenspan 1969). However, this problem has only received scant attention in global geometries exhibiting the so-called topographic beta effect (which strongly modifies the geostrophic flows, e.g. Greenspan 1968). Exploring rotating turbulence thus deserves further work using global models.

The incompressible Navier–Stokes equation is commonly adopted to explore the turbulence driven by mechanical forcings, together with the no-slip boundary conditions (NS-BC). The latter are appropriate to model the flow dynamics in the presence of a rigid boundary (e.g. the solid interface between a liquid core and a solid overlying mantle in planetary interiors). However, the range of parameters that is accessible to global simulations with NS-BC is severely limited, in particular for the Ekman number $E$ (which crucially controls the dynamics of rapidly rotating flows). Typical values in natural systems are $E \leq {O}(10^{-12})$, whereas direct numerical simulations (DNS) and laboratory experiments of mechanically driven rotating turbulence can only reach much larger values $E \gtrsim 10^{-6}$ (e.g. Grannan et al. 2017; Le Reun et al. 2019). As a consequence, the Ekman boundary layer is often a prominent feature in the models (whereas the smallness of $E$ in planetary systems suggests that viscosity should rather play a minor dynamical role), and its resolution requires considerable computational resources when $E$ is lowered. Moreover, the overestimated viscous torque at the boundary can also largely inhibit the fluid response to mechanical forcings (which is primarily driven by the shape deformation of the fluid boundary, combined with non-stationary effects due to the possibly oscillatory angular velocity of the container). Therefore, different modelling approaches are worth considering to simulate such flows at more realistic parameters for planetary applications.

One natural way to avoid the physical and computational disadvantages of NS-BC is to employ stress-free boundary conditions (SF-BC). A thin outer Ekman boundary layer is still present for stress-free boundaries (e.g. Livermore, Bailey & Hollerbach 2016), but its dynamical role is expected to be less important because the boundary-layer flow is much weaker in amplitude than the bulk flow (e.g. Rieutord 1992). Moreover, SF-BC are also commonly employed in astrophysical modelling since they are often believed to yield similar results to those obtained with a realistic free surface (Barker 2016a). However, SF-BC have scarcely been used in spheres and ellipsoids because of mathematical difficulties. The most serious one is related to angular momentum conservation. Angular momentum can indeed be arbitrary in axisymmetric geometries, leading to spurious solutions on long time scales (e.g. Jones et al. 2011; Guermond et al. 2013). The usefulness of SF-BC for simulating rotating flows in ellipsoids has thus been questioned, but we believe that this mathematical set-up deserves further analysis.

In this paper, we thus revisit the influence of SF-BC for rotating ellipsoids using asymptotic analysis when $E \ll 1$ and targeted numerical simulations. The paper is organised as follows. The model is presented in § 2 and applied to precessing ellipsoids in § 3. The results are discussed in § 4, and we end the paper in § 5.

2. Mathematical modelling

2.1. Fluid dynamic equations

We consider a fluid-filled ellipsoid of uniform density and volume $V$, which is assumed to co-rotate with the surrounding mantle at the angular velocity $\boldsymbol {\varOmega }_c (t) = \varOmega _0 [ \boldsymbol {\varOmega } + \boldsymbol {\delta }(t) ]$ with respect to the inertial frame ($\boldsymbol {\delta }(t)$ being the time-dependent departure from the steady global rotation $\boldsymbol {\varOmega }$ along the unit vector $\boldsymbol {1}_\varOmega = \boldsymbol {\varOmega }/|\boldsymbol {\varOmega }|$). To have a tractable mathematical problem, we seek mechanically driven flows in the mantle reference frame in which the ellipsoidal boundary $S$ is steady and $\boldsymbol {\delta }(t) \neq \boldsymbol {0}$. This set-up allows us to model flows driven by precession or librations, which have already received consideration using NS-BC (e.g. Zhang, Chan & Liao 2012,  2014; Noir & Cébron 2013; Vantieghem, Cébron & Noir 2015). We non-dimensionalise the problem using $\varOmega _0^{-1}$ as the time scale, and a typical length $R$ as the length scale (which is here arbitrary). Considering a Newtonian fluid of uniform kinematic viscosity $\nu$, the dimensionless equations for the velocity $\boldsymbol {v}$ are

(2.1a)\begin{gather} \partial_t \boldsymbol{v} + ( \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} ) \boldsymbol{v} + 2 \boldsymbol{\varOmega}_c \times \boldsymbol{v} ={-} \boldsymbol{\nabla} p + 2 E \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\epsilon}(\boldsymbol{v}) + \boldsymbol{r} \times \mathrm{d}_t \boldsymbol{\delta}, \end{gather}

(2.1b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} = 0, \end{gather}

where $\boldsymbol {r}$ is the position vector, $\boldsymbol {\epsilon } (\boldsymbol {v}) =(1/2) [ \boldsymbol {\nabla } \boldsymbol {v} + (\boldsymbol {\nabla } \boldsymbol {v})^\top ]$ is the strain-rate tensor and $E = \nu /(\varOmega _0 R^2)$ is the Ekman number. The ellipsoidal geometry, which is assumed to be steady in the mantle frame, is given by the dimensionless equation

(2.2)\begin{equation} (x/a)^2 + (y/b)^2 + (z/c)^2 = 1, \end{equation}

where $[a,b,c]$ are the (dimensionless) ellipsoidal semi-axes and $[x,y,z]$ are the Cartesian coordinates. In the following, axisymmetric geometries refer to ellipsoids with a revolution symmetry axis (i.e. when either $a=b, b=c$ or $a=c$). Finally, spheroids will refer to the particular axisymmetric geometries for which the revolution symmetry axis is aligned with the rotation axis (with $a=b$ and $\boldsymbol {\varOmega } \propto \boldsymbol {1}_z$ in this study). We aim to consider the SF-BC given in the mantle frame by

(2.3a,b)\begin{equation} \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{1}_n |_S = 0, \quad [ \boldsymbol{\epsilon} (\boldsymbol{v}) \boldsymbol{\cdot} \boldsymbol{1}_n] \times \boldsymbol{1}_n |_S = \boldsymbol{0}, \end{equation}

where $\boldsymbol {1}_n$ is the outward normal unit vector at the boundary, instead of the NS-BC

(2.4a,b)\begin{equation} \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{1}_n |_S = 0, \quad \boldsymbol{v} \times \boldsymbol{1}_n |_S = \boldsymbol{0}. \end{equation}

It is obvious from SF-BC (2.3) and NS-BC (2.4) that the tangential velocity at the boundary will differ between the two cases (since the flow is allowed to freely slip on the boundary with the SF-BC). One may thus wonder in which circumstances the above conditions will lead to similar flows in the bulk (i.e. far from the boundary region).

A necessary condition is that the mechanical forcings can sustain flows against viscous dissipation for the two BC in the mantle frame. This is evidenced by the conservation equation for the volume-averaged kinetic energy $E_k$. In a frame where the fluid boundary is steady, it is given by (e.g. equation (5) in Wu & Roberts 2009)

(2.5)\begin{equation} \mathrm{d}_t E_k = \int_V \boldsymbol{v} \boldsymbol{\cdot} [ \boldsymbol{r} \times \mathrm{d}_t \boldsymbol{\delta} ] \, \mathrm{d} V + 2 E \left( \int_S \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\mathcal{T}} \, \mathrm{d}S - \mathcal{D}_\nu \right),\end{equation}

where $\boldsymbol {\mathcal {T}} = \boldsymbol {\epsilon }(\boldsymbol {v}) \boldsymbol {\cdot } \boldsymbol {1}_n$ is the surface traction and $\mathcal {D}_v \geq 0$ is a volume-averaged viscous dissipation (for both the NS-BC and SF-BC). For a velocity satisfying the no-penetration condition such that $\boldsymbol {v} = (\boldsymbol {v} \boldsymbol {\cdot } \boldsymbol {1}_n) \boldsymbol {1}_n - \boldsymbol {1}_n \times (\boldsymbol {1}_n \times \boldsymbol {v}) = - \boldsymbol {1}_n \times (\boldsymbol {1}_n \times \boldsymbol {v})$, the surface integral can actually be written as

(2.6)\begin{equation} \int_S \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\mathcal{T}} \, \mathrm{d}S ={-} \int_S \boldsymbol{\mathcal{T}} \boldsymbol{\cdot} [ \boldsymbol{1}_n \times (\boldsymbol{1}_n \times \boldsymbol{v}) ] \, \mathrm{d} S ={-} \int_S [ \boldsymbol{\mathcal{T}} \times \boldsymbol{1}_n ] \boldsymbol{\cdot} [\boldsymbol{v} \times \boldsymbol{1}_n ] \, \mathrm{d}S, \end{equation}

where we have used a property of the scalar triple product to obtain the last expression. Thus, the above surface integral exactly vanishes for both SF-BC (2.3) and NS-BC (2.4) in the mantle frame. Then, equation (2.5) shows that we can have $\mathrm {d}_t E_k \geq 0$ for both SF-BC and NS-BC if the mechanical forcings are oscillatory in the mantle frame (i.e. when $\mathrm {d}_t \boldsymbol {\delta } \neq \boldsymbol {0}$). Harmonic mechanical forcings, such as precession or librations, can thus sustain flows against viscous dissipation in the mantle frame (even with the SF-BC). Note that a very different conclusion is obtained for steady forcings, such as precession viewed in the frame of precession for spheroidal geometries (Lorenzani & Tilgner 2003; Wu & Roberts 2009). We indeed have $\mathrm {d}_t E_k < 0$ at every time for the SF-BC in the precession frame, whereas precession could sustain non-vanishing flows against viscous dissipation for the NS-BC (since $\boldsymbol {v} \times \boldsymbol {1}_n |_S \neq \boldsymbol {0}$ for a no-slip boundary in the precession frame, e.g. Cébron et al. 2019). In the following, we will only investigate the dynamics driven by oscillatory forcings in the mantle frame with SF-BC.

2.2. Angular momentum

The angular momentum $\boldsymbol {L} = \int _V \boldsymbol {r} \times \boldsymbol {v} \, \mathrm {d}V$ of the flow plays a central dynamical role for mechanically driven flows in ellipsoids. Actually, the Cartesian components of the angular momentum $\boldsymbol {L} = (L_x, L_y, L_z)^\top$ are exactly given for incompressible flows by

(2.7a)\begin{gather} L_x = \int_V ( y v_z - z v_y ) \, \mathrm{d} V = \int_V ( \boldsymbol{1}_x \times \boldsymbol{r} + \boldsymbol{\nabla} \varPsi_x ) \boldsymbol{\cdot} \boldsymbol{v} \, \mathrm{d}V, \end{gather}

(2.7b)\begin{gather}L_y = \int_V ( z v_x - x v_z ) \, \mathrm{d} V = \int_V ( \boldsymbol{1}_y \times \boldsymbol{r} + \boldsymbol{\nabla} \varPsi_y ) \boldsymbol{\cdot} \boldsymbol{v} \, \mathrm{d}V, \end{gather}

(2.7c)\begin{gather}L_x = \int_V ( x v_y - y v_x ) \, \mathrm{d} V = \int_V ( \boldsymbol{1}_z \times \boldsymbol{r} + \boldsymbol{\nabla} \varPsi_z ) \boldsymbol{\cdot} \boldsymbol{v} \, \mathrm{d}V, \end{gather}

where $[\varPsi _x,\varPsi _y,\varPsi _z]$ are arbitrary scalar potentials if $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v} = 0$ and if the flow obeys the no-penetration BC in rigid ellipsoids. The scalar potentials are thus often discarded to simply express the angular momentum as projections onto the solid-body rotations $\boldsymbol {1}_i \times \boldsymbol {r}$ (e.g. Guermond et al. 2013). Yet, the solid-body rotations are not admissible flow solutions in non-spherical geometries (even without viscosity), since they do not satisfy the no-penetration condition.

A more appropriate definition of the angular momentum for incompressible flows is thus given in ellipsoids by

(2.8)\begin{equation} \boldsymbol{L} \boldsymbol{\cdot} \boldsymbol{1}_i = \int_V \boldsymbol{e}_i \boldsymbol{\cdot} \boldsymbol{v} \, \mathrm{d}V,\end{equation}

where $\{\boldsymbol {e}_i \}_{i\in \{x,y,z\}}$ is the set of uniform-vorticity (flow) elements defined by

(2.9a–c)\begin{equation} \boldsymbol{e}_i = \boldsymbol{1}_i \times \boldsymbol{r} + \boldsymbol{\nabla} \varPsi_i, \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{e}_i = 0, \quad \boldsymbol{e}_i \boldsymbol{\cdot} \boldsymbol{1}_n |_S = 0. \end{equation}

The scalar functions $\varPsi _i$ allow the elements $\boldsymbol {e}_i$ to satisfy the no-penetration condition. In ellipsoidal geometries, they are explicitly given by (e.g. Noir & Cébron 2013)

(2.10a–c)\begin{equation} \varPsi_x = \frac{c^2-b^2}{b^2+c^2} yz, \quad \varPsi_y = \frac{a^2-c^2}{a^2+c^2} xz, \quad \varPsi_z = \frac{b^2-a^2}{a^2 + b^2} xy. \end{equation}

It is worth noting that definiton (2.8) is purely kinematic. It thus remains valid in the presence of additional effects, for instance without global rotation or with magnetic effects (e.g. Gerick et al. 2020). Moreover, this definition can also be generalised for compressible flows under the anelastic approximation (see Appendix A). Consequently, we can always rigorously expand incompressible velocity fields in ellipsoids as

(2.11a,b)\begin{equation} \boldsymbol{v}(\boldsymbol{r},t) = \boldsymbol{U}(\boldsymbol{r},t) + \boldsymbol{v}_f (\boldsymbol{r},t), \quad \int_V \boldsymbol{r} \times \boldsymbol{v}_f \, \mathrm{d} V = \boldsymbol{0},\end{equation}

where the uniform-vorticity flow $\boldsymbol {U}$ carrying the angular momentum is given by

(2.12a,b)\begin{equation} \boldsymbol{U}(\boldsymbol{r},t) = \omega_x (t) \boldsymbol{e}_x (\boldsymbol{r}) + \omega_y (t) \boldsymbol{e}_y (\boldsymbol{r}) + \omega_z (t) \boldsymbol{e}_z (\boldsymbol{r}), \quad \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{1}_n |_S = 0, \end{equation}

and with the effective rotation vector of the fluid $\boldsymbol {\omega } (t) =(\omega _x (t), \omega _y (t), \omega _z (t))^\top$. The velocity $\boldsymbol {v}_f$, which does not carry angular momentum by definition since $\int _V \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {v}_f \, \mathrm {d} V=0$, contains bulks flows of higher spatial complexity (e.g. flow instabilities or turbulence) and also viscous structures (e.g. the Ekman boundary layer, Rieutord 1992). The Cartesian components of $\boldsymbol {L}$ are then exactly given by

(2.13a,b)\begin{equation} \boldsymbol{L} =\mathcal{L}^{{-}1} \boldsymbol{\omega}, \quad \mathcal{L}^{{-}1} = \frac{16{\rm \pi}}{15} abc\ \mathrm{diag} \left[ \frac{b^2c^2}{b^2+c^2}, \frac{a^2c^2}{a^2+c^2}, \frac{a^2b^2}{a^2+b^2} \right]. \end{equation}

Finally, the time evolution of the angular momentum (or equivalently that of $\boldsymbol {\omega }$) is affected by viscosity through the action of the viscous torque $\boldsymbol {\varGamma }_\nu$ on long time scales. We have for example $\boldsymbol {\varGamma }_\nu = \boldsymbol {0}$ in spheres, such that angular momentum has to be conserved for uniformly rotating fluids in the inertial frame (e.g. Jones et al. 2011). To clarify the dynamical role of SF-BC in ellipsoids, it is worth computing the viscous torque.

2.3. Viscous torque in stress-free ellipsoids

Because of definition (2.8), the Cartesian components of the viscous torque $\boldsymbol {\varGamma }_\nu = (\boldsymbol {\varGamma }_{\nu } \boldsymbol {\cdot } \boldsymbol {1}_x, \boldsymbol {\varGamma }_{\nu } \boldsymbol {\cdot } \boldsymbol {1}_y, \boldsymbol {\varGamma }_{\nu } \boldsymbol {\cdot } \boldsymbol {1}_z)^\top$ are exactly given for SF-BC (2.3) by

(2.14)\begin{equation} \boldsymbol{\varGamma}_{\nu} \boldsymbol{\cdot} \boldsymbol{1}_i = 2 E \int_V \boldsymbol{e}_i \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\epsilon}(\boldsymbol{v}) \, \mathrm{d} V ={-} 2E \int_V \boldsymbol{\epsilon}(\boldsymbol{e}_i) : \boldsymbol{\epsilon}(\boldsymbol{v}) \, \mathrm{d} V,\end{equation}

where we have used integration by parts and the decomposition $\boldsymbol {e}_i = (\boldsymbol {1}_n \boldsymbol {\cdot } \boldsymbol {e}_i) \boldsymbol {1}_n - \boldsymbol {1}_n \times (\boldsymbol {1}_n \times \boldsymbol {e}_i) = - \boldsymbol {1}_n \times (\boldsymbol {1}_n \times \boldsymbol {e}_i)$ to cancel out the surface integral for SF-BC (e.g. see the proof of proposition 2.1 in Guermond et al. 2013). We recover from the formula that $\boldsymbol {\varGamma }_\nu = \boldsymbol {0}$ in spheres since $\boldsymbol {\epsilon }(\boldsymbol {e}_i)$ exactly vanishes when $\boldsymbol {e}_i$ is a solid-body rotation, but we also obtain that $\boldsymbol {\varGamma }_\nu \neq \boldsymbol {0}$ in triaxial geometries (because $\boldsymbol {\epsilon }(\boldsymbol {e}_i) \neq 0$ when $a\neq b \neq c$). Moreover, it shows that $\boldsymbol {\varGamma }_{\nu } \boldsymbol {\cdot } \boldsymbol {1}_i =0$ when the Cartesian vector $\boldsymbol {1}_i$ is an axis of revolution of the geometry (irrespective of the fluid global rotation, as $\boldsymbol {e}_i$ is then a solid-body rotation).

We can now inspect the long-term evolution of angular momentum since pathological behaviours have been reported in some axisymmetric configurations (Guermond et al. 2013). To illustrate this behaviour, we expand the angular momentum as $\boldsymbol {L}=\boldsymbol {L}_{0} + \boldsymbol {L}_{1}$, where $\boldsymbol {L}_{0}$ is the angular momentum of a dynamical solution of the problem and $\boldsymbol {L}_{1}$ is a modification of $\boldsymbol {L}_{0}$ associated with an additional uniform-vorticity flow. The time evolution of $\boldsymbol {L}_{1}$ is then given in the rotating frame by (e.g. Roberts & Aurnou 2012)

(2.15)\begin{equation} \mathrm{d}_t \boldsymbol{L}_{1} + \boldsymbol{\varOmega}_c \times \boldsymbol{L}_{1} = \boldsymbol{\varGamma}_{p,{1}} + \boldsymbol{\varGamma}_{\nu,{1}},\end{equation}

where $\boldsymbol {\varGamma }_{p,{1}} = \int _S p_{1} \boldsymbol {1}_n \times \boldsymbol {r} \, \mathrm {d} S$ is the pressure torque and $\boldsymbol {\varGamma }_{\nu,{1}}$ is the viscous torque. Since the viscous and pressure torques are non-zero when $a \neq b \neq c$, equation (2.15) shows that the angular momentum is affected by viscosity in triaxial ellipsoids. The situation is possibly different in axisymmetric geometries. If the fluid is not globally rotating (i.e. when $\boldsymbol {\varOmega } = \boldsymbol {0}$), then the component $\boldsymbol {L}_{1} \boldsymbol {\cdot } \boldsymbol {1}_i$ carried by the uniform-vorticity element $\boldsymbol {e}_i$ is arbitrary when $\boldsymbol {1}_i$ is a revolution symmetry axis (since $\boldsymbol {\varGamma }_{p,{1}} \boldsymbol {\cdot } \boldsymbol {1}_i = \boldsymbol {\varGamma }_{\nu,{1}} \boldsymbol {\cdot } \boldsymbol {1}_i = 0$). Similarly, if the fluid is globally rotating along the revolution symmetry axis $\boldsymbol {1}_i$, then the perturbation angular momentum $\boldsymbol {L}_{1} \propto \boldsymbol {1}_i$ is arbitrary (it will depend on the initial conditions, e.g. as shown in Guermond et al. 2013).

The two situations are illustrated numerically in figure 1 for a spheroid $a=b$ subject to the precession forcing (see its definition below in § 3). We have performed DNS using the standard finite-element method as implemented in the commercial software comsol. The latter has already been employed to simulate precession-driven flows in ellipsoids with NS-BC (e.g. Noir & Cébron 2013) and can also account for SF-BC (e.g. for tidal flows in Cébron et al. 2013). The geometry is modelled by an unstructured mesh with tetrahedral elements in the bulk, surrounded by a boundary-layer mesh (made of prism elements) to ensure the convergence of the thin Ekman layer. We have employed Lagrange elements $P2$–$P3$ (i.e. quadratic for the pressure field and cubic for the velocity field). The total number of degrees of freedom ranges between $3 \times 10^5$ and $5 \times 10^5$, such that every targeted simulation took a few days to run in parallel on a cluster (to investigate the long-term evolution of $\boldsymbol {L}$). We observe that the axial angular momentum $L_z$ does not converge in time for the considered stress-free spheroid (it is still growing or decaying even after several viscous time scales) if either the fluid is non-rotating in average as in panel (a) or $\boldsymbol {\varOmega } \propto \boldsymbol {1}_z$ as in panel (b). However, a definitive conclusion should not be drawn for every axisymmetric geometry. The situation is indeed different if the global rotation is not aligned with the revolution axis, since the three components of the angular momentum should be strongly coupled in (2.15) for such configurations (even if $\boldsymbol {\varGamma }_\nu \boldsymbol {\cdot } \boldsymbol {1}_i = 0$, see § 3).

Figure 1. Non-convergence of the angular momentum $L_z$ in DNS after several viscous time units. Precession forcing given by definition (3.1) with $P_x = 10^{-2}$ in stress-free spheroids ($a=b=1, c=0.95$). At $t=0, [\omega _x,\omega _y]$ are chosen to match asymptotic solution (3.7). (a) DNS at $Po=-1$ for the two values of the Ekman number $E=5 \times 10^{-3}$ (e.g. as considered in Wu & Roberts 2009) and $E=5 \times 10^{-4}$. At $t=0, \omega _z \approx 0$ for the two simulations. (b) DNS at $Po=-1.8$ and $E=5\times 10^{-4}$ for $\omega _z \approx 0$ (top panel) and $\omega _z = 0.1$ (bottom panel) at $t=0$. (a) Non-rotating $\boldsymbol {\varOmega } \simeq \boldsymbol {0}$ and (b) Rotating $\boldsymbol {\varOmega } \propto \boldsymbol {1}_z$.

3. Application to precession-driven flows

We consider precession-driven flows in ellipsoids, which have only received scant attention with SF-BC (Lorenzani & Tilgner 2003; Wu & Roberts 2009; Guermond et al. 2013). We work in the mantle frame rotating with respect to the inertial frame at the dimensionless angular velocity (e.g. Noir & Cébron 2013)

(3.1a,b)\begin{equation} \boldsymbol{\varOmega}_c (t) = \underbrace{(1 + P_z) \boldsymbol{1}_z}_{\boldsymbol{\varOmega}} + \boldsymbol{\delta} (t), \quad \boldsymbol{\delta}(t) = P_x [ \cos(t) \boldsymbol{1}_x - \sin(t) \boldsymbol{1}_y ], \end{equation}

with $P_x = Po\, \sin(\alpha)$ and $P_z = Po \, \cos(\alpha)$, where $Po = \varOmega_p/\varOmega_0$ is the Poincaré number ($\varOmega_p$ being the angular velocity of precession and $\varOmega_0$ that of the mantle) and $\alpha$ is the angle of precession measured from $\boldsymbol{1}_z$. Because the Poincaré force $\boldsymbol {r}\times \mathrm {d}_t \boldsymbol {\delta }$ is linear in Cartesian coordinates, the primary response of the fluid is a laminar uniform-vorticity flow (e.g. Noir & Cébron 2013; Kida 2020), on top of which secondary flows and turbulence can develop. For analytical progress, we expand the velocity field as $\boldsymbol {v}=\boldsymbol {v}_{0}+\boldsymbol {v}_{1}$, where $\boldsymbol {v}_{0}$ is the primary forced flow (which is mainly of uniform vorticity) and $\boldsymbol {v}_{1}$ represents small-amplitude additional flows such that $|\boldsymbol {v}_{1}| \ll |\boldsymbol {v}_{0}|$. We first seek analytical solutions of the primary flow in § 3.1, which are compared with DNS in § 3.2. Then, we explore the flow instabilities $\boldsymbol {v}_{1}$ growing upon the forced flow in § 3.3.

3.1. Laminar forced flows

The forced laminar flows, which have been explored for a long time after the seminal work of Poincaré (1910), can be obtained using boundary-layer theory (BLT) in the low-viscosity regime $E \ll 1$ for SF-BC. To do so, we seek $\boldsymbol {v}_{0}$ as

(3.2)\begin{equation} \boldsymbol{v}_{0}(\boldsymbol{r},t) \simeq \underbrace{\omega_x (t) \boldsymbol{e}_x + \omega_y (t) \boldsymbol{e}_y + \omega_z (t) \boldsymbol{e}_z}_{\boldsymbol{U}(\boldsymbol{r},t)} + E^{1/2} \tilde{\boldsymbol{U}}(\boldsymbol{r},t), \end{equation}

where $\boldsymbol {U}(\boldsymbol {r},t)$ is a forced uniform-vorticity flow carrying angular momentum, and $\tilde {\boldsymbol {U}}(\boldsymbol {r},t)$ is the viscous flow within the boundary layer at the leading order in $E^{1/2}$ (e.g. Rieutord 1992). A direct consequence of asymptotic expansion (3.2) is that the bulk flow for SF-BC can be determined without explicitly solving for the boundary-layer flow (since the latter has an amplitude that is $E^{1/2}$ smaller than the bulk flow amplitude). The exact viscous torque given by formula (2.14) can then be approximated as

(3.3)\begin{equation} \boldsymbol{\varGamma}_\nu \simeq{-}\frac{16 {\rm \pi}}{3}abc \,E\, \mathrm{diag} \left[ \frac{(b^2-c^2)^2}{(b^2+c^2)^2}, \frac{(a^2-c^2)^2}{(a^2+c^2)^2}, \frac{(a^2-b^2)^2}{(a^2+b^2)^2} \right] \boldsymbol{\omega}.\end{equation}

The viscous flow $E^{1/2} \tilde {\boldsymbol {U}}$ in expansion (3.2) has a contribution of amplitude ${O}(E^{3/2})$ to the viscous torque (since $|\boldsymbol {\epsilon } (\tilde {\boldsymbol {U}})| = {O}(E^{-1/2})$ and the volume scales as ${O}(E^{1/2})$ within the Ekman layer), which can be neglected compared with expression (3.3) in the asymptotic regime $E \ll 1$. We recover from formula (3.3) that $\boldsymbol {\varGamma }_{\nu } \boldsymbol {\cdot } \boldsymbol {1}_i =0$ when the Cartesian vector $\boldsymbol {1}_i$ is a revolution symmetry axis, but also that the three components of the viscous torque are non-zero when $a \neq b \neq c$. Then, the momentum equation reduces to

(3.4)\begin{equation} \mathrm{d}_t \boldsymbol{\omega} - [( \boldsymbol{\omega} + \boldsymbol{\varOmega}_c ) \boldsymbol{\cdot} \boldsymbol{\nabla} ] \boldsymbol{U} ={-} \mathrm{d}_t \boldsymbol{\varOmega}_c + \mathcal{L} \boldsymbol{\varGamma}_\nu, \end{equation}

where $\boldsymbol {\varGamma }_\nu$ is the viscous torque given by formula (3.3) and $\mathcal {L}$ is the matrix given by the inverse of expression (2.13b). The approximated viscous term is thus

(3.5)\begin{equation} \mathcal{L} \boldsymbol{\varGamma}_\nu ={-}5 E \,\mathrm{diag} \left[ \frac{(b/c-c/b)^2}{b^2+c^2}, \frac{(a/c-c/a)^2}{a^2+c^2}, \frac{(a/b-b/a)^2}{a^2+b^2} \right] \boldsymbol{\omega}.\end{equation}

Equations (3.4) and (3.5) extend the asymptotic viscous model of Noir & Cébron (2013) to stress-free ellipsoids, but we remind the reader that this stress-free model is not valid in spheres (since the angular momentum would be arbitrary in spheres because of $\boldsymbol {\varGamma }_\nu = \boldsymbol {0}$). The close similarity between the no-slip and stress-free cases, for which only the expression of the viscous term in (3.4) differs, suggests that the same interior solution should be approached when $E \to 0$ in no-slip and stress-free ellipsoids.

Precession is often characterised by $|P_x| \ll 1$ in planetary liquid cores (e.g. Noir & Cébron 2013). Hence, we seek asymptotic solutions of (3.4) in powers of $P_x$ as

(3.6)\begin{equation} \boldsymbol{\omega} (t) = \boldsymbol{\omega}^{(0)} (t) + P_x \boldsymbol{\omega}^{(1)} (t) + P_x^2 \boldsymbol{\omega}^{(2)} (t) + \cdots . \end{equation}

Since the mean rotation axis is $\boldsymbol {\varOmega } \propto \boldsymbol {1}_z$ when $|P_x| \ll 1$, we assume that $a \neq b$ (to avoid the pathological situations outlined in § 2 for the angular momentum conservation). The zeroth-order solution $\boldsymbol {\omega }^{(0)} (t)$ corresponds to a decaying transient when $t\to \infty$ (because of viscosity). We thus discard $\boldsymbol {\omega }^{(0)} (t)$ in the following and solve the first-order problem in $P_x$. In the regime of vanishing viscosity $E \to 0$, we obtain the first-order solution

(3.7a,b)\begin{equation} \omega_x^{(1)} (t) \simeq{-} \frac{1+[1+ P_z] A_1}{1-[1 + P_z]^2 \lambda_{so}^2} \cos(t), \quad \omega_y^{(1)} (t) \simeq \frac{1+[1+P_z] B_2}{1-[1 + P_z]^2 \lambda_{so}^2} \sin(t), \end{equation}

and $\omega _z^{(1)} (t) \to 0$, with $A_1 = 2 a^2/(a^2+c^2), B_2 = {2 b^2}/(b^2+c^2)$ and $\lambda _{so}=\sqrt {A_1 B_2}$. We have finally to compute the second-order solution $\boldsymbol {\omega }^{(2)}$, accounting for weakly nonlinear interactions in the viscous interior, to estimate the axial angular velocity (since it is undefined at the first order). An analytical solution can be obtained when $E \neq 0$, showing that $\boldsymbol {\omega }^{(2)} = \omega _z^{(2)} \boldsymbol {1}_z$, but the general expression of $\omega _z^{(2)}$ is too lengthy to be given here. In the regime of vanishing viscosity $E \to 0$, it simplifies into

(3.8a)\begin{equation} \omega_z^{(2)} (t) = \frac{c^2}{4 \mathcal{D}_2^2} [\bar{\omega}_z^{(2)} + \delta \omega_z^{(2)} \cos(2t)] \end{equation}

with the denominator $\mathcal {D}_2 = a^2b^2 \tilde {P}_z (P_z+1/2 ) -c^2 (a^2+b^2+c^2 )/4$ and $\tilde {P}_z = P_z + 3/2$, where the amplitude of the mean geostrophic flow is given by

(3.8b)\begin{align} \bar{\omega}_z^{(2)} ={-} \left(\frac{c^2}{2} + a^2 \tilde{P}_z \right) \left(\frac{c^2}{2} + b^2 \tilde{P}_z \right) \frac{a^2+b^2}{(a/b-b/a)^2} \left[\left(\frac{a}{c}-\frac{c}{a} \right)^2 + \left(\frac{b}{c}-\frac{c}{b} \right)^2 \right] \end{align}

and that of the oscillatory component by $\delta \omega _z^{(2)} = (P_z+1)(a^2-b^2) ( a^2b^2 {\tilde {P}_z}^2 - c^4/4 )$. It is worth noting that the mean geostrophic flow $\bar {\omega }_z^{(2)}$ has an amplitude that is independent of $E$ in the vanishing regime $E \to 0$, which is somehow similar to the mean geostrophic flows driven by nonlinear boundary-layer interactions for NS-BC (e.g. Cébron et al. 2021).

A striking property of the asymptotic solution is that it exhibits two inviscid direct resonances, which occur when the common denominator in expressions (3.7a,b) vanishes at the two resonant values $Po^\pm$ given by $\lambda _{so} [ 1 + Po^\pm \cos (\alpha ) ] = \pm 1$. The resonance associated with $Po^+$ actually corresponds to the inviscid resonance initially predicted by Poincaré (1910), which has been observed for no-slip boundaries (e.g. Vormann & Hansen 2018; Nobili et al. 2021; Burmann & Noir 2022). However, the second resonance at $Po^-$ is new, although precession-driven flows have been explored for more than a century in triaxial ellipsoids (e.g. Poincaré 1910; Noir & Cébron 2013).

3.2. Numerical simulations

We have checked that the analytic expressions are in excellent agreement with the numerical integration of the exact uniform-vorticity model (3.4) when $E \to 0$ (not shown). Yet, it remains to confirm the validity of the asymptotic solutions against DNS with SF-BC. We first show in figure 2 the time evolution of the rotation vector $\boldsymbol {\omega } (t)$ in the DNS (performed with comsol, as explained in § 2). We illustrate the DNS at $P_x=10^{-2}$ with $Po=-1.8$ and $E=5 \times 10^{-4}$, in the particular axisymmetric geometry $a=1.5$ and $b=c=1$ (other parameters yield similar results, not shown). The fluid angular velocity $\boldsymbol {\omega }$ has been computed in the DNS using either the volume-averaged vorticity or formula (2.13a) after having computed the angular momentum. Both methods are found to be in excellent quantitative agreement for the SF-BC (as observed in the figure). For such an axisymmetric geometry, we may naively think (before any computation) that the long-term evolution of $\omega _x$ (or equivalently that of $L_x$) is unconstrained due to the vanishing component of the viscous torque $\boldsymbol {\varGamma }_\nu \boldsymbol {\cdot } \boldsymbol {1}_x = 0$ according to formula (3.3). We observe that $\omega _x$ initially displays a complicated transient (panel a), which dies out because of viscosity as expected from the asymptotic theory. Then, it converges towards a well-defined oscillatory state after a few viscous time scales (i.e. when $E t \gg 1$ in dimensionless units). The total angular velocity $\boldsymbol {\omega }$, which exhibits no long-term spurious dynamics (panel b), has a small amplitude compared with the mean rotation axis of the fluid $\boldsymbol {\varOmega }=\boldsymbol {1}_z$ with respect to the inertial frame. We have checked that the final state is robust, as it is recovered by varying the numerical resolution and adopting different initial conditions for a few values of $Po$ and $E$ (although multiple solutions may exist close to the inviscid resonances, as shown for sufficiently small Ekman numbers with NS-BC in Cébron 2015).

Figure 2. DNS of precessing ellipsoids with SF-BC at $Po=-1.8, E=5 \times 10^{-4}$ and $P_x = Po \sin (\alpha ) = 10^{-2}$. Axisymmetric geometry $a=1.5$ and $b=c=1$. Time evolution of the Cartesian component $\boldsymbol {\omega } \boldsymbol {\cdot } \boldsymbol {1}_x$ (a) and absolute value $|\boldsymbol {\omega }|$ (b) of the angular velocity, computed in the DNS either from the volume-averaged vorticity as $\boldsymbol {\omega } = (1/2) \int _V \boldsymbol {\nabla } \times \boldsymbol {v} \, \mathrm {d}V$ or using the angular momentum as $\boldsymbol {\omega } = \mathcal {L} \boldsymbol {L}$ using expression (2.13).

The comparison between the asymptotic results and the DNS is further illustrated in figure 3, still considering the illustrative axisymmetric geometry $a=1.5$ and $b=c=1$ (other geometries with $a \neq b$ give again similar results, not shown). The DNS are in excellent quantitative agreement with the asymptotic solution, although the latter has been obtained assuming $E \to 0$, for both the time-averaged and the instantaneous angular velocity (see panel b after seven viscous time scales). We also have checked that $\delta \omega _z^{(2)}$ is accurately recovered in the DNS (not shown). The observed excellent quantitative agreement with theoretical precession-driven flows has not been obtained using NS-BC in ellipsoids, both in DNS (e.g. Noir & Cébron 2013) and laboratory experiments (e.g. Nobili et al. 2021; Burmann & Noir 2022). Finally, the DNS also confirm the physical existence of the two inviscid resonances of solutions (3.7).

Figure 3. Precession-driven flows (SF-BC) at $P_x = Po \sin (\alpha ) = 10^{-2}$ for $a=1.5$ and $b=c=1$. Comparison between asymptotic solution (3.7) and DNS at $E= 5 \times 10^{-4}$. (a) Time-averaged angular velocity $\epsilon = \overline {|\boldsymbol {\omega }|}$ as a function of $Po$. The fluid is not globally rotating when $Po \simeq -1$ if $|P_x| \ll 1$ (grey area). Vertical dashed lines show the two resonances of asymptotic solutions (3.7) at $Po^\pm = - \sqrt {P_x^2 + (1/\lambda _{so}\mp 1)^2}$. Teal vertical line shows the region $|Po|<10^{-2}$ where no $\alpha$ can satisfy $P_x=10^{-2}$. (b) Value of $|\boldsymbol {\omega }|$ as a function of the re-scaled time $E t$ at $Po=-1.8$.

3.3. Asymptotic theory of flow instabilities

The forced laminar flow $\boldsymbol {U} (\boldsymbol {r},t)$, given by (3.4) when $E \ll 1$, can be destabilised by various hydrodynamic instabilities in ellipsoids. Precession-driven instabilities are classified either as viscously driven if they only exist when $E \neq 0$, or as inertial if they survive when ${E = 0}$. Viscous instabilities exist in no-slip spheres, such as boundary-layer instabilities (e.g. Lorenzani & Tilgner 2001; Buffett 2021) or the conical-shear instability (e.g. Lin, Marti & Noir 2015; Cébron et al. 2019). On the contrary, the inertial instabilities only exist in non-spherical geometries (e.g. Kerswell 1993; Wu & Roberts 2011; Vidal & Cébron 2017). In the following, we extend the prior inviscid linear analyses of the inertial instabilities, which all considered precession at $\alpha ={\rm \pi} /2$ and in the precession frame (i.e. only for spheroids), to account for the SF-BC and the time-dependent background flow (3.7) in the mantle frame. To do so, we expand the governing equations with respect to $\boldsymbol {U}$ (discarding the small-amplitude viscous flow $E^{1/2} \tilde {\boldsymbol {U}}$ in the bulk, which is negligible when $E \ll 1$ as found in the DNS). The perturbation velocity $\boldsymbol {v}_{1}$, which is assumed to be of small amplitude compared with $\boldsymbol {U}$, is governed in the mantle frame by

(3.9a)\begin{gather} \partial_t \boldsymbol{v}_{1} + 2 \boldsymbol{\varOmega}_c \times \boldsymbol{v}_{1} = \boldsymbol{\mathcal{L}} (\boldsymbol{v}_{1}) + 2 E \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\epsilon}(\boldsymbol{v}_{1}) - \boldsymbol{\nabla} p, \end{gather}

(3.9b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}_{1} = 0, \end{gather}

with the linearised advection operator $\boldsymbol {\mathcal {L}} (\boldsymbol {a}) = -( \boldsymbol {a} \boldsymbol {\cdot } \boldsymbol {\nabla } ) \boldsymbol {U} - ( \boldsymbol {U} \boldsymbol {\cdot } \boldsymbol {\nabla } ) \boldsymbol {a}$. The perturbation velocity $\boldsymbol {v}_{1}$ then satisfies the SF-BC (Mason & Kerswell 2002; Wu & Roberts 2009)

(3.10a,b)\begin{equation} \boldsymbol{v}_{1} \boldsymbol{\cdot} \boldsymbol{1}_n |_S = 0, \quad [ \boldsymbol{\epsilon} (\boldsymbol{v}_{1}) \boldsymbol{\cdot} \boldsymbol{1}_n] \times \boldsymbol{1}_n |_S = \boldsymbol{0}. \end{equation}

To explore the low-viscosity regime $E \ll 1$, which is difficult to probe using DNS, we develop an asymptotic model. We seek $\boldsymbol {v}_{1}$ using BLT as (e.g. Rieutord 1992)

(3.11a–c)$$\begin{gather} \boldsymbol{v}_{1} (\boldsymbol{r},t) \simeq \boldsymbol{u} (\boldsymbol{r},t) + E^{1/2} \tilde{\boldsymbol{u}} (\boldsymbol{r},t), \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \quad \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{1}_n |_S = 0, \end{gather}$$

where $\boldsymbol {u}(\boldsymbol {r},t)$ represents the inviscid bulk flow and $\tilde {\boldsymbol {u}}(\boldsymbol {r},t)$ is the leading-order viscous flow within the Ekman layer to satisfy SF-BC (3.10). Because the boundary-layer flow has an amplitude that is $E^{1/2}$ smaller than the bulk flow amplitude, SF-BC strongly weaken the viscous instabilities in ellipsoids. In particular, the critical shear layers spawned by the Ekman layer at the critical latitudes are almost suppressed in stress-free ellipsoids without an inner core (Tilgner 1999). Consequently, the inertial instabilities triggered in the (nearly) inviscid bulk are expected to be largely favoured in stress-free ellipsoids (compared with viscous instabilities).

To solve problem (3.11), we introduce the finite-dimensional polynomial vector space $\boldsymbol {\mathcal {V}}_n$ spawned by the global real-valued incompressible elements $\{ \boldsymbol {u}_k\}$, made of Cartesian monomials $x^i y^j z^k$ of maximum degree $i+j+k \leq n$ and satisfying the no-penetration BC (e.g. Vidal, Su & Cébron 2020; Vidal & Cébron 2021a). Such vector elements are indeed known to form a complete basis for smooth velocity fields in ellipsoids when $n \to \infty$ (e.g. Lebovitz 1989; Backus & Rieutord 2017). Then, we seek the bulk flow using the Galerkin expansion (written using Einstein's convention)

(3.12a–c)\begin{equation} \boldsymbol{u}(\boldsymbol{r},t) = \alpha_k (t) \boldsymbol{u}_k (\boldsymbol{r}), \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}_k = 0, \quad \boldsymbol{u}_k \boldsymbol{\cdot} \boldsymbol{1}_n |_S = 0, \end{equation}

where $\boldsymbol {\alpha } = (\alpha _1, \alpha _2, \ldots, \alpha _N)^\top$ is the state vector of the modal coefficients. The number of elements $N$ for a given maximum degree $n$ in expansion (3.12) is $N=n(n+1)(2n+7)/6$. In practice, we truncate the polynomial expansion at the maximum degree $n$, substitute the truncated expansion into (3.9) and, finally, project the resulting equations onto every basis element $\boldsymbol {u}_i$ to minimise the residual with respect to the real-valued inner product defined by $\langle \boldsymbol {a}, \boldsymbol {b} \rangle _V = \int _V \boldsymbol {a} \boldsymbol {\cdot } \boldsymbol {b} \, \mathrm {d} V$. The governing equations then reduce to

(3.13)\begin{equation} \boldsymbol{M} \,\mathrm{d}_t \boldsymbol{\alpha} = ( \boldsymbol{L} - \boldsymbol{C} - \boldsymbol{D} ) \boldsymbol{\alpha},\end{equation}

where $\boldsymbol {M}_{ij} = \langle \boldsymbol {u}_i, \boldsymbol {u}_j \rangle _V$ is the mass matrix, $\boldsymbol {C}_{ij} = \langle \boldsymbol {u}_i, 2 \boldsymbol {\varOmega }_c \times \boldsymbol {u}_j \rangle _V$ represents the Coriolis force, $\boldsymbol {L}_{ij} = \langle \boldsymbol {u}_i, \boldsymbol {\mathcal {L}} (\boldsymbol {u}_j) \rangle _V$ is the matrix representing the linearised advection terms and the viscous matrix $\boldsymbol {D}$ is given by (after integration by parts)

(3.14)\begin{equation} \boldsymbol{D}_{ij} = 2 E \int_V \boldsymbol{\epsilon} (\boldsymbol{u}_i) : \boldsymbol{\epsilon} (\boldsymbol{u}_j) \, \mathrm{d} V,\end{equation}

in which we have enforced SF-BC (3.10) in the projection to simplify the integration (e.g. Guermond et al. 2013). As already noticed for the forced flow, a useful consequence of expansion (3.11) is that the bulk flow $\boldsymbol {u}$ can be determined in (3.13) without an explicit solution of $\tilde {\boldsymbol {u}}$ for SF-BC. This has also been reported for asymptotic models of thermal convection or waves in rotating stress-free spheres (Liao, Zhang & Earnshaw 2001; Zhang & Liao 2004). This is a noticeable difference from asymptotic models using NS-BC, which require a matching between the boundary-layer flow and the interior solution (which are of the same order of magnitude, e.g. Zhang, Liao & Busse 2007; Zhang et al. 2014).

Since asymptotic solution (3.7) is periodic of period $T=2{\rm \pi}$, we investigate the linear stability using Floquet theory. We first compute the eigenvalues $\chi$ of the monodromy matrix $\boldsymbol {\varPhi }(2{\rm \pi} )$ given by

(3.15a,b)\begin{equation} \boldsymbol{M} \, \mathrm{d}_t \boldsymbol{\varPhi} = ( \boldsymbol{L} - \boldsymbol{C} - \boldsymbol{D} ) \boldsymbol{\varPhi}, \quad \boldsymbol{\varPhi}(0) = \boldsymbol{I}, \end{equation}

where $\boldsymbol {I}$ is the identity matrix. Then, we compute the complex-valued Lyapunov exponents as $\mu = (1/T) \log \chi$ whose real part ${\rm Re}(\mu ) = \sigma$ is the growth rate of the instability. As initially noticed by Kerswell (1993) and Wu & Roberts (2011), the finite-dimensional polynomial is left invariant by the linear operator in the momentum equation, that is $\boldsymbol {\mathcal {L}}(\boldsymbol {\mathcal {V}}_n) \in \boldsymbol {\mathcal {V}}_n$. Therefore, we can construct exact polynomial solutions of (3.9) giving sufficient conditions for linear instability in the inviscid regime $E=0$.

We show in figure 4(a) the results of the linear inviscid stability analysis at $P_x = 10^{-2}$. We have numerically solved (3.15) using a fourth-order Runge–Kutta solver and standard linear algebra routines. As in Kerswell (1993) and Wu & Roberts (2011), there are no instabilities associated with the linear elements $n=1$. The first instabilities, which are here associated with the quadratic modes with $n=2$, only occur near the resonance at $Po^+$. When $n$ is increased, additional tongues of inertial (topographic) instabilities appear with a growth rate scaling in the inviscid regime as

(3.16)\begin{equation} \sigma_{topo} = {O}(\epsilon \eta),\end{equation}

where $\epsilon = \overline {|\boldsymbol {\omega }|}$ is the mean value of the differential rotation between the fluid and the mantle and $\eta = a^2/c^2 - 1$ is the polar flattening. The numerical prefactor is found to be $\sigma _{topo} /(\epsilon \eta ) \approx 0.1$ when $n \leq 20$ (as shown in the figure). Moreover, when $n \to \infty$, the growth is expected to approach the upper bound given in the unbounded short-wavelength approximation (Kerswell 1993). This shows that the forced laminar flow is generically unstable to short-wavelength perturbations without viscosity.

Figure 4. Growth rate $\sigma$ of the inertial (topographic) instabilities growing upon flow (3.7) at $P_x=Po \sin (\alpha )=10^{-2}$, as a function of $Po$ (using sampled values). Teal vertical line shows the interval $|Po|<10^{-2}$ in which no $\alpha$ can satisfy $P_x=10^{-2}$. The fluid is not globally rotating near $Po \simeq -1$ when $|P_x| \ll 1$ (grey area). (a) Inviscid growth rate for various degrees $n$ of the global modes. Dashed black curve is obtained in the unbounded short-wavelength limit (Kerswell 1993). (b) Viscous effects. Dotted blue line shows the upper bound of the inviscid growth rate. Olive coloured area shows the unstable region for SF-BC at $E=3 \times 10^{-6}$, and thick red line shows the unstable zone for SF-BC at $E=5 \times 10^{-4}$ (both computed at $n=20$). Purple coloured curves show viscous growth rate (3.17) for NS-BC, with $K \in [4,10]$ to account for the Ekman damping of the large-scale modes (see figure 5).

However, the short-wavelength modes are more damped by viscosity than the large-scale ones. Consequently, viscous effects will select the allowable unstable modes for a given value of the Ekman number. To show this, we have explored the linear stability including viscous damping in figure 4(b). At $E=5 \times 10^{-4}$, the forced flow is only unstable in extremely thin tongues near the two resonances at $Po^\pm$ for the SF-BC. This is consistent with the absence of instabilities in the DNS performed at $E=5 \times 10^{-4}$ (see figure 3). More challenging DNS with SF-BC at smaller values $E = {O}(10^{-6})$, which are beyond the scope of the present paper, could allow us to obtain instabilities for values of $Po$ in a larger interval. Finally, it is also useful to compare the stability of the forced flow with SF-BC and NS-BC. A proper asymptotic theory for the no-slip case, rooted in the BLT of the inertial modes (e.g. Greenspan 1968), will be considered elsewhere. Nonetheless, an upper bound for the viscous growth rate of the inertial instabilities can be estimated as

(3.17)\begin{equation} \sigma_{topo} \approx 0.1 \epsilon \eta - K \sqrt{E [ 1 + P_z]},\end{equation}

assuming that the fluid is rotating on average at $1+P_z$ in the mantle frame. Here, the numerical prefactor $K = 4\unicode{x2013}10$ heuristically accounts for the Ekman damping of the large-scale flow structures with NS-BC (see figure 5). For the small value $E = 3 \times 10^{-6}$, we observe that the forced flow at $P_x = 10^{-2}$ would be mainly stable with NS-BC (except near the resonance $Po^+$), whereas it would be unstable for other values of $Po$ with SF-BC. Therefore, the figure clearly illustrates that adopting SF-BC (instead of NS-BC) can be useful to explore the turbulence driven by inertial instabilities in the bulk of the fluid.

Figure 5. Viscous decay rates $\tau _k$ of the inertial modes of maximum polynomial degree $n=20$, as a function of the inviscid eigenfrequency $\lambda _k \neq 0$. Axisymmetric ellipsoid with semi-axes $a=1.5$ and $b=c=1$ rotating at the angular frequency $\boldsymbol {\varOmega }=\boldsymbol {1}_z$. (a) Complex-valued $\tau _k$ for NS-BC. Colour bar shows the normalised imaginary part ${\rm Im}(\tau _k)$. (b) Real-valued $\tau _k$ (i.e. ${\rm Im}(\tau _k)=0$) given by formula (4.6) for SF-BC.

4. Discussion

4.1. Physical insight from the Coriolis eigenmodes

We have illustrated with the case of precession-driven flows that the long-term evolution of angular momentum is damped by viscosity in triaxial ellipsoids. Similarly, viscosity affects the angular momentum in axisymmetric ellipsoids if the mean rotation axis $\boldsymbol {\varOmega }$ is not aligned with the revolution symmetry axis (even if $\boldsymbol {\varGamma }_i \boldsymbol {\cdot } \boldsymbol {1}_i = 0$ in such geometries, where $\boldsymbol {1}_i$ is the revolution axis along one of the principal semi-axes). Asymptotic analysis offers a physical understanding of why the cases $\boldsymbol {\varOmega } \propto \boldsymbol {1}_i$ and $\boldsymbol {\varOmega } \not \propto \boldsymbol {1}_i$ strongly differ in axisymmetric geometries.

When $E \ll 1$, the solutions of (2.1a) and (2.1b) in stress-free or no-slip ellipsoids can be rigorously expanded onto a combination of the inviscid eigenmodes of the (steady) Coriolis operator given by (e.g. Backus & Rieutord 2017)

(4.1a–c)\begin{equation} \mathrm{i} \lambda_k \boldsymbol{\nabla} \times \boldsymbol{Q}_k ={-}2 \boldsymbol{\nabla} \times ( \boldsymbol{\varOmega} \times \boldsymbol{Q}_k), \quad \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{Q}_k = 0, \quad \boldsymbol{Q}_k \boldsymbol{\cdot} \boldsymbol{1}_n |_S = 0, \end{equation}

where $[\lambda _k, \boldsymbol {Q}_k (\boldsymbol {r})]$ is the $k$th eigenvalue–eigenfunction pair. Only three of these eigenmodes carry a non-zero angular momentum in ellipsoids (by virtue of the orthogonality of the eigenmodes, see Ivers 2017), namely the spin over mode $\boldsymbol {Q}_{so}$, its complex conjugate $\boldsymbol {Q}_{so}^{\dagger}$ and the zero-frequency geostrophic mode $\boldsymbol {Q}_{sup}$ associated with axial (differential) rotation along $\boldsymbol {\varOmega }$. Because these three modes are uniform-vorticity flows such as $\boldsymbol {Q}_k = \omega _{k,x} \boldsymbol {e}_x + \omega _{k,y} \boldsymbol {e}_y + \omega _{k,z} \boldsymbol {e}_z$, they are given by the matrix eigenvalue problem

(4.2)\begin{equation} \begin{pmatrix} 0 & {2 a^2 \varOmega_z}/{(a^2+c^2)} & -{2 a^2 \varOmega_y}/{(a^2+b^2)} \\ -{2 b^2 \varOmega_z}/{(b^2+c^2)} & 0 & {2 b^2 \varOmega_x}/{(a^2+b^2)} \\ {2 c^2 \varOmega_y}/{(b^2+c^2)} & -{2 c^2 \varOmega_x}/{(a^2+c^2)} & 0 \\ \end{pmatrix} \boldsymbol{\omega}_k = \mathrm{i} \lambda_k \boldsymbol{\omega}_k\end{equation}

with $\boldsymbol {\varOmega } = (\varOmega _x, \varOmega _y, \varOmega _z)^\top$, where the rotation vector $\boldsymbol {\omega }_k = (\omega _{k,x},\omega _{k,y},\omega _{k,z})^\top$ of the eigenmode $\boldsymbol {Q}_k$ is given by the $k$th eigenvector of matrix (4.2). Consequently, the uniform-vorticity components $\omega _i (t) \boldsymbol {e}_i$ of the flow in expansion (2.11a,b) are not mutually independent in rotating ellipsoids but, instead, are tied to the dynamics of these modes. More precisely, the equatorial components of the angular momentum $\boldsymbol {L} \times \boldsymbol {1}_\varOmega$ are coupled through the dynamics of the two spin-over modes. Similarly, the axial angular momemtum $\boldsymbol {L} \boldsymbol {\cdot } \boldsymbol {1}_\varOmega$ (related to the fluid spin-up) is piloted by the dynamics of the geostrophic mode $\boldsymbol {Q}_{sup}$.

From a physical viewpoint, whether viscosity affects the long-term evolution of angular momentum or not is thus deeply rooted in the viscous dynamics of these three eigenmodes. We can quantify how viscosity impacts the inviscid eigenmodes by estimating the global viscous decay rates $\tau _k$ of the Coriolis modes as

(4.3)\begin{equation} \partial_{t} \boldsymbol{Q}_k |_{t=0} \simeq \tau_k \boldsymbol{Q}_k. \end{equation}

For NS-BC, $\tau _k$ is a complex-valued quantity with a real part ${\rm Re} (\tau _k) \leq 0$ representing the volume-averaged viscous decay rate, and an imaginary part ${\rm Im}(\tau _k)$ characterising the frequency shift due to viscous effects (e.g. Greenspan 1968). Typical values are illustrated in figure 5 for a particular ellipsoidal geometry. It has also been recognised for a long time that, for NS-BC (2.4), the viscous torque in the mantle frame is related to the viscous damping of these three eigenmodes (e.g. Rochester 1976). In no-slip spherical geometries, it is given by (see formula (35) in Rochester 1976)

(4.4)\begin{equation} \boldsymbol{\varGamma}_\nu \propto E^{1/2} [ {\rm Re}(\tau_{so}) \boldsymbol{\omega}_\perp{-} {\rm Im}(\tau_{so}) \boldsymbol{1}_z \times \boldsymbol{\omega}_\perp{+} \tau_{sup} \omega_z \boldsymbol{1}_z], \end{equation}

at the leading order in $E$ (assuming $\boldsymbol {\varOmega } = \boldsymbol {1}_z$), where $\boldsymbol {\omega } = \boldsymbol {\omega }_\perp + \omega _z \boldsymbol {1}_z = (\omega _x, \omega _y, \omega _z)^\top$ is the uniform vorticity of the forced flow. Note that similar expressions have been later rediscovered for the particular case of precession as viewed in the precession frame (e.g. Noir et al. 2003; Noir & Cébron 2013). Formula (4.4) clearly shows that the equatorial components $\boldsymbol {L} \times \boldsymbol {1}_z$ are damped by viscosity when ${\rm Re}(\tau _{so}) \neq 0$ and, similarly, $\tau _{sup}\neq 0$ (which is a real number for this mode) ensures that the axial angular momentum $\boldsymbol {L} \boldsymbol {\cdot } \boldsymbol {1}_z$ is affected by viscosity. Since ${\rm Re}(\tau _{so}) \neq 0$ and $\tau _{sup}\neq 0$ in no-slip spheres and ellipsoids, we have $\boldsymbol {\varGamma }_\nu \neq \boldsymbol {0}$ from formula (4.4) such that the angular momentum is affected by viscosity on long time scales for the NS-BC.

Similar reasoning can be applied to the stress-free rotating case. It can be shown that leading-order viscous torque (3.3) depends on the viscous decay rates $[\tau _{so},\tau _{sup}]$ for the SF-BC (not given here, since it vainly makes the expression more complex because a full description of the viscous cross-interactions between $\boldsymbol {Q}_{so}$ and $\boldsymbol {Q}_{sup}$ is required contrary to the no-slip case). We can thus get physical insight into formula (3.3) by computing the viscous decay rates for SF-BC. To do so, we expand the velocity as $\boldsymbol {Q}_k + E^{1/2} \tilde {\boldsymbol {Q}}_k$ (Rieutord 1992), where $\tilde {\boldsymbol {Q}}_k$ is the boundary-layer flow such that $\boldsymbol {Q}_k + E^{1/2} \tilde {\boldsymbol {Q}}_k$ satisfies SF-BC (2.3). The viscous decay rate for SF-BC is then given at the leading order in $E$ by (e.g. Liao et al. 2001)

(4.5)\begin{equation} \tau_k \int_V |\boldsymbol{Q}_k|^2 \, \mathrm{d} V = E \int_V \boldsymbol{Q}_k^{\dagger} \boldsymbol{\cdot} \nabla^2 (\boldsymbol{Q}_k + E^{1/2} \tilde{\boldsymbol{Q}}_k) \, \mathrm{d}V. \end{equation}

Contrary to the no-slip case (for which the boundary-layer flow is of the same order of magnitude as the inviscid flow, e.g. Greenspan 1968), an explicit solution of $\tilde {\boldsymbol {Q}}_k$ for SF-BC is not required to estimate $\tau _k$ in (4.5). Indeed, the representative volume-averaged viscous decay rate of all the eigenmodes is given at leading order in $E$ for our SF-BC by (e.g. Rieutord & Zahn 1997)

(4.6)\begin{equation} \tau_k \int_V |\boldsymbol{Q}_k|^2 \, \mathrm{d} V ={-} 2 E \int_V \boldsymbol{\epsilon} (\boldsymbol{Q}_k) : \boldsymbol{\epsilon} (\boldsymbol{Q}_k^{\dagger}) \, \mathrm{d} V. \end{equation}

Expression (4.6) generalises formula (3.14) in Liao et al. (2001), which is only valid for spheres (see Appendix B), to triaxial ellipsoids. Since the right-hand side of (4.5) is real, we have $\tau _k \leq 0$ for SF-BC. Consequently, there is no viscous correction of the inviscid eigenfrequency $\lambda _k$ at the leading order in $E$ for SF-BC (as initially reported in Liao et al. 2001). Formula (4.6) is illustrated in figure 5 for a particular configuration. We recover from the formula that $\tau _{so}=\tau _{sup}=0$ in spherical geometries (since $\boldsymbol {Q}_{so}$ and $\boldsymbol {Q}_{sup}$ are exact solid-body rotations in spheres), which agrees with the fact that $\boldsymbol {\varGamma }_\nu = \boldsymbol {0}$ in spheres (e.g. Jones et al. 2011).

Explicit expressions of $\tau _{so}$ and $\tau _{sup}$ can be obtained for the uniform-vorticity modes in ellipsoids, because the eigenvectors $[\boldsymbol {\omega }_{so},\boldsymbol {\omega }_{sup}]$ of matrix (4.2) can be analytically obtained. The analytical formula of $\tau _{so}$, which is too lengthy to be given here, shows that $\tau _{so} \neq 0$ in every non-spherical geometry. The mathematical reason is that the spin-over mode $\boldsymbol {Q}_{so}$ is no longer a solid-body rotation in ellipsoids (i.e. $\boldsymbol {\epsilon } (\boldsymbol {Q}_k)$ is non-zero for the spin-over mode in ellipsoids). Thus, from a physical viewpoint, a non-zero boundary-layer flow $\tilde {\boldsymbol {Q}}_{so}$ is required to match the SF-BC within a thin Ekman boundary layer. Since the spin-over mode is damped by viscosity in ellipsoids, the equatorial angular momentum $\boldsymbol {L} \times \boldsymbol {1}_\varOmega$ is affected by viscosity on long time scales (even in axisymmetric ellipsoids). After little algebra, the decay rate $\tau _{sup}$ is explicitly given by

(4.7)\begin{equation} \frac{\tau_{sup}}{E} \int_V |\boldsymbol{Q}_{sup}|^2 \, \mathrm{d}V ={-}\frac{16 {\rm \pi}}{3} abc [ \varOmega_x^2 (b^2-c^2)^2 + \varOmega_y^2 (a^2-c^2)^2 + \varOmega_z^2 (a^2-b^2)^2 ], \end{equation}

where the axial geostrophic mode is $\boldsymbol {Q}_{sup} = \omega _{{sup},x} \boldsymbol {e}_x + \omega _{{sup},y} \boldsymbol {e}_y + \omega _{{sup},z} \boldsymbol {e}_z$ with $\omega _{{sup},x} = \varOmega _x (b^2 + c^2), \omega _{{sup},y} = \varOmega _y (a^2 + c^2)$ and $\omega _{{sup},z} = \varOmega _z (a^2 + b^2)$. Formula (4.7) shows that $\tau _{sup} \neq 0$ when $a \neq b \neq c$, illustrating that the axial geostrophic mode is damped by viscosity in triaxial geometries. Therefore, the physical reason why $\boldsymbol {\varGamma }_\nu \neq \boldsymbol {0}$ in triaxial stress-free ellipsoids is that the spin-over and geostrophic modes are damped by viscosity (as evidenced by the non-zero decay rates $\tau _{so} \neq 0$ and $\tau _{sup} \neq 0$ in such geometries). Moreover, formula (4.7) shows that $\tau _{sup}=0$ when $\boldsymbol {\varOmega }$ is an axis of revolution of the geometry (i.e. when $\boldsymbol {\varOmega } \propto \boldsymbol {1}_x$ if $b=c, \boldsymbol {\varOmega } \propto \boldsymbol {1}_y$ if $a=c$, or $\boldsymbol {\varOmega } \propto \boldsymbol {1}_z$ if $a=b$). The axial geostrophic mode is thus unaffected by viscous dissipation, which explains why the long-term evolution of $\boldsymbol {L} \boldsymbol {\cdot } \boldsymbol {1}_\varOmega$ is physically unconstrained in such pathological configurations. This was the situation previously considered for precession-driven flows in spheroids (Lorenzani & Tilgner 2003; Wu & Roberts 2009; Guermond et al. 2013). Yet, the conclusion is not valid for every axisymmetric geometry with global rotation. Indeed, we have $\tau _{sup}\neq 0$ in axisymmetric geometries if $\boldsymbol {\varOmega }$ is not the revolution symmetry axis (such that $\boldsymbol {L} \boldsymbol {\cdot } \boldsymbol {1}_\varOmega$ will be damped by viscosity).

The BLT of Coriolis eigenmodes has thus explained why the long-term angular momentum evolution is damped by viscosity in triaxial geometries, but also in axisymmetric ellipsoids if the mean rotation axis $\boldsymbol {\varOmega }$ is not the revolution symmetry axis.

4.2. Resonance conditions for mechanical forcings

A key property of the primary uniform-vorticity flow is its ability to enter in direct resonance with the precession forcing (as evidenced by the divergent amplitude of the asymptotic solution). A direct resonance requires a close spatial and temporal matching between the Poincaré force and the flow response (e.g. Greenspan 1968). The spatial matching is ensured by the fact that both the Poincaré force and the forced uniform-vorticity flow are linear in the Cartesian coordinates. Heuristically, the temporal resonance condition requires that the frequency $\omega _p$ of the forcing (for monochromatic forcings) must be equal (or close) to the angular frequency $f$ of the forced flow in the mantle frame, which gives $f = \pm \omega _p$. The latter condition generally predicts the existence of two resonances for mechanically driven flows in ellipsoids (if the spatial resonance conditions are satisfied). A quick inspection of (3.4) shows that the uniform-vorticity dynamics roughly corresponds to that of a harmonic oscillator driven by the Poincaré force in the inviscid regime $E=0$. Consequently, direct resonances occur when the forced flow corresponds to a free oscillatory eigenmode of the unforced system, namely the spin-over mode $\boldsymbol {Q}_{so}$ such that $f \propto \lambda _{so}$ (up to a normalisation prefactor). For this reason, longitudinal librations (which only directly excite the zero-frequency geostrophic mode) do not exhibit any inviscid resonance in spheres (e.g. Zhang et al. 2013) or ellipsoids. On the contrary, latitudinal librations can trigger the spin-over mode and the corresponding forced laminar flow exhibits two inviscid resonances occurring at $\lambda_{so} = \pm \, \omega_p$ in non-spherical geometries (Zhang et al. 2012; Vantieghem et al. 2015), where $\omega _p$ is the libration angular frequency. Similarly, a second resonance has already been found for the interaction between tides and precession in triaxial ellipsoids (Cébron, Le Bars & Meunier 2010). A second resonance for pure precession is thus also expected in ellipsoids from simple theoretical arguments. Assuming that the forced uniform-vorticity flow is oscillating in the mantle frame at the effective angular frequency $f \simeq [1 + P_z ] \lambda _{so}$ when $|P_x| \ll 1$, the temporal resonance condition predicts two direct resonances for precession at the resonant Poincaré numbers $Po^\pm$ given by

(4.8)\begin{equation} 1 + Po^{{\pm}} \cos (\alpha) ={\pm} 1/\lambda_{so},\end{equation}

where $\lambda _{so} = 2 ab/\sqrt {(a^2+c^2)(b^2+c^2)}$ is here the eigenfrequency of the spin-over mode in (4.2) with $\boldsymbol {\varOmega }=\boldsymbol {1}_z$ (see also formula 3.21 in Vantieghem 2014). The above condition is exactly the resonance condition of asymptotic solution (3.7). The two resonances at $[Po^-, Po^+]$ are thus robust features of precession-driven flows, but it remains to elucidate why the second resonance at $Po^-$ has not been reported before.

We have numerically solved (3.4) in time to explore the behaviour of the solutions near the double resonances in figure 6. The two resonances at $Po^\pm$ are continuously shifted when $b$ is varied and, at $b=1$, the resonant value $Po^-$ differs from $Po^+$ as observed in panel (b). This directly results from condition (4.8), which predicts that the two resonances are linked by $[Po^+ + Po^-] \cos (\alpha ) = -2$. This clearly shows that the two direct resonances do not merge together in ellipsoids. We further explore the behaviour near $Po^-$ in figure 7. We have fixed the resonant value $Po^-$ at its value given in figure 3 for $a=1.5$ and $b=c=1$ and, then, adjusted the polar axis $c$ to maintain the resonance at $Po^-$ for different values of $a$. We observe that the width of the resonance peak decreases when $a\to b$ (panel a). This is a purely inviscid feature of the asymptotic solution, which is recovered in the DNS. The particular case $a=b$ is not formally defined for SF-BC, but it can be approached by decreasing $a-b$ (panel b). The amplitude of the stress-free solution at $Po=Po^-$ is limited by the viscosity and approaches, when $a \to b$, the inviscid Poincaré solution for $a=b$. The differential rotation $\epsilon$ of the inviscid Poincaré solution is given by (assuming $\omega _z = 1$, see Appendix B in Wu & Roberts 2011)

(4.9)\begin{equation} \epsilon = \left| \frac{P_x (2 + \eta)}{\eta+2(1+\eta)P_z}\right|, \end{equation}

which is non-divergent when $Po=Po^-$. This agrees with a lengthy mathematical analysis of the behaviour near the inviscid resonances (not given here), which shows that the second inviscid resonance at $Po^-$ disappears in spheroids with $a=b$ contrary to the other resonance at $Po^+$ (e.g. Busse 1968; Noir & Cébron 2013).

Figure 6. Double resonance at $Po^\pm$ of the forced precession-driven flow in ellipsoids for SF-BC with $a=1$ and $c=0.9$ (values of $b$ given in the legend). Time-averaged differential rotation $\epsilon = \overline {|\boldsymbol {\omega }|}$ of numerical solutions of (3.4) at $E=10^{-3}$ and small precession angle $\alpha = 3^\circ$. Vertical dashed lines show $Po^{\pm }$ predicted by (4.8) at $b=1$. (a) $Po^{-}$ and (b) $Po^+$.

Figure 7. Behaviour of the forced flow near the second resonance $Po^{-}$ in stress-free ellipsoids with $b=1$ and $P_x=10^{-2}$. The resonant value is fixed at the value $Po^-$ obtained with $a=1.5$ and $c=1$ (as in figure 3). To maintain a fixed resonance when $a$ is varied, the polar axis is given by $c= 0.5 [-2 a^2 - 2 + 2 \sqrt {-32 a^2 \varDelta ^{1/2} + 1 + a^4 + 2 (8 \varDelta + 7) a^2}]^{1/2}$ with $\varDelta = (Po^{-} - P_x)(Po^{-} + P_x)$. In the two panels, the dashed teal line shows the expected inviscid value from Poincaré solution (4.9) for $a=b=1$. (a) Comparison between DNS at $E=5\times 10^{-4}$ and asymptotic solution (3.7). (b) Numerical solutions of (3.4) for SF-BC at $Po=Po^-$ and $E=10^{-3}$.

We have thus understood why precession-driven flows are subject to two inviscid resonances in triaxial ellipsoids, which occur at the resonant Poincaré numbers $Po^\pm$ given by (4.8) when $|P_x| \ll 1$. Since the two resonances are inviscid features of the forced flow in ellipsoids, they exist for both SF-BC and NS-BC. The second resonance actually disappears in spheroidal geometries $a=b$ (i.e. its amplitude is vanishing), which explains why previous works in spheroids have not observed it (e.g. Cébron 2015; Nobili et al. 2021). Previous studies in triaxial geometries (e.g. Noir & Cébron 2013; Burmann & Noir 2022) have also overlooked it, because it usually occurs at $|Po^-| \gg |Po^+|$.

4.3. Implications for DNS

We have shown that the long-term evolution of angular momentum is affected by viscosity, due to the existence of an Ekman boundary layer in rapidly rotating ellipsoids. The uniform-vorticity elements carrying angular momentum in expansion (2.11a,b) do not indeed satisfy the SF-BC in triaxial geometries. Thus, they are associated with an Ekman boundary layer to match the boundary conditions. This is a noticeable difference with the more usual spherical geometry, in which $\boldsymbol {\varGamma }_\nu = \boldsymbol {0}$ (e.g. Jones et al. 2011). The Ekman boundary layer in ellipsoids is clearly observed in figure 8. Its typical thickness is still ${O}(E^{1/2})$ but, contrary to the case of NS-BC, the amplitude of the boundary-layer flow is ${O}(E^{1/2})$ smaller than the bulk flow amplitude (in agreement with Rieutord 1992).

Figure 8. DNS of precession-driven flow with SF-BC at $Po=-1.8$, $P_x = Po \sin (\alpha) = 10^{-2}$ and $E=5 \times 10^{-4}$. Axisymmetric geometry $a=1.5$ and $b=c=1$. Normalised velocity $\boldsymbol{v}_f/(\mathcal{U} E^{1/2})$, as defined in expansion (2.11a,b), at time $t=39\,530$ where $\mathcal{U}=0.0129$ is the maximum of $|\boldsymbol{v}_f|$. (a) Three-dimensional rendering of the velocity magnitude using a linear scale. (b) Axial velocity component as a function of $z$ along the $c$-axis.

This could have implications for numerical studies using stress-free boundaries. A numerical strategy has to be employed to ensure the conservation of angular momentum in spherical codes (e.g. Jones et al. 2011). This is no longer necessary in triaxial ellipsoids since $\boldsymbol {\varGamma }_\nu \neq \boldsymbol {0}$ (albeit such a strategy may be considered to ensure the conservation of the axial angular momentum if the mean rotation axis is an axis of revolution symmetry, as proposed in Guermond et al. 2013). However, for the moderate values of the Ekman number achievable in DNS, the flow within the Ekman layer will modify the value of the viscous torque (which pilots the long-term evolution of angular momentum). Indeed, instead of using expression (2.14), the viscous torque is usually computed with the surface integral $\boldsymbol {\varGamma }_\nu = 2 E \int _S \boldsymbol {r} \times ( \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\epsilon } ) \, \mathrm {d} S$ as

(4.10)\begin{equation} \boldsymbol{\varGamma}_\nu = 2 E \int_S \boldsymbol{r} \times \boldsymbol{\mathcal{T}} \, \mathrm{d} S = 2 E \int_S \boldsymbol{r} \times [ (\boldsymbol{\mathcal{T}} \boldsymbol{\cdot} \boldsymbol{1}_n) \boldsymbol{1}_n ] \, \mathrm{d} S,\end{equation}

in which we have used formula (9) in Rochester (1962) for a symmetric tensor to obtain the first equality and, then, have written the surface traction as $\boldsymbol {\mathcal {T}} = (\boldsymbol {\mathcal {T}} \boldsymbol {\cdot } \boldsymbol {1}_n) \boldsymbol {1}_n - \boldsymbol {1}_n \times (\boldsymbol {1}_n \times \boldsymbol {\mathcal {T}}) = (\boldsymbol {\mathcal {T}} \boldsymbol {\cdot } \boldsymbol {1}_n) \boldsymbol {1}_n$ on the boundary for SF-BC (2.3). Formula (4.10) shows that the normal component of the surface traction, which is non-zero in the presence of an Ekman layer in stress-free ellipsoids, contributes to the viscous torque. Hence, numerical and local approximations of SF-BC (2.3) have no reasons to yield a vanishing torque component in formula (4.10) for axisymmetric ellipsoids if the boundary layer is not sufficiently resolved (as observed in some DNS, not shown). Using a refined boundary-layer mesh may thus be required to properly describe the Ekman layer in ellipsoids and ensure sufficient torque accuracy (which can be used to check the numerical convergence).

4.4. Scaling laws

Despite the existence of a thin Ekman layer, we believe that adopting SF-BC in global simulations is useful to probe bulk mechanisms that can be hampered by viscous effects when NS-BC are employed. The case of precession is illuminating in this respect. Indeed, the laminar precession-driven flow can be destabilised by several hydrodynamic instabilities in no-slip ellipsoids, such as the inertial (topographic) instabilities outlined in § 3.3 and the conical-shear instability (CSI). The former are due to the ellipticity of the boundary and survive in the inviscid regime $E = 0$. On the contrary, the CSI is a parametric instability existing because of the viscous conical shear layers spawned from the Ekman layer at the critical latitudes (Lin et al. 2015). In addition, precession also often triggers boundary-layer instabilities within the Ekman layer for NS-BC (e.g. Lorenzani & Tilgner 2001; Cébron et al. 2019; Buffett 2021). A comprehensive study of these instabilities deserves further work, but we can estimate their relevance as follows. As outlined in § 3.3, the typical inviscid growth rate of the precession-driven inertial instabilities is given by formula (3.16) for the large-scale modes. For the CSI, the growth rate in full spheres and ellipsoids is given by Lin et al. (2015) and Horimoto, Katayama & Goto (2020) as

(4.11)\begin{equation} \sigma_{CSI} = {O}(\epsilon E^{1/5}).\end{equation}

Quantitatively, a necessary condition for the existence of the two instabilities is that growth rates (3.16) and (4.11) are larger than the viscous damping. For the NS-BC, this damping is mainly due to the Ekman layer and its amplitude is of the order ${O}(E^{1/2})$ (Greenspan 1968). Actually, it appears that large-scale inertial instabilities are difficult to obtain for the moderately small values of the Ekman number usually considered in experiments or DNS (as outlined in figure 4).

A linear analysis is, however, not sufficient to determine the physical relevance of these instabilities. In particular, scaling laws are worth finding to estimate the strength of the precession-driven flows driven by such instabilities. Indeed, the inertial instabilities have presumably a saturation amplitude almost independent of the Ekman number (as found for the turbulence driven by tidal instabilities, e.g. Grannan et al. 2017), whereas the CSI amplitude could decrease when $E \to 0$ as the instability results from viscous effects. A rigorous description of the nonlinear regimes requires dedicated simulations, but the saturation amplitudes can be crudely estimated using simple order-of-magnitude arguments (which have already proven useful for tidal flows, e.g. in Barker & Lithwick 2013; Barker 2016b). We assume that the flow amplitude $\mathcal {U}$ resulting from the primary instability grows until secondary instabilities, characterised by the growth rate $\sigma _{sec}$, become strong enough to prevent further growth of the primary instability. A saturated turbulent regime would then be obtained when $\mathcal {U} \sim \sigma _{sec} \ell$, where $\ell$ is a characteristic length scale of the primary unstable flow. The nonlinear saturation of the inertial (topographic) instabilities would thus be given by (in dimensionless units)

(4.12)\begin{equation} \mathcal{U}_{topo} = {O}(\epsilon \eta),\end{equation}

with $\ell \sim 1$ for a large-scale instability. A good agreement with the above scaling law has been found using DNS in shearing periodic boxes (Barker 2016b), but the scaling law might be different for short-wavelength instabilities with $\ell \ll 1$. Using the same reasoning for the CSI, the relevant length scale is likely the width of the critical shear layer $\ell \sim E^{1/5}$ (Lin et al. 2015). Assuming that the CSI is limited by secondary CSI within the critical shear layers, we obtain the (dimensionless) scaling law

(4.13)\begin{equation} \mathcal{U}_{CSI} = {O}(\epsilon E^{2/5}).\end{equation}

We compare in figure 9(a) the above scaling law with previously published DNS in no-slip full spheres (Lin et al. 2015; Cébron et al. 2019). Considering the full sphere geometry allows us to discard the possible CSI resulting from the inner boundary (which would give a different scaling). However, even in the full sphere, identifying the instability mechanism is difficult due to the competition between the CSI and the boundary-layer instabilities. Moreover, due to the non-trivial dependence of the two viscously driven instabilities with the forcing parameters, it is unlikely that a single scaling law could fully describe the entire simulation dataset. The onset distance is indeed difficult to estimate (e.g. see figure 6 in Cébron et al. 2019). Besides, the simulations may not be in the asymptotic regime $E \ll 1$. Despite such uncertainties, a fairly good agreement is found between the DNS and scaling law (4.13) sufficiently far from the onset. This suggests that the CSI was present in the nonlinear regime and that its saturation amplitude obeys scaling law (4.13) for sufficiently small values of the Ekman number.

Figure 9. (a) Comparison between scaling law (4.13) and DNS for $|Po| < 0.1$ in no-slip full spheres from the database of figure 7 in Cébron et al. (2019), with $E_f = E/|1+Po| \simeq E$ when $|Po| \ll 1$. Colour bar indicates $\log _{10} |Po|$. Grey area shows the scaling law $\mathcal {U} E^{-2/5} = (6.5 \pm 1.5) (\epsilon -\epsilon _c)$ and dashed line $\mathcal {U} E^{-2/5} = 6.5 (\epsilon -\epsilon _c)$, where $\epsilon _c \approx 7 E_f^{3/10}$ is an estimate of the onset value (see equation (17) in Cébron et al. 2019). (b) Competition between the inertial instabilities and the CSI in precessing ellipsoids. Empty blue squares ${\square }$ (blue) show conditions for which inertial (topographic) instabilities are expected, and red crosses ${\times }$ (red) indicate where the CSI is expected or observed. Grey area shows $\eta \propto E^{2/5}$ with a (unknown) prefactor chosen in the range $[1,100]$, in which we expect $\mathcal {U}_{topo} \sim \mathcal {U}_{CSI}$ when both instabilities exist. Hatched area is the region where $\sigma _{topo} \gtrsim \sigma _{CSI}$ (if the two instabilities coexist). White area is the region where $\mathcal {U}_{topo} \ll \mathcal {U}_{CSI}$. Estimates for the early Moon and Earth taken from Appendix C with $\eta \approx 2 f$.

Finally, the comparison between scaling laws (4.12) and (4.13) shows that the inertial instability would have a larger amplitude than the CSI when $\eta \gg E^{2/5}$ (if the two instabilities were simultaneously triggered). The resulting regime diagram is illustrated in figure 9(b), using planetary estimates given in Appendix C. Precession-driven inertial instabilities may only have been excited in the primitive liquid cores of the Earth and Moon, whereas the CSI is expected to be present (respectively absent) in the core of the Moon (respectively the Earth) during its whole history (Lin et al. 2015; Landeau et al. 2022). In the early Moon, the inertial instabilities may have dominated the CSI in flow amplitude (although the CSI may have had a larger growth rate than the inertial instabilities according to previous formulas, not shown). Therefore, the inertial instabilities may actually be more relevant than the CSI for some planetary conditions (although they have not been convincingly observed yet in experiments, e.g. Nobili et al. 2021; Burmann & Noir 2022). This could be key for the generation of planetary magnetic fields, as initially postulated for the geodynamo (Malkus 1968). Preliminary estimates of the dynamo capability of the precession-driven instabilities, obtained using (speculative) order-of-magnitude arguments, are given in Appendix C.

5. Conclusion

5.1. Summary

We have investigated precession-driven flows in stress-free ellipsoids, using asymptotic analysis and targeted DNS. We have developed a reduced model for SF-BC to determine the forced uniform-vorticity flows, which carry angular momentum. We have shown that angular momentum is affected on long time scales by viscosity in triaxial ellipsoids, but also in axisymmetric geometries if the mean rotation axis is not a revolution symmetry axis. This is a noticeable difference from spherical geometries, in which angular momentum is unaffected by viscosity. The fundamental reason is that the flows carrying a non-zero angular momentum in ellipsoids are associated with an Ekman boundary layer in rotating ellipsoids. From a numerical viewpoint, a boundary-layer mesh may be necessary to get numerical convergence of the angular momentum in rotating ellipsoids. We also have obtained the analytical solution of the time-dependent laminar flow forced by precession in the mantle frame, which is valid for planetary parameters and triaxial geometries. The comparison with the DNS has shown that, even for moderately small values of the Ekman number, the forced laminar flow in the DNS converges to the asymptotic solution in the vanishing viscosity regime. Moreover, we have uncovered a second (inviscid) resonance of the forced laminar flow in triaxial ellipsoids.

Then, we have explored the inertial instabilities growing upon the forced laminar flow in the bulk, which survive in the inviscid regime $E=0$. We have shown that these instabilities could be more easily observed in stress-free ellipsoids than in no-slip ones (at least for the moderate values $E \gtrsim 10^{-6}$ considered in DNS). We have finally proposed scaling laws for the velocity amplitude of the inertial instabilities and of the CSI, which are in good agreement with previous DNS. The comparison between the two scaling laws confirms that replacing NS-BC with SF-BC in the mantle frame could be useful to directly probe scenarios of bulk turbulence in the low-viscosity regime (which are of interest for planetary modelling).

5.2. Perspectives

Despite the presence of a thin Ekman boundary layer, we believe that SF-BC are relevant for global models of mechanically driven flows. The stress-free model could be used to investigate the saturated flows driven by the inertial (topographic) instabilities in precession ellipsoids and, then, their dynamo capability for planetary applications (as outlined in Appendix C). Stress-free models could indeed shed new light on alternative mechanisms giving birth to dynamo fields in planetary interiors. For instance, the past dynamo of the Moon may have been driven by precession (e.g. Dwyer, Stevenson & Nimmo 2011). Yet, previous numerical investigations of precession-driven dynamos failed to reproduce large-scale magnetic fields in spherical geometries (Cébron et al. 2019). This could result from the fact that the turbulence was driven in those simulations by viscous flows (e.g. the CSI or boundary-layer instabilities), which may be negligible in amplitude compared with the turbulence driven by the inertial (topographic) instabilities in the early Moon (as discussed in § 4.4). This hypothesis could be tested in simulations using stress-free ellipsoids. Similarly, energetic arguments suggest that the dynamo of the early Earth may have been sustained by tidal flows (Landeau et al. 2022). However, the associated fluid dynamics remains to be quantitatively studied to go beyond prior proof-of-concept simulations (Reddy et al. 2018; Vidal et al. 2018). Precessing stress-free ellipsoids are also relevant for short-period hot Jupiters (Barker 2016b), or gaseous planets with a big moon outside the equatorial plane (e.g. the Neptune/Triton pair, Wicht & Tilgner 2010).

Finally, SF-BC could also be used to revisit the long-standing problem associated with the generation of geostrophic flows in rotating fluids (Greenspan 1969). Nonlinear interactions within the Ekman boundary layers for NS-BC (e.g. Busse 1968; Cébron et al. 2021) or in the bulk through the action of the Reynolds stresses (e.g. Zhang & Liao 2004; Livermore et al. 2016), are usually invoked, but geostrophic flows can also result from bulk turbulence. However, it remains unclear whether two- or three-dimensional rotating bulk turbulence is established in natural systems (e.g. Le Reun et al. 2019). This fundamental problem has been attacked in cylindrical or plane-layer geometries (e.g. Kerswell 1999; Brunet, Gallet & Cortet 2020; Le Reun et al. 2020). The latter geometries are, however, not directly relevant for planetary modelling, due to the absence of the so-called topographic beta effect that strongly modifies the geostrophic flows in spheres and ellipsoids (e.g. Greenspan 1968). We believe that using SF-BC opens the way for new fundamental studies dealing with the interplay between waves and geostrophic flows in global geometries.