The present numerical study aims at shedding light on the mechanism underlying the precessional instability in a sphere. Precessional instabilities in the form of parametric resonance due to topographic coupling have been reported in a spheroidal geometry both analytically and numerically. We show that such parametric resonances can also develop in spherical geometry due to the conical shear layers driven by the Ekman pumping singularities at the critical latitudes. Scaling considerations lead to a stability criterion of the form |*P _{o}*| >

*O*(

*E*

^{4/5}), where

*P*represents the Poincaré number and

_{o}*E*represents the Ekman number. The predicted threshold is consistent with our numerical simulations as well as previous experimental results. When the precessional forcing is supercriticial, our simulations show evidence of an inverse cascade, i.e., small scale flows merging into large scale cyclones with a retrograde drift. Finally, it is shown that this instability mechanism may be relevant to precessing celestial bodies such as the earth and earth’s moon.

## I. INTRODUCTION

Precession corresponds to the gyroscopic motion of a rotating object due to a torque orthogonal to its spin axis. When an internal liquid layer is present, such as the liquid core or the subsurface ocean of a planet or the fuel tank of a spinning spacecraft, energy can be dissipated in the liquid due to the induced flows. In line with this idea, it has been proposed that orbital perturbations, in particular precession and nutation, could be used to probe the interior of the planet^{1–3} and possibly generate magnetic fields through the so-called dynamo process.^{4–8}

It is well-established that the primary response in the bulk of a rapidly rotating fluid cavity subject to precession is of uniform vorticity, i.e., the fluid rotates along an axis tilted with respect to the mean rotation axis of the container.^{2,9–12} Although viscosity plays an important role mostly in the viscous boundary layer next to the solid wall of the cavity, it has been shown that singularities at the so-called critical latitudes excite conical shear layers in the interior that are superimposed on the uniform vorticity flow.^{13–16} In addition, weakly non-linear interactions in these singular regions of the boundary layer generate steady geostrophic shears,^{9} which were observed both numerically^{17–19} and experimentally.^{10,20–25}

In a pioneering piece of experimental work, Malkus^{20} showed that precession driven flows can become unstable and even create space filling turbulence. These initial findings were re-confirmed by later experiments.^{10,21,23,26} A theoretical justification was first proposed by Kerswell,^{27} who argued that the precession driven uniform vorticity flow in a spheroidal cavity is inertially unstable due to the constant strain field exerted by the solid wall. Indeed, in a spheroid, the uniform vorticity flow is a superposition of a solid body rotation and a gradient flow, which is required to fulfil the non-penetration condition at the wall. This gradient flow can be decomposed into two components, one leading to elliptically deformed streamlines and the other to a shear of the centres of the streamlines. Both parts can interact with a pair of free inertial modes of the rotating cavity through parametric couplings. This mechanism has been confirmed numerically^{28} for oblate spheroids, while similar inertial instabilities were reported experimentally in cylindrical precessing tanks.^{29–31}

Both numerical simulations^{18,32,33} and laboratory experiments^{23,26} in spherical geometries have shown very rich dynamics ranging from laminar flows to fully developed turbulence. However, in the case of a spherical cavity, the uniform vorticity solution reduces to a purely solid body rotation, preventing the aforementioned inertial instability to develop. Hence, instabilities in a sphere can only originate from the viscous correction to the solid body rotation flow.^{34}

The present study aims at shedding light on the onset of the unstable flows in a precessing sphere. Our numerical results show that the shear in the conical structures spawned from the critical latitudes is key to the destabilization process, echoing the shear of the centres of the streamlines in the case of a spheroidal cavity. We conjecture that the conical shear layers couple with two free inertial modes leading to a parametric resonance. The known scalings of the conical shear layers allow us to derive a stability criterion for the onset of the unstable flow that is in agreement with our numerical simulations. Furthermore, we show that for large enough forcing, small scale vortices merge into large scale cyclonic vortices with a retrograde drift.

The paper is organized as follows. We first introduce the mathematical background and the numerical model in Sec. II, and then, we present our numerical results in Sec. III. Using heuristic arguments, we derive a scaling for the onset of the observed instability in Sec. IV. Finally, we discuss our findings in the context of planetary dynamics in Sec. V.

## II. MATHEMATICAL BACKGROUND AND NUMERICAL METHOD

We consider a sphere of radius *R* filled with a homogeneous and incompressible fluid of density *ρ* and kinematic viscosity *ν*. The sphere rotates at $ \Omega o = \Omega o k \u02c6 $ and precesses at $ \Omega p = \Omega p k \u02c6 p $, where $ k \u02c6 $ and $ k \u02c6 p $ are unit vectors along the spin and precession axes, respectively (Figure 1). Hereinafter, we refer to the mantle frame as the frame attached to the container, the precession frame as the frame rotating at **Ω**_{p} and the fluid frame, the frame attached to the solid body rotation of the fluid. In this paper, all the numerical calculations are implemented in the mantle frame, while the precession frame and the fluid frame are used from time to time to discuss the results.

### A. Equations

Using the radius *R* as the length scale and $ \Omega o \u2212 1 $ as the time scale, the dimensionless Navier-Stokes equations describing the fluid velocity ** u** and pressure

*p*in the mantle frame become

^{35}

where the time-dependent precession vector $ k \u02c6 p $ is given by

Here, $ ( \u2009 i \u02c6 , j \u02c6 , k \u02c6 ) $ are unit vectors of Cartesian coordinates (*x*, *y*, *z*) whose *z*-axis is along the rotation vector $ k \u02c6 $. In Eq. (1), the last term on the left hand side is the Coriolis force due to the rotation and precession, the last term on the right hand side is the so-called Poincaré acceleration, and *p* is the reduced pressure.

In the present study, the angle *α _{p}* between the spin and precession axes is kept constant and equal to 60°. The evolution of the system is then governed by two dimensionless parameters, the Poincaré number

*P*and the Ekman number

_{o}*E*,

which measure the dimensionless rate of precession and the ratio of the typical viscous force and Coriolis force, respectively. Negative (positive) values of *P _{o}* correspond to retrograde (prograde) precession; we consider only retrograde precession in the present study.

### B. Numerical solver

Equations (1) and (2), together with no slip boundary condition ** u** = 0 on the wall, are solved using the fully spectral code developed by Marti.

^{36}Using the so-called toroidal/poloidal decomposition in a spherical coordinate system (

*r*,

*θ*,

*ϕ*)

the velocity field is then represented by two scalar fields *T* and *P*. The incompressibility condition is automatically satisfied by such a decomposition. The scalar fields *T* and *P* are then expanded as

where $ Y l m ( \theta , \varphi ) $ are the spherical harmonics of degree *l* and order *m* and $ W n l ( r ) $ are the so-called Worland polynomials. For a given harmonic degree *l*, the Worland polynomials are a combination of the *r ^{l}* factor and the Jacobi polynomials $ P n ( \alpha , \beta ) ( x ) $, i.e., $ W n l ( r ) = r l P n \u2212 1 / 2 , l \u2212 1 / 2 ( 2 r 2 \u2212 1 ) $, which exactly satisfy the parity and regularity at the origin of the sphere.

^{37}

By using the decomposition from Eqs. (5)–(7), the radial component of the ∇ × Eq. (1) and ∇ × ∇ × Eq. (1) yield decoupled evolution equations for the spectral coefficients $ t l , n m $ and $ p l , n m $. The time integration is implemented using a second order predictor-corrector scheme.

The numerical code has been widely benchmarked in several contexts including that of precession driven flows.^{32,38} We use a typical truncation of the spectral expansion up to *N* = 63, *L* = *M* = 127 at moderate Ekman numbers (*E* ⩾ 3 × 10^{−5}). For a few calculations at lower Ekman numbers (*E* ⩽ 10^{−5}), we truncated as high as *N* = 255, *L* = 511, and *M* = 255 to well resolve the thin boundary layer and instabilities.

### C. Derived quantities

In this subsection, we introduce some useful quantities derived from our numerical simulations performed in the mantle frame.

The solution to the linearized equations (1) and (2) is symmetric around the origin due to the parity of the precessional forcing^{32,35} (see also in Appendix A), namely, ** u**(−

**) = −**

*r***(**

*u***), where**

*r***is the position vector. Any breaking of the symmetry must be caused by an instability. Hence, it is of interest to decompose the total velocity into its symmetric and antisymmetric parts, i.e.,**

*r***=**

*u*

*u*_{s}+

*u*_{a}with

Such a decomposition can easily be performed in the spectral space by taking advantage of the parity of the spherical harmonics. Specifically, the symmetric part includes only odd degree *l* in the toroidal field *T* and even degree *l* in the poloidal field *P*, while the antisymmetric part is the other way round.

Following the symmetry decomposition, we define the total energy, symmetric energy, and antisymmetric energy as

The growth of an anti-symmetric energy component is a sufficient condition to identify the development of an instability. However, it should be noted that a symmetry-breaking is not a necessary condition, indeed instabilities with negligible anti-symmetric component have been reported in numerical studies.^{28} In all cases, the instabilities will lead to time variations of the energy in the system that is otherwise steady.

The main observations from previous experimental studies are reported in the precession frame;^{10,20,21} it is therefore natural to provide the reader with an equivalent point of view in our numerical simulations. To that end, when needed, after we perform our simulations in the mantle frame, velocities are simply transformed in the precession frame as

Anticipating the peculiar role of the conical shear layers which are coaxial with the rotation axis of the fluid,^{10,21} we extract the rotation vector of the fluid *ω*_{F} from the mean vorticity in the bulk^{18}

where < ⋅ > denotes averaging in the fluid volume excluding a thin boundary layer (10*E*^{1/2}).

Once *ω*_{F} is determined, the spectral coefficients with respect to the system of coordinates aligned with the rotation axis of the container are converted into a new set associated with the coordinate system aligned with *ω*_{F}, from which the velocities are reconstructed. This transformation is made using a Matlab subroutine^{39} based on the formula in Ref. 40. The angle *α _{f}* between the rotation axes of the container and the fluid is defined as $cos \alpha f = k \u02c6 \u22c5 \omega F $. Finally, $\epsilon =| k \u02c6 \u2212 \omega F |$ represents the differential rotation between the container and the bulk of the fluid.

## III. RESULTS

### A. Base flow

It is now well established that the base flow in a weakly precessing sphere is in the form of a solid body rotation along an axis inclined to the rotation axis of the container.^{9,16,35} Figure 2 represents the angle *α _{f}* between the rotation axes of the container and the fluid as a function of the Poincaré number

*P*and Ekman number

_{o}*E*. We observe a quantitative agreement between our simulations (symbols) and the theoretical predictions

^{9}(solid lines) that have been thoroughly validated.

^{10,18}We note in particular that the scaling

*α*∝

_{f}*E*

^{−1/2}at fixed

*P*expected from the resonance at |

_{o}*P*|cos

_{o}*α*≪

_{p}*E*

^{1/2}is well recovered.

^{9}Figure 2 also includes data points for which an instability occurs (▴), showing that even in these cases the asymptotic theory of Busse

^{9}remains valid at first order.

A dominant feature of the viscous corrections to the solid body rotation flow is the conical shear layers spawned from the critical latitudes. These structures, coaxial with *ω*_{F}, correspond to the so-called characteristic surfaces of the hyperbolic inertial wave equation, i.e., the unforced inviscid Navier-Stokes equations in the rotating frame. The evidence of such oblique shear layers has been reported in numerical and experimental studies,^{17,19,22} with a typical velocity of *O*(ε*E*^{1/5}) over a typical width of *O*(*E*^{1/5}).^{14,16,41} A typical example at *P _{o}* = − 1.0 × 10

^{−4}and

*E*= 1.0 × 10

^{−6}is represented in Figure 3, showing the conical shear layers in the fluid frame. As we shall see later on in this paper, they play an important role in the destabilization mechanism of the flow.

### B. The parametric instability regime

As we increase the precession rate at a given Ekman number *E*, the base flow becomes unstable. A typical example just above the threshold is presented in Figs. 4–6 for *P _{o}* = − 7.0 × 10

^{−3}and

*E*= 3.0 × 10

^{−5}. Figure 4 represents the time evolution of the antisymmetric energy

*Ek*in the mantle frame. After a transient stage, 0 <

_{a}*t*/2

*π*< 200, the antisymmetric energy

*Ek*grows exponentially until saturation is reached at

_{a}*t*/2

*π*∼ 1000.

Figure 5 shows a snapshot (*t*/2*π* = 973, just before the saturation) of the antisymmetric velocities *u _{ra}* and

*u*

_{θa}in the fluid frame, on a spherical shell at 10

*E*

^{1/2}below the surface (a) and (d), in a meridional cross section (b) and (e), and in the equatorial plane perpendicular to

*ω*_{F}(c) and (f). The antisymmetric flow is mostly confined between two cylinders (0.6 <

*s*< 0.8) co-axial with

*ω*_{F}, which correspond to latitudes higher than the critical latitudes in the Ekman layer at ±30°. In the equatorial plane perpendicular to

*ω*_{F}, the antisymmetric velocities are dominated by

*m*= 17 for

_{F}*u*and

_{ra}*m*= 18 for

_{F}*u*

_{θa}, where

*m*is the azimuthal wavenumber with respect to

_{F}

*ω*_{F}. This is consistent with the symmetry of the velocity components, indeed for an antisymmetric flow,

Hence, only modes with odd *m _{F}* contribute to

*u*and modes with even

_{ra}*m*to

_{F}*u*

_{θa}in the equatorial plane.

Figure 6 shows the contribution of individual *m _{F}* modes to the symmetric and antisymmetric energy at the same instant and parameters as in Figure 5. The symmetric energy

*Ek*is dominated by the

_{s}*m*= 1 component, whereas the antisymmetric energy is dominated by the

_{F}*m*= 17 and

_{F}*m*= 18 components. These observations in the spectrum are consistent with the observed velocities in Figure 5. The less significant peaks at higher wavenumbers may be attributed to secondary bifurcations.

_{F}^{34}

In Figure 7(a), we extract time series of the velocities corresponding to *m _{F}* = 17 and

*m*= 18 at a fixed position in the fluid frame during the growth phase of the antisymmetric energy. The associated Discrete Fourier Transform (DFT) is shown in Figure 7(b) where frequencies are normalized by |

_{F}

*ω*_{F}|. The animation (Movie 1 in the supplementary material

^{50}) of the velocities in the equatorial plane indicates that the

*m*= 17 mode is prograde and the

_{F}*m*= 18 mode is retrograde. So the frequencies of the two modes are

_{F}*ω*

_{17}/|

*ω*_{F}| = − 0.34 and

*ω*

_{18}/|

*ω*_{F}| = 0.67, respectively, which satisfy

*ω*

_{18}/|

*ω*_{F}| −

*ω*

_{17}/|

*ω*_{F}| ≈ 1.0. We also see some small amplitude high frequency fluctuations in the time series which may be due to uncertainties in the determination of

*ω*_{F}.

The observed azimuthal wavenumbers and frequencies are consistent with a parametric resonance mechanism similar to the shear instability described by Kerswell^{27} in precessing spheroids. Such an instability can arise when the background flow couples with two free inertial modes that satisfy the so-called parametric resonant conditions, *ω*_{2} − *ω*_{1} = 1.0, *m*_{2} − *m*_{1} = 1, and *l*_{2} = *l*_{1}, where *ω* is the eigen-frequency, *m* is the azimuthal wavenumber, and *l* is the spherical harmonic degree of the inertial modes (see Appendix C). To fully characterize the possible inertial modes interacting in our simulations, we calculate the eigen-frequencies and the velocity structures of the inviscid inertial modes in a sphere following the analytical approach in Ref. 42 (see also Appendix C). Among all possible inertial modes, we only have to consider the antisymmetric ones (odd *l*) with azimuthal wavenumber *m* = 17 and *m* = 18 based on our observations. In Figure 8, we plot the eigen-frequencies *ω*_{17} + 1.0 and *ω*_{18} as a function of *l* for inertial modes with *m* = 17 (×) and *m* = 18 (∘). We note that there are several combinations nearly satisfying *ω*_{18} − *ω*_{17} = 1.0, i.e., collocated symbols in Figure 8. However, only one combination (in the dashed ellipse), with *l* = 27, is found to match the observed frequency (the dashed line).

Figure 9 compares the velocity structures of the observed unstable modes (a) and (b) in the numerical simulations with the identified inviscid inertial modes (c) and (d) matching the observations. We observe that the inviscid inertial modes have the same structure in the unstable region 0.6 < *s* < 0.8 as in our numerical simulations. In the most outer region, the agreement is more qualitative as we do not observe significant flows in our simulations while the inviscid modes exhibit a well defined pattern with significant velocities. This discrepancy may be attributed to viscous effect that is more influential in this region.

### C. Transition to turbulence

As we increase the precession rate well above threshold, the flow exhibits an inverse cascade of the vorticity. This is illustrated in Figure 10, which shows a sequence of snapshots of the axial vorticity in the equatorial plane and the axisymmetric azimuthal velocity in a meridional section in the fluid frame at *P _{o}* = − 1.35 × 10

^{−2}and

*E*= 3.0 × 10

^{−5}. At the beginning (Figure 10(a)), we see small scale vortices similar to the flow observed during the growth of the parametric instability in Figure 5. The small scale structures then start to merge into large scale elongated cyclonic structures closer to the rotation axis of the fluid (Figures 10(b)–10(d)). In addition, we observed a retrograde (westward) drift of the pattern, as seen in the animation (Movie 2 in the supplementary material

^{50}). These observations are consistent with an inverse cascade of energy characteristic of two dimensional turbulence. Although not yet formally established, the retrograde drift seems to result from the conservation of angular momentum.

In the case of Figure 10, the large scale cyclones break down to small scales and form again (Figures 10(e) and 10(f)); in some cases, the cyclones can be sustained for more than 200 rotation periods of the container. Similar large scale vortices were also observed experimentally in a precessing cylinder^{43} and numerically in rotating Rayleigh Bénard convection,^{44} yet the underlying merging mechanism remains poorly understood.

### D. Stability diagram

Figure 11 shows the regime diagram in the (*P _{o}*,

*E*)-parameter space accessible in our numerical study. The flow is characterized as stable (open circles) when it is centrosymmetric and has a steady total kinetic energy after the transient stage. In contrast, symmetry-broken flows (filled triangles) or centrosymmetric flows but with time-varying total kinetic energy (open triangles) indicate the presence of an instability. We note that although the unstable flows are antisymmetric for the most part, we observe a few cases where the velocity field remains centrosymmetric, as for instance at

*P*= − 0.06 and

_{o}*E*= 10

^{−4}. We also observe that the flow in a precessing sphere can be stable when the precession is sufficiently strong, as previously reported in laboratory experiments.

^{23,26}In these cases, the fluid axis is nearly aligned with the precession vector.

The solid line in Figure 11 represents the lower instability threshold of laboratory experiments in a precessing sphere.^{45} We can see that our numerical simulations are in good agreement with the experimental results particularly when *E* ≤ 10^{−4}. Finally, both numerical and experimental results are consistent with a critical Poincaré number scaling as |*P _{o}*|sin

*α*=

_{p}*O*(

*E*

^{4/5}) that we shall now derive on the basis of a parametric instability mechanism.

## IV. THE CONICAL SHEAR-DRIVEN PARAMETRIC INSTABILITY (CSI)

Our observations are consistent with a parametric instability mechanism similar to the shear instability described by Kerswell.^{27} However, in contrast with the spheroidal geometry, the shear cannot result from topographic effects. Instead, we argue that the conical shear layers spawned from the critical latitudes in the boundary layers can also induce a parametric instability. Indeed, the conical shear layers driven by the linear viscous interactions in the boundary layer are *m _{F}* = 1 and

*ω*/|

*ω*_{F}| = 1.0 in the fluid frame and steady in the precession frame. Thus, they satisfy the parametric resonant conditions with the observed modes in our simulations. Hereinafter, we will refer to this instability as a CSI, for Conical-Shear-Instability.

Figure 12(a) shows the streamlines in the frame of precession just before the growth of the unstable modes illustrating the distortion along the conical shear surfaces, here represented in light grey. Figure 12(b) shows the traces of the conical shear layers (left half) and velocity vectors (right half) in the meridional plane (**Ω**_{o}, ** ω**).

Following the well established theory,^{46,47} we shall now derive the scaling for the onset condition. For an exact resonance, the growth rate of the parametric instability can be written as

where *κ*_{1} and *κ*_{2} are the real decay rate factors of two inertial modes. The interaction parameters between the conical shear *u*_{shear} and two inertial modes *u*_{1}, *u*_{2} are given by

with <** A**,

**> = ∭**

*B*_{V}

*A*^{∗}⋅

**d**

*B**V*and

*A*^{∗}is the complex conjugate of

**.**

*A*A complete derivation of the growth rate requires evaluation of the volume integrals in Eqs. (14) and (15). Although an explicit expression of *u*_{shear} was derived recently,^{16} a set of partial differential equations have to be solved numerically to get *u*_{shear}. Nevertheless, a scaling law for the onset of the parametric instability can be established using heuristic arguments. It has been demonstrated that the amplitude of the conical shear layers scales as |*u*_{shear}| = *O*(ε*E*^{1/5}) over a width *O*(*E*^{1/5}), where ε represents the differential rotation between the fluid and the surrounding solid shell.^{14,16,19} Therefore, we estimate |∇ × *u*_{shear}| = *O*(ε*E*^{1/5}/*E*^{1/5}) = *O*(ε), while the integration volume is proportional to *E*^{1/5}. It follows

leading to

In a precessing sphere, the differential rotation ε is a function of *P _{o}*,

*α*, and

_{p}*E*. At low precession rate, i.e., |

*P*|sin

_{o}*α*≪

_{p}*E*

^{1/2}, the direct resonance mechanism between the tilt-over mode and the precessional forcing leads to ε =

*O*(|

*P*|sin

_{o}*α*/

_{p}*E*

^{1/2}).

^{9,16}This scaling is also confirmed by our numerical simulations in Figure 2. Thus, the lower bound for the threshold of the CSI in a sphere can be written as

This scaling is in quantitative agreement with our numerical simulations (Figure 11) and previous experimental results.^{45}

The conical shear layers can also be excited in precessing spheroidal cavities^{48} and thus the CSI should be induced as well. In this geometry, we must calculate the differential rotation ε following Busse^{9} which also depends on the ellipticity *η* of a spheroid and then apply the instability criterion of Eq. (18).

## V. DISCUSSION AND CONCLUDING REMARKS

The present numerical study investigates the stability of precession-driven flows in a full sphere at moderate to low Ekman numbers. It is found that at low precession rate the flow is of uniform vorticity with a viscous correction superimposed. As the precession rate increases, the internal conical shear layers driven by the linear interactions in the boundary layer induce a parametric instability. The threshold conditions have been established using heuristic arguments leading to a scaling |*P _{o}*|sin

*α*=

_{p}*O*(

*E*

^{4/5}) in quantitative agreement with both our numerical simulations and former experimental observations. At onset, the inertial modes involved in the destabilization mechanism concentrate the energy in an annular region between two cylinders coaxial with the rotation axis of the fluid (0.6 <

*s*< 0.8). Above threshold, the flow evolves from its onset geometry with relatively large azimuthal wavenumbers to a simpler structure with few cyclonic vortices, closer to axis of rotation of the fluid and traveling in a westward direction.

Our results may allow us to shed light on the long standing ill-understood wave-like instabilities reported by Vanyo *et al.*^{21} almost 20 yr ago in a precessing spheroid with ellipticity *η* = 1/100. In one of their experiments, dye was injected next to the boundary layer of a precessing fluid cavity. Pictures taken over more than one hour of experimentation reveal travelling wave-like structures organized on a circular path concentric with the rotation axis of the fluid. Figure 13 compares the observations from dye injection by Vanyo *et al.*^{21} with the vorticity in our numerical simulations. We argue that the spiraling structure observed by Vanyo *et al.*^{21} can be explained as a pair of traveling inertial modes excited through a parametric resonance. In Vanyo’s experiments, which used a spheroid, three mechanisms may generate a parametric instability: the two mechanisms proposed by Kerswell^{27} and the conical shear layer driven one investigated in the present paper. To establish the dominant mechanism, we first calculate the differential rotation ε using Busse’s theory^{9} in the spheroid of Vanyo *et al.*,^{21} and then we use Eq. (18) for the CSI and the derivation of Kerswell^{27} for the two other mechanisms. The results summarized in Table I clearly show that the experimental setup of Vanyo *et al.*^{21} is more prone to CSI rather than the classical topography driven parametric instabilities of Kerswell.^{27} In the original experiment in precessing spheroid by Malkus,^{20} the wave-like instability was also reported (see Figure 2(b) of Ref. Ref. 20). Based on our estimates in Table I, both the topographic effect and the conical shear layers may lead to the parametric instability in the experiment by Malkus.^{20}

. | Parameter . | Malkus . | Vanyo . | Earth . | Moon . |
---|---|---|---|---|---|

Ellipticity | η | 4.0 × 10^{−2} | 1.0 × 10^{−2} | ∼2.5 × 10^{−3} | ∼10^{−5} |

Ekman number | E | 1.0 × 10^{−5} | 8.0 × 10^{−7} | ∼10^{−14} | ∼10^{−12} |

Differential rotation | ε | 3.3 × 10^{−1} | 2.0 × 10^{−2} | ∼1.7 × 10^{−5} | ∼3.0 × 10^{−2} |

Inviscid growth rate | |||||

Elliptical instability | O(ε^{2}η) | ∼4.4 × 10^{−3} | ∼4.0 × 10^{−6} | ∼7.2 × 10^{−13} | ∼9.0 × 10^{−9} |

Shear instability | O(εη) | ∼1.3 × 10^{−2} | ∼2.0 × 10^{−4} | ∼4.3 × 10^{−8} | ∼3 × 10^{−7} |

CSI | O(εE^{1/5}) | ∼3.3 × 10^{−2} | ∼1.2 × 10^{−3} | ∼2.7 × 10^{−8} | ∼1.2 × 10^{−4} |

Viscous decay rate | |||||

Viscous damping | O(E^{1/2}) | ∼3.2 × 10^{−3} | ∼8.9 × 10^{−4} | ∼10^{−7} | ∼10^{−6} |

. | Parameter . | Malkus . | Vanyo . | Earth . | Moon . |
---|---|---|---|---|---|

Ellipticity | η | 4.0 × 10^{−2} | 1.0 × 10^{−2} | ∼2.5 × 10^{−3} | ∼10^{−5} |

Ekman number | E | 1.0 × 10^{−5} | 8.0 × 10^{−7} | ∼10^{−14} | ∼10^{−12} |

Differential rotation | ε | 3.3 × 10^{−1} | 2.0 × 10^{−2} | ∼1.7 × 10^{−5} | ∼3.0 × 10^{−2} |

Inviscid growth rate | |||||

Elliptical instability | O(ε^{2}η) | ∼4.4 × 10^{−3} | ∼4.0 × 10^{−6} | ∼7.2 × 10^{−13} | ∼9.0 × 10^{−9} |

Shear instability | O(εη) | ∼1.3 × 10^{−2} | ∼2.0 × 10^{−4} | ∼4.3 × 10^{−8} | ∼3 × 10^{−7} |

CSI | O(εE^{1/5}) | ∼3.3 × 10^{−2} | ∼1.2 × 10^{−3} | ∼2.7 × 10^{−8} | ∼1.2 × 10^{−4} |

Viscous decay rate | |||||

Viscous damping | O(E^{1/2}) | ∼3.2 × 10^{−3} | ∼8.9 × 10^{−4} | ∼10^{−7} | ∼10^{−6} |

Finally, we shall discuss our results in a planetary context. Using a hydrostatic model for the shape of the Lunar’s Core-Mantle boundary (CMB), the differential rotation ε due to the 18.6 yr precession has been estimated to be of order 3% of the mantle rotation.^{3,11} In line with this conclusion, Williams *et al.*^{3} argued that the resulting flow should be turbulent in order to account for the large dissipation inferred from the Lunar Laser Ranging (LLR) measurements. Yet, if one considers the two classical precessional instabilities of Kerswell^{27} due to the polar ellipticity, the inviscid growth rate remains smaller than the viscous decay rate leading to a stable flow. In contrast, considering the effects of the conical shear layers, we obtain an inviscid growth rate 100 times larger than the viscous decay rate assuming that the Ekman number is around 10^{−12} (Table I). These highly super critical conditions are likely to be associated with more complex flows than the one reported here at moderate Ekman number simulations, leading to significant dissipation as observed through inversion of the LLR time series.

In contrast with the moon, the earth’s polar flattening (*η* ≈ 1/400 on the CMB^{21}) and the large period of precession (26 000 yr) result in a liquid core strongly coupled to the mantle. As a consequence, the differential rotation is small at present day, ε ≈ 1.7 × 10^{−5}. As seen in Table I, all three mechanisms may lead to a stable precession-driven flow in the outer core assuming the Ekman number of the outer core is around 10^{−14}. In addition, it has been shown that another set of conical shear layers emanate from the inner core boundary,^{17,41} with an amplitude of *O*(ε*E*^{1/6}) over a thickness of *O*(*E*^{1/3}). In principle, these conical shear layers could also participate to a CSI. A stability criterion can then be derived using the same heuristic argument presented herein leading to ε = *O*(*E*^{1/3}). Hence, CSI in the earth’s liquid core may occur if the Ekman number of the outer core is below 10^{−15}.

In planetary condition, however, several other parameters have to be considered to draw a more reliable conclusion such as the effect of stratification and magnetic fields, as well as the interaction of precession driven flows with other sources of motions such as thermo-chemical convection.

## Acknowledgments

We would like to acknowledge Andrew Jackson for very useful discussions and suggestions on this study. We thank Andreas Fichtner for computational assistance. Simulations were run on Swiss National Supercomputing Center (CSCS) under the Project No. s369. Figure 13(b) is reproduced from the Figure. 5(b) in Ref. 21, for which we thank the Oxford University Press and the original authors for permission to reprint it. Y.L. and J.N. are supported by ERC Grant No. 247303 (MFECE) at ETH Zurich. P.M. is supported by the National Science Foundation EAR CSEDI Grant No. #1067944.

### APPENDIX A: SYMMETRY OF THE BASE FLOW

The linear viscous solution of the governing Eqs. (1) and (2) is symmetric around the origin, namely, ** u**(

**) = −**

*r***(−**

*u***). Neglecting the nonlinear term in Eq. (1) and changing the sign of the position vector, we have**

*r*and

and

where

and

are symmetric and antisymmetric parts, respectively. We can see that only the symmetric solution is forced by precession.

### APPENDIX B: BASE FLOW IN THE MANTLE FRAME AND PRECESSION FRAME

In order to provide different views of precession driven base flow as suggested by one referee, we show the velocities in the mantle frame and the precession frame in Fig. 14 for *P _{o}* = − 1.0 × 10

^{−4}and

*E*= 1.0 × 10

^{−6}. In the mantle frame, the base flow is a combination of the Poincaré mode, i.e., a solid body rotation about an axis in the equatorial plane, and the secondary flow due to the viscous correction. The secondary flow is dominated by the conical shear layers spawned from the critical latitudes, which is hidden behind the solid body rotation. Note that the Poincaré mode is a travelling inertial mode in the mantle frame. In the view of the precession frame, we see both global rotation of the container and the solid body rotation around an axis in the equatorial plane (the Poincaré mode), leading to a tilted rotation axis of the fluid with respect to the rotation axis of the container. The flow is steady in the precession frame. Again, the secondary flow due to the viscous correction is hidden beneath the solid body rotation.

### APPENDIX C: INERTIAL MODES IN A SPHERE

Inertial modes are solutions of the so-called Poincaré equation which can be written as^{42}

in the cylindrical coordinates (*s*, *ϕ*, *z*) with the no-penetration boundary condition

on the spherical surface *s*^{2} + *z*^{2} = 1. Solutions can be found using separation of variables and they are

where *γ* = 0 if (*l* − *m*) is even and *γ* = 1 if (*l* − *m*) is odd, and *x _{j}* are the $N= 1 2 ( l \u2212 m \u2212 \gamma ) $ zeros of the Legendre polynomial $ P l m ( x ) $ of degree

*l*and order

*m*. The eigen-frequency

*ω*is the

_{klm}*k*th solution of the transcendental equation

For each given integer of *l* and *m*, there are *l* − *m* solutions if *m* ≠ 0 and *l* − 1 solutions if *m* = 0. Each eigen-frequency *ω _{klm}* and the associated eigenfunction are identified by three indexes (

*k*,

*l*,

*m*). Here, we only consider

*m*⩾ 0 and the solution can be extended to

*m*< 0 by the relationship

*ω*(

*k*,

*l*, −

*m*) = −

*ω*(

*k*,

*l*,

*m*).

Finally, the velocity field of the each eigenmode is obtained from

The explicit expression of velocities are given by Zhang *et al.*^{49}