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

David Cébron, Raphaël Laguerre, Jerome Noir, Nathanaël Schaeffer

<figure><img src="image_1.png" /></figure>

## To cite this version:

David Cébron, Raphaël Laguerre, Jerome Noir, Nathanaël Schaeffer. Precessing spherical shells: flows, dissipation, dynamo and the lunar core. Geophysical Journal International, 2019, 219 (Supplement_1), pp.S34-S57. 10.1093/gji/ggz037. insu-01874075v2

HAL Id: insu-01874075

https://insu.hal.science/insu-01874075v2

Submitted on 5 Feb 2019

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
Precessing spherical shells: flows, dissipation, dynamo and
the lunar core

D. Cébron¹*, R. Laguerre², J. Noir³, N. Schaeffer¹†

¹ Université Grenoble Alpes, CNRS, ISTerre, Grenoble, France

² Royal Observatory of Belgium, Brussels, Belgium

³ Institute of Geophysics, ETH Zurich, Zurich, Switzerland

January 17, 2019

Abstract

Precession of planets or moons affects internal liquid layers by driving flows, instabilities and possibly dynamos. The energy dissipated by these phenomena can influence orbital pa- rameters such as the planet’s spin rate. However, there is no systematic study of these flows in the spherical shell geometry relevant for planets, and the lack of scaling law prevents convinc- ing extrapolation to celestial bodies.

We have run more than 900 simulations of fluid spherical shells affected by precession, to systematically study basic flows, instabilities, turbulence, and magnetic field generation. We observe no significant effects of the inner core on the onset of the instabilities. We obtain an analytical estimate of the viscous dissipation, mostly due to boundary layer friction in our simulations. We propose theoretical onsets for hydrodynamic instabilities, and document the intensity of turbulent fluctuations.

We extend previous precession dynamo studies towards lower viscosities, at the limits of
today's computers. In the low viscosity regime, precession dynamos rely on the presence of
large-scale vortices, and the surface magnetic fields are dominated by small scales. Interest-
ingly, intermittent and self-killing dynamos are observed. Our results suggest that large-scale
planetary magnetic fields are unlikely to be produced by a precession-driven dynamo in a
spherical core. But this question remains open as planetary cores are not exactly spherical,
and thus the coupling between the fluid and the boundary does not vanish in the relevant limit
of small viscosity. Moreover, the fully turbulent dissipation regime has not yet been reached
in simulations.

Our results suggest that the melted lunar core has been in a turbulent state throughout its
history. Furthermore, in the view of recent experimental results, we propose updated formulas
predicting the fluid mean rotation vector and the associated dissipation in both the laminar and
the turbulent regimes.

**1 Introduction**

The origin of the magnetic fields of planets and stars is attributed to the dynamo mechanism. It is commonly
thought that most of the dynamos are powered by compositional and thermal convection in the liquid
part of these objects. Nevertheless, this scenario is sometimes difficult to apply. This is for instance
the case for the early Moon, for which the intensity of the magnetic field generated by convection might
not be sufficient (Stegman et al., 2003) or the Earth, where recent estimates of thermal and electrical

<footnote>* david.cebron@univ-grenoble-alpes.fr</footnote>

<footnote>† nathanael.schaeffer@univ-grenoble-alpes.fr</footnote>
conductivity of liquid iron imply that convection would be far less efficient than previously thought (Pozzo et al., 2012). Mechanical forcings constitute then alternative ways to sustain dynamo action (Le Bars et al., 2015), as shown numerically for libration (Wu & Roberts, 2013), tides (Cébron & Hollerbach, 2014; Vidal et al., 2018) or precession. The present study focuses on precession, which has already been demonstrated numerically to be able to grow a magnetic field in spherical shells (Tilgner, 2005, 2007), full spheres (Lin et al., 2016), cylinders (Nore et al., 2011; Cappanera et al., 2016; Giesecke et al., 2018) and cubes (Goepfert & Tilgner, 2016, 2018). Hence, the possibility of a precession driven dynamo in the liquid core of the Earth (Kerswell, 1996) or the Moon (Dwyer et al., 2011) cannot be excluded. However, current numerical simulations operate at viscosities many orders of magnitude higher than natural dynamos. The present work aims at shedding some light on the consequences of precession in spheres, including dissipation and magnetic field generation. To this end, we make extensive use of numerical simulations pushing down the viscosity to the limits of current supercomputers.

A rotating solid object is said to precess when its rotation axis itself rotates about a secondary axis that is fixed in an inertial frame of reference. The first theoretical studies of precession considered an inviscid fluid (Hough, 1895; Sloudsky, 1895; Poincaré, 1910). Assuming a uniform vorticity, they obtained a solution for the spheroid, called Poincaré flow, given by the sum of a solid body rotation and a potential flow. However, the Poincaré solution is modified by the existence of boundary layers, and some strong internal shear layers are also created in the bulk of the flow (Stewartson & Roberts, 1963). In 1968, Busse took into account these viscous effects as a correction to the inviscid flow in a spheroid, by considering carefully the Ekman layer and its critical regions (Busse, 1968; Zhang et al., 2010). Based on these works, Cébron et al. (2010) and Noir & Cébron (2013) have proposed models for the flow forced in precessing triaxial ellipsoids. Beyond this correction approach, the complete viscous solution (including the fine description of all viscous layers) has been obtained for the sphere in the two limit cases of a weak (Kida, 2011) and strong (Kida, 2018) precession rates.

When the precession forcing is large enough compared to viscous effects, instabilities can occur in precessing spherical containers, destabilizing the Poincaré flow (e.g. Hollerbach et al. (2013)). First, the Ekman layers can be destabilized (Lorenzani, 2001) through standard Ekman layer instabilities (Lingwood, 1997; Faller, 1991). In this case, the instability remains localized near the boundaries. Second, the whole Poincaré flow can be destabilized, leading to a volume turbulence : this is the precessional instability (Malkus, 1968). It has been argued by Lorenzani (2001), and more recently by Lin et al. (2015), that the conical shears spawned at the critical latitudes can couple non-linearly with pairs of inertial modes, leading to the Conical Shear Instability (CSI) of Lin et al. (2015).

In this work, we will study the influence of the presence of a solid inner core on the precessional insta-
bility. We thus consider a spinning and precessing spherical shell filled with a conducting fluid. Adopting
the approach used for ellipsoids by Noir & Cébron (2013), we obtain an explicit expression of the fluid
rotation vector in the presence of an inner core. We then derive hydrodynamical stability criteria for the
CSI involving conical shear layers spawned by the outer and the inner spherical shell. Finally, based on
our hydrodynamic simulations, we investigate precession driven dynamos in different flow regimes. The
pioneering results obtained by Tilgner (2005) are extended to smaller viscosities, in the hope to reach an
asymptotic regime relevant for planetary fluid layers.

Our paper is organized as follow. Section 2 presents the governing equations of the problem and a brief description of the numerical method used to solve the equations. A reduced model for the base flow in spherical shells is then derived and compared to numerical simulations in Section 3.1, while transition to unstable flows is studied in Section 3.2. Precession driven dynamos are studied in Section 4. Finally, we apply our findings to the Moon (§5) and draw some conclusions (§6).

**2 Description of the problem and mathematical background**

We consider an incompressible Newtonian fluid of density ρ, kinematic viscosity ν, electrical conductivity
γ, and magnetic permeability μ, enclosed in a spherical shell of outer radius R and inner radius Rᵢ. In the
following, the outer boundary is also named CMB, for core-mantle boundary, whereas the inner one is also
called ICB, for inner core boundary. When present, the ICB rotation vector is assumed to be the same than
the CMB one (thus precessing with the mantle). The cavity rotates with an angular velocity Ωs = Ωŝz and
<figure><img src="image_3.png" /><figcaption>Figure 1: Schematic description of the problem.</figcaption></figure>

is in precession at $\Omega_p = \Omega_p \hat{k}$, with the unit vector $\hat{k}$ and $\hat{k} \cdot \hat{z} = \cos \alpha$. We define the frame of precession as the frame of reference precessing at $\Omega_p$, in which we construct a Cartesian coordinate $(\hat{x}, \hat{y}, \hat{z})$ system centered on the sphere, with $\hat{z}$ along $\Omega_s$, $\hat{x}$ such that $\Omega_p$ is in the plane xOz and $\hat{y} = \hat{z} \times \hat{x}$ (Fig. 1).

## 2.1 Mathematical formulation

Defining $\Omega_o = \Omega_s + \Omega_p$, we choose $\Omega_o^{-1}$ as the unit of time such that it remains relevant in both limits of large $\Omega_s$ or large $\Omega_p$ (see also Goepfert & Tilgner, 2016). We choose $R$ and $R\Omega_o\sqrt{\mu\rho}$ as the respective units of length and magnetic field. In the frame of reference precessing at $\Omega_p$, the dimensionless governing equations take the form:

$$ \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} = -\nabla p + E \nabla^2 \mathbf{u} - 2 \frac{Po}{1+Po} \hat{\mathbf{k}} \times \mathbf{u} + (\nabla \times \mathbf{B}) \times \mathbf{B}, \quad (1) $$

$$ \nabla \cdot \mathbf{u} = 0 \qquad (2) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \frac{E}{Pm} \nabla^2 \mathbf{B} + \nabla \times (\mathbf{u} \times \mathbf{B}), \quad (3) $$

$$ \nabla \cdot \mathbf{B} = 0, \qquad (4) $$

where $p$ is the reduced pressure accounting for centrifugal forces (in our simulations, $p$ is eliminated by taking the curl of equation 1). The four dimensionless parameters controlling the dynamics of the system are the Ekman, magnetic Prandtl, Poincaré numbers and aspect ratio, respectively defined by $E = \nu/(\Omega_o R^2)$, $Pm = \mu\gamma\nu$, $Po = \Omega_p/\Omega_s$, and $\eta = R_i/R$.

Note that with our choice of time scale, the instantaneous rotation vector of the spherical container is $\tilde{\Omega}_s = (1 + Po)^{-1}\hat{z}$. Hence, our modified Ekman number is related to the more classical definition based on $\Omega_s$, $E' = \nu/(\Omega_s R^2) = E(1 + Po)$. Note also that, in this study, our convention is to consider positive $Po$ and $\alpha > 90^\circ$ for retrograde precession, while others consider negative values of $Po$ with $\alpha \in [0, 90^\circ]$ (e.g. Dwyer et al., 2011).
<figure><img src="image_4.png" /><figcaption>Figure 2: Parameter space covered by numerical studies and for typical planetary cores (opposite corners show here the current Earth's and Moon's liquid cores).</figcaption></figure>

We also use the spherical coordinate system (r, θ, φ), r being the radial distance, θ the colatitude, and φ the azimuthal angle (with Ox the axis of zero longitude).

## 2.2 Numerical approach

We impose no-slip boundary conditions at both boundaries, i.e. $u = \tilde{\Omega}_s \times r$ in our precessing frame of reference. We consider insulating boundary conditions for the magnetic field at the CMB and an inner core electrical conductivity equal to that of the fluid. The problem is solved using the XSHELLS code (freely available at https://bitbucket.org/nschaeff/xshells). This high performance, parallel Navier-Stokes solver works in spherical coordinates using a toroidal-poloidal decomposition and a pseudo-spectral approach. Note that with this approach, equations (2) and (4) are automatically satisfied. The spherical harmonic transforms are performed using the efficient SHTns library (Schaeffer, 2013) while the radial direction is discretized with second order finite differences. XSHELLS has been benchmarked on convective dynamo problems with or without a solid inner-core (Marti et al., 2014; Matsui et al., 2016). In the following, the so-called hydrodynamic simulations do not take into account the magnetic field (B), whereas our so-called magnetohydrodynamic (MHD) simulations solve the full system (1)-(4). The range of hydrodynamic parameters investigated in this study are summarized in Fig. 2 and compared to the typical values expected for planetary cores. In addition we explore the range $0.2 < Pm < 20$ in our MHD simulations.

We checked on a few cases that computing in the mantle frame (i.e. the frame co-rotating with the solid shell) gives the same result as in the precession frame. Note that in both frames, the flow is dominated by a strong solid-body rotation (along a third axis). Such large advection speeds are known to cause accuracy issues that can lead to suppression of instabilities (e.g. Springel, 2010). The time-step in XSHELLS is adjusted to ensure stability, but we checked accuracy by dividing the time-step by 3 to 7 on several cases. We found that when instabilities are saturated, stable time steps are small enough to ensure accurate results. However, we have noticed that for one specific case very close to the onset of instability, the growth rate was biased towards stability. As not all runs could be checked, it is not impossible that a few such cases are still included in our results, but it would not change the conclusion drawn in this paper.
In the following, we call hydrodynamic cases those with no magnetic field or with magnetic energy lower than $10^{-16}$ to ensure no perturbation of the flow by the magnetic field.

# 3 Hydrodynamics

## 3.1 Laminar Base Flow

The primary flow forced by precession in a sphere is mainly a tilted solid body rotation, a flow of uniform vorticity (Poincaré, 1910). In a spherical container, the direction and amplitude of the fluid rotation vector are governed by a balance between the viscous torque at the core-mantle boundary and the gyroscopic torque resulting from the precession of the liquid core (Busse, 1968; Noir et al., 2003). We briefly recall in section A.1 of appendix A the derivation of Busse (1968). A limitation of this general formulation, accounting for the Ekman boundary layer action, arises from the implicit nature of the final equation. Indeed, while approximate expressions can be obtained in certain limits (e.g. Boisson et al. (2012)), we cannot derive a general analytical explicit solution. In the context of a precessing ellipsoid, Noir & Cébron (2013) proposed an alternative to the torque balance, using a simpler ad-hoc viscous term. Using this successful approach in a spherical shell, we obtain an explicit expression of the dimensionless fluid rotation vector $\Omega$ in the frame of precession (see section A.2 of appendix A for details):

$$ \Omega_x = \frac{1}{1 + Po} \frac{[\lambda_i + \chi \cos \alpha] \chi \sin \alpha}{\chi(\chi + 2\lambda_i \cos \alpha) + |\lambda|^2}, \quad (5) $$

$$ \Omega_y = -\frac{1}{1 + Po} \frac{\chi \lambda_r \sin(\alpha)}{\chi(\chi + 2\lambda_i \cos \alpha) + |\lambda|^2}, \quad (6) $$

$$ \Omega_z = \frac{1}{1 + Po} \frac{\chi(\chi \cos^2 \alpha + 2\lambda_i \cos \alpha) + |\lambda|^2}{\chi(\chi + 2\lambda_i \cos \alpha) + |\lambda|^2}, \quad (7) $$

where $\chi = Po/\sqrt{E}$. Note that, in the sphere, there is always a single solution for the basic flow $\Omega$, whereas multiple solutions can be obtained in spheroids or ellipsoids (Noir et al., 2003; Cébron, 2015; Vormann & Hansen, 2018).

In equations (5)-(7), the complex viscous damping coefficient $\lambda = \lambda_r + i\lambda_i$ of the spin-over mode is the sum of the (real) viscous damping $\lambda_r$ and the (real) viscous correction to the inviscid frequency $\lambda_i$ of the spin-over mode (37), such that $|\lambda|^2 = \lambda_r^2 + \lambda_i^2$. For a full sphere, Greenspan et al. (1968) derived an expression for $\lambda$ at the order $O(\sqrt{E})$. To account for the presence of an inner core, one can use

$$ \lambda = \lambda_{(\eta=0)} \frac{1 + \eta^4}{1 - \eta^5} \quad (8) $$

for a no-slip inner core (Rieutord, 2001). Moreover, viscous corrections of $\lambda$ should be considered for finite values of $E$ (see appendix A.1). In the following, we account for these various corrections by calculating $\lambda$ using equation (38).

Following Lin et al. (2015), we note $\epsilon$ the differential rotation $\epsilon = ||\Omega - \tilde{\Omega}_s||$ between the fluid and the cavity. From the approximated solution of uniform vorticity (5)-(7) we derive the following analytical expression

$$ \epsilon = \left| \frac{\chi \sin \alpha}{1 + Po} \right| \frac{1}{\sqrt{\chi(\chi + 2\lambda_i \cos \alpha) + |\lambda|^2}}. \quad (9) $$

The equations and the forcing being centro-symmetric, we calculate the uniform vorticity from the symmetric toroidal energy $E_{s,T}$ of the flow that superimposes on the rotation of the cavity,

$$ E_{s,T} = \frac{1}{2} \int \boldsymbol{u}_{s,T}^2 \, dV, \quad (10) $$

with $\boldsymbol{u}_{s,T} = [\boldsymbol{u}_{m,T}(\boldsymbol{r}) - \boldsymbol{u}_{m,T}(-\boldsymbol{r})]/2$ and $\boldsymbol{u}_{m,T} = \boldsymbol{u}_T - \boldsymbol{\Omega}_s \times \boldsymbol{r}$, where $\boldsymbol{u}_T$ is the toroidal component of the velocity field such that $\boldsymbol{u}_T = \nabla \times (\boldsymbol{T}\boldsymbol{r})$. Thus, the (uniform vorticity) differential rotation is calculated
<figure><img src="image_5.png" /><figcaption>Figure 3: Comparison between the differential rotation $\epsilon$ of the reduced model (eq. 9 – dashed line) with all 866 simulations with the required diagnostics (including dynamos). Filled and open symbols represent stable and unstable flows, respectively. The inset shows the distribution of relative error between reduced model and simulations for all 866 runs (red) and for the 720 ones not affected by magnetic field (black).</figcaption></figure>

using

$$ \epsilon = \sqrt{\frac{2E_{s,T}}{I_c}}, \qquad (11) $$

where $I_c = 8\pi(1 - \eta^5)/15$ is the moment of inertia (per unit of mass) of the spherical shell enclosing the fluid. Fig. 3 shows that equation (9) is in quantitative agreement with the uniform vorticity component deduced from the energy in our hydrodynamic simulations, which validates both our reduced model and the estimated differential rotation in the numerics from the toroidal symmetric energy. However, relative deviations up to 10% remain between the measured and predicted $\epsilon$ (see error distribution in Fig. 3). The distribution of deviations does not change much when keeping only the lowest viscosity or only the stable simulations. Thus, these deviations are either due to the way we measure the differential rotation, or to the approximation used to obtain equation 9 (see §A.2). Finer measures of the differential rotation may reduce these deviations, but were unfortunately not available from our runs. While the origin of these deviations is thus unclear, differential rotation of most cases are accurately predicted by equation 9 and measured using the toroidal symmetric energy by equation 11.

In addition to the uniform vorticity flow, a secondary viscous circulation will develop in the interior due to the Ekman pumping at the ICB and at the CMB. In contrast with the classical uniform thickness of the planar Ekman boundary layer, oscillatory motions in a sphere result in local discontinuities at some critical latitudes, propagating in the interior along cones aligned with the axis of rotation (Fig. 4). For an oscillatory motion at an angular frequency $\omega$, it corresponds to the colatitude $\alpha_c$ such that $\cos \alpha_c = \omega/2$. For precession, $\omega = 1$, which gives $\alpha_c = 60^\circ$, i.e. a critical latitude of 30° (Kerswell, 1995; Hollerbach & Kerswell, 1995; Stewartson & Roberts, 1963; Noir et al., 2001). Both the ICB and the CMB generate oblique shear layers. Based on the work of Stewartson & Roberts (1963), Noir et al. (2001) corrected the predictions of Kerswell (1995) and obtained the following scalings for the width and strength of the flow
<figure><img src="image_6.png" /><figcaption>Figure 4: Schematic of the conical shear layers spawn from the critical colatitude $\alpha_c$ given by $\cos \alpha_c = 1/2$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 5: Contour plot of an instantaneous azimuthal average of the kinetic energy in the fluid frame for (a) $\eta = 0.01$, (b) $\eta = 0.3$, (c) $\eta = 0.5$ and (d) $\eta = 0.7$. In each case the flow is stable. $\alpha = 120^\circ$, $Po = 1 \times 10^{-3}$ and $E = 3.0 \times 10^{-5}$. Contours range from 0 (white) to 0.028 (dark red).</figcaption></figure>
in these oblique shear layers:

$$ \delta_{CMB} \propto E^{1/5}, \quad U_{CMB} \propto \epsilon E^{1/5}. \qquad (12) $$

and,

$$ \delta_{ICB} \propto E^{1/3}, \quad U_{ICB} \propto \eta \epsilon E^{1/6}, \qquad (13) $$

The subscript refers to the region at the origin of the internal structures (Fig. 4).

To identify these structures in our numerical simulations, we look at the flow from a frame of reference attached to the mean rotation of the fluid. Fig. 5 represents the azimuthal average kinetic energy in a meridian plan for different inner core sizes, clearly exhibiting conical shear layers. The ones spawned from the ICB critical latitude being more intense, they quickly dominate as the inner core radius increases.

## 3.2 Hydrodynamic instabilities

We track the instability onset by looking for non-zero anti-symmetric energy (Lorenzani, 2001; Lin et al., 2015):

$$ E_a = \frac{1}{2} \int u_a^2 dV, \qquad (14) $$

where $u_a = [u(r) + u(-r)]/2$. Although centro-symmetric unstable flows exist (Hollerbach et al., 2013), they are however limited to a narrow range of parameters with moderate to large values of $E$ (Lin et al., 2015). Disregarding these possible instabilities, we focus on unstable flows which break the centro-symmetry, i.e. with $E_a \neq 0$.

It has been recently argued by Lin et al. (2015) and Lorenzani (2001) that the oscillating conical shear layers, originating from the CMB Ekman boundary layers, can couple non-linearly with two inertial modes, $u_1$ and $u_2$, leading to a parametric resonance, the so-called Conical Shear Instability (CSI, Lin et al., 2015). For precession, the two free inertial modes are subject to the following selective rules (Kerswell, 2002):

$$ \begin{aligned} \omega_1 \pm \omega_2 &= 1, \\ m_1 \pm m_2 &= 1, \\ l_1 \pm l_2 &= 1. \end{aligned} \qquad (15) $$

where $\omega_{1,2}$ and $m_{1,2}$ are the frequencies and azimuthal wave numbers of the two free inertial modes, respectively. $l_{1,2}$ is the degree of the Legendre polynomial characterizing the latitudinal complexity.

Based on the scaling of the oblique shear layers emanating from the CMB (Fig. 4), Lin et al. (2015) proposed that the onset of the CSI is governed by a critical value of the differential rotation $\epsilon$, scaling as

$$ \epsilon_c \propto E^{3/10}, \qquad (16) $$

in agreement with their numerical simulations in a full sphere as well as the experimental results from Goto et al. (2014). The same argument can be used in the spherical shell to derive a criteria for the onset of a CSI driven by the oblique shear layers emanating from the CMB as well as from the ICB.

For clarity, we name CSI-ICB and CSI-CMB, the parametric instabilities of the conical shear layers spawned from the Ekman boundary layer of the ICB and CMB, respectively. Adopting the same approach as Lin et al. (2015), with the scaling of shear layers shown in Fig. 4, we obtain

$$ \epsilon_c = K_{CMB} E^{3/10}, \qquad (17) $$

for the CSI-CMB and

$$ \epsilon_c = K_{ICB} \frac{E^{1/3}}{\eta} \qquad (18) $$

for the CSI-ICB, where $K_{CMB}$, $K_{ICB}$ are two constants.

In addition, the CMB and ICB Ekman boundary layers may be unstable to a local shear instability. The onset of this boundary-layer instability is characterized by the local Reynolds number $Re_{bl} = v\delta/\nu \approx 55$
<figure><img src="image_8.png" /><figcaption>Figure 6: Stability criteria, where filled and open symbols represent stable and unstable flows, respectively. Top: for CMB related instabilities, the dashed and solid lines represent respectively the theoretical CSI (eq. 17) and boundary layer (eq. 19) stability criteria. Bottom: same criteria for the ICB related instabilities (eq. 18 and 20). We have set $K_{CMB} = K_{ICB} = 8$ and $K_{BL} = 55$. Anomalous points (ie stable points above the lines and unstable points below) are plotted on top. Only the 473 points with $Po < 0.3$ and low magnetic energy are shown. Note that keeping only $Po < 0.1$ removes the two anomalous stable points near $E = 3 \times 10^{-5}$, which have been carefully checked to be stable even when the time-step is reduced (see §2.2).</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 7: Two ways to collapse the turbulent fluctuation velocity $U_{rms} = \sqrt{4E_a/V}$ of the 84 simulations with $Po < 0.1$, $E \le 10^{-4}$ and low magnetic energy. The dashed lines are tentative and speculative laws that explain some of our data; they are not supported by theory and should not be used for extrapolation.</figcaption></figure>
(e.g. Lorenzani, 2001; Sous et al., 2013) based on the Ekman layer thickness $\delta = \sqrt{\nu/\Omega} = \sqrt{ER}$ and the maximum (differential) tangential velocity $v$ at the edge of the boundary layers, with $v = \epsilon$ at the CMB and $v = \eta\epsilon$ at the ICB. The stability criteria for the ICB and CMB read

$$Re_{CMB} = \frac{\epsilon}{\sqrt{E}} > K_{BL}, \qquad (19)$$

$$Re_{ICB} = \frac{\eta\epsilon}{\sqrt{E}} > K_{BL}, \qquad (20)$$

with $K_{BL} \approx 55$ (Lorenzani, 2001; Sous et al., 2013), and the associated instabilities will be noted respectively BL-CMB and BL-ICB. In all cases, the outer boundary will become unstable first.

We arbitrarily distinguish the stable and unstable cases by $E_a < 10^{-10}$ and $E_a > 10^{-10}$, respectively. To unravel the underlying destabilizing mechanism, we represent our results in an $(\epsilon, E)$ parameter space against the above mentioned onset criterion for the CSI-CMB and CSI-ICB (Fig. 6). We cover a wide range of parameters, $\eta = 0, 0.1, 0.35, 0.7$, $Po < 0.3$ and $E \le 10^{-3}$. Contrarily to the CMB-related criteria, the ICB-related ones do not separate stable from unstable points, and we thus conclude that the first instability is due to the CMB. Fig. 6 further suggests that in our numerical simulations the first instability is a CSI-CMB at moderate to large Ekman numbers. Meanwhile, below $E = 3 \times 10^{-5}$, our theory predicts that the boundary layer is unstable before the CSI-CMB. One can notice the robustness of the CMB-related instability criteria for all inner core radii investigated up to $\eta = 0.7$.

In order to distinguish clearly between the two mechanisms one should carry out numerical simulations in the range $10^{-8} < E < 10^{-6}$, a range of values still hardly accessible. Further considerations regarding the prevalence of a possible viscous boundary layer at the ICB will be discussed in the last section in a geophysical context at very low Ekman numbers.

As a proxy for the turbulent fluctuation level around the mean flow, let us define the root-mean-square fluctuation velocity as $U_{rms} = \sqrt{4E_a/V}$, where $V = 4/3\pi(1 - \eta^3)$ is the fluid volume and the factor 4 assumes equipartition between symmetric and anti-symmetric turbulent fluctuations. While the onset of instability is reasonably well captured by the BL-CMB at low viscosity ($E \lesssim 3 \times 10^{-5}$), it is more difficult to understand the amplitude of the saturated turbulent velocity $U_{rms}$. Figure 7 shows two ways we found to collapse our data, including a viscosity-free law (Fig. 7b). None of them is fully convincing and further calculations at lower Ekman numbers are clearly necessary to uncover a saturation scaling law. Nevertheless, some systematic behavior is captured in this figure with rather low amplitude fluctuations $\epsilon\sqrt{E} \lesssim U_{rms} < \epsilon/10$ in the planetary parameter range $\epsilon \ll 1, E \ll 1$.

## 3.3 Instability flow structure

Fig. 8 represents the axial velocity during the growth phase of the instability in the system for three different inner core sizes ($\eta = 0.01, 0.3, 0.7$) and for the following control parameters, $\alpha = 120^{\circ}$, $E = 3.0 \times 10^{-5}$, $Po = 8 \times 10^{-3}$ except for the largest inner core for which the instability is detected only for $Po = 8.5 \times 10^{-3}$. Not surprisingly, the smallest inner core (Fig. 8a) is comparable to the full sphere case of Lin et al. (2015) with two inertial modes of wave numbers $m = 17$ and $m = 18$ developing in the outer part of the fluid domain. As we increase the volume of the inner core, the modes involved in the parametric resonance remains high order near onset (Fig. 8b and c) and tend to develop in regions above and below the inner core, again exhibiting pairs of inertial modes satisfying the parametric resonant conditions. These results near onset suggest that a CSI-CMB mechanism is operating at this low value of the Ekman number.

Fig. 9 shows the total anti-symmetric kinetic energy and the anti-symmetric kinetic energy for each azimuthal wave numbers $m$ for increasing Poincaré numbers, from $Po = 7 \times 10^{-3}$ just above onset to about 2 times critical at $Po = 1.3 \times 10^{-2}$ for $\eta = 0.01$. Near onset, the parametric instability remains at a saturated state (Fig. 9(a)). As shown by Fig. 9(b), a very small increase of the forcing at $Po = 8.5 \times 10^{-3}$ leads to a quasi periodic behaviour of the system with a typical period of order $T = 500$ (this behaviour is reminiscent of the resonant collapses observed e.g. in Lin et al., 2014). Despite the modest increase in $Po$, we observe an anti-symmetric energy at saturation that is an order of magnitude larger, yet pairs of modes in parametric resonance can be clearly identified supporting a CSI-CMB underlying mechanism.

In Fig. 9(c), we further increase the precession rate to $Po = 1.3 \times 10^{-2}$. We do not observe a clear initial growth of any particular modes but intermittent states with chaotic and quasi steady phases appear.
<figure><img src="image_10.png" /><figcaption>Figure 8: Isosurfaces of the anti-symmetric energy in the mean fluid rotating frame frame for (a) $\eta = 0.01$, $Po = 8 \times 10^{-3}$ (b) $\eta = 0.3$, $Po = 8 \times 10^{-3}$ and (c) $\eta = 0.7$, $Po = 8.5 \times 10^{-3}$. In each case $\alpha = 120^{\circ}$ and $E = 3.0 \times 10^{-5}$. Color correspond to positive (black) and negative (white) axial velocity. The snapshots are taken during the initial growth phase of the instability.</figcaption></figure>

At $Po = 1.3 \times 10^{-2}$ the dynamic of system is no longer quasi periodic in time, rapid fluctuations are observed together with periods of stable energy of modes with $m = 3$ as for instance between $t \sim 1250$ and $t \sim 2100$. The system alternates between phases in which $E_a$ is concentrated in the $m = 3$ azimuthal mode, contrasting with phases during which $E_a$ is distributed over a wider range of $m$. In Fig. 10, we show that three cyclonic large-scale vortices (LSV) are seen in the phases where $E_a$ is concentrated in $m = 3$, while they are absent in the other phases.

Finally, we increase $\eta$ from 0.3 to 0.7 at $E = 3.0 \times 10^{-5}$, $Po = 8.5 \times 10^{-3}$, $Po = 1.3 \times 10^{-2}$ and $Po = 2 \times 10^{-2}$, to investigate the influence of the inner core on the dynamics above the onset. At $Po = 8.5 \times 10^{-3}$ we observe a similar dynamics as for $\eta = 0.01$ with quasi-periods increasing with $\eta$, the mode structures are qualitatively similar to the full sphere as seen on Fig. 8. While at $\eta = 0.1, 0.3$ we could still observe LSV, although with a shorter life time, they completely disappear for $\eta > 0.3$ for $Po = 1.3 \times 10^{-2}$ and $Po = 2 \times 10^{-2}$. At this Poincaré numbers with moderate to large inner cores, the flow exhibits small scale structures with rapid temporal variations. In contrast with the mode-coupling regime, the typical time scale of the energy fluctuations decreases with increasing inner core size. These calculations quickly become computationally challenging as the Poincaré number is increased. Since an extensive survey of the parameter space to characterise the onset of the LSV is beyond the scope of this paper, we did not explore the higher $Po$ range where the LSV may be driven even in the presence of a large inner core.

## 3.4 Energy Dissipation

The total viscous dissipation is given by

$$D_{\nu} = E \int (\nabla \times \mathbf{u}_{m})^{2} dV \qquad (21)$$

with $\mathbf{u}_m = \mathbf{u} - \mathbf{\Omega}_s \times \mathbf{r}$ the velocity field in the mantle frame such that $\mathbf{u}_m = 0$ at the boundaries. This dissipation arises, in absence of instability, purely through viscous friction in the boundary layers at the inner and outer walls and in the oblique shear layers in the bulk.

First, we consider the oblique shear layers in the bulk. To calculate the associated dissipation, one can restrict the volume integral to the shear layers in equation (21). For the conical shear layer originating from the CMB we obtain $D_{\nu} \sim \epsilon^2 E^{6/5}$ using the scalings (12), whereas for the one originating from the ICB we find $D_{\nu} \sim \epsilon^2 \eta^2 E$ using the scalings (13). The viscous dissipation due to boundary friction is given by $D_{\nu}^{BL} = \Gamma_{\nu} \cdot (\Omega - \tilde{\Omega}_s)$ where $\Gamma_{\nu}$ is the associated viscous torque. Having shown in section 3.1 that the reduced model performs equally well than the model of Busse (1968) but provides explicit expression for
<figure><img src="image_11.png" /><figcaption>Figure 9: Time evolution of the anti-symmetric kinetic energy of each azimuthal mode $m$ in the fluid frame (color map) and of the total anti-symmetric kinetic energy (White solid line). Fixed $\eta = 0.01$, $\alpha = 120^\circ$ and $E = 3.0 \times 10^{-5}$. The Poincaré number increases from (a) $Po = 7 \times 10^{-3}$, (b) $Po = 8.5 \times 10^{-3}$ and (c) $Po = 1.3 \times 10^{-2}$, with dashed gray lines corresponding to the times of figure 10.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: Snapshot of an equatorial cross-section of the axial vorticity $\omega_z$ in the fluid frame for $\eta = 0.01$, $Po = 1.3 \times 10^{-2}$, $\alpha = 120^\circ$ and $E = 3.0 \times 10^{-5}$, same as Fig. 9. (a) $t = 822$, (b) $t = 1142$ and (c) $t = 1567$</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11: Total viscous dissipation $D_ν$ obtained in simulations (symbols) and laminar dissipation $D_ν^{lam}$ (solid line) predicted by equation (24), as a function of parameter $x_o \simeq 2.62 E^{1/2}/Po$. Filled and open symbols represent respectively stable and unstable flows for the 573 cases with $Po < 0.3$ and $E \le 10^{-3}$, including MHD runs (magenta symbols, magnetic energy $Em > 10^{-16}$). The inset shows the relative difference between $D_ν$ and $D_ν^{lam}$ (with the MHD runs in black).</figcaption></figure>

$\epsilon$, we use equation (43) for $\Gamma_ν$ to obtain the laminar dissipation

$$D_{\nu}^{BL, lam} = -I_{c} \epsilon^{2} \lambda_{r} \sqrt{E}. \qquad (22)$$

Note that $D_ν^{BL, lam}$ only differs by a factor $\Omega^{1/2}$ (with $\Omega \lesssim 1$, see appendix A) from the one obtained with the viscous torque (42) of the Busse (1968) model. For small Ekman numbers $E \ll 1$, the dissipation in the oblique shear layers is thus negligible, and the total laminar dissipation $D_ν^{lam}$ reduces to $D_ν^{BL, lam}$. Substituting $\epsilon$ with its expression (9) into equation (22) leads to

$$D_{\nu}^{lam} = \frac{8\pi(1 - \eta^5)|\lambda_r|\sqrt{E}}{15|1 + Po|^2} \frac{|\chi \sin \alpha|^2}{|\chi(\chi + 2\lambda_i \cos \alpha) + \lambda_r^2 + \lambda_i^2|} \qquad (23)$$

Defining $x_o = |\lambda|E^{1/2}/Po = |\lambda|/\chi$ and noticing that the contribution of $\lambda_i$ is always negligible, expression (23) can be rewritten as

$$D_{\nu}^{lam} = \frac{2x_o}{1 + x_o^2} D_{max}^{lam}, \qquad (24)$$

with

$$D_{max}^{lam} = \frac{4\pi (1 - \eta^5) |\lambda_r| Po \sin^2 \alpha}{15 (1 + Po)^2 |\lambda|} \qquad (25)$$

$D_{max}^{lam}$ is the maximum laminar non-dimensional viscous dissipation obtained for $x_o \simeq 1$, which is independent of $E$.
<figure><img src="image_14.png" /><figcaption>Figure 12: Total viscous dissipation $D_ν$ in simulations normalized by the laminar dissipation $D_ν^{BL,lam}$ computed from equation (22) for the 172 simulations with $Po < 0.05$ and $E \le 10^{-4}$ and low magnetic energy. Filled and open symbols represent stable and unstable flows, respectively. For $Re_{bl} \gtrsim 100$ we see the transition to the turbulent regime characterized by a larger dissipation.</figcaption></figure>

Fig. 11 presents $D_ν/D_{max}^{lam}$ from our numerical simulations together with the laminar estimate (24). It shows that the laminar boundary-layer dissipation given by the explicit expression (24) captures the main contribution to viscous dissipation in our simulations, with less than 30% error. In particular, the maximum viscous dissipation is reached in our low Ekman and low $Po$ simulations for $x_o \simeq 1$ that is $Po \simeq 2.6 E^{1/2}$. Note that the current Earth would be at $x_o \sim 1$, i.e. near the maximum of dissipation, whereas the current Moon would be far from this maximum, at $x_o \sim 10^{-3}$.

Beyond the laminar boundary dissipation, we can compare our results with estimates of turbulent dissipation. Kerswell (1996) derived an upper bound for viscous dissipation which reads in our notations:

$$D_{\nu}(1 + Po)^3 \le \begin{cases} 8\pi\chi^2/525 & \text{for } \chi = Po/\sqrt{E} \lesssim 3 \\ 0.43 & \text{for } \chi = Po/\sqrt{E} \gtrsim 3 \end{cases} \qquad (26)$$

All our simulations have a dissipation lower than this upper bound, sometimes several order of magnitude lower, even for unstable flows.

More recently, the turbulent boundary-layer dissipation has been investigated by Sous et al. (2013) showing that above a critical value of the boundary layer Reynolds number $Re_{bl} = \epsilon E^{-1/2}$ of order 150, the dissipation becomes significantly larger than the laminar one. Concentrating on simulations with $E \le 10^{-4}$ and $Po < 0.05$, figure 12 shows the ratio of the total measured dissipation $D_ν$ to laminar dissipation $D_ν^{BL,lam}$ given by expression 22, as a function of the local Reynolds number. Stable cases have their viscous dissipation accurately predicted by the laminar boundary friction. However, shortly after the onset of the instability, for $Re_{bl} \gtrsim 100$ we observe an increase of the dissipation of up to 20% compared to the laminar model, in qualitative agreement with Sous et al. (2013), albeit an order of magnitude smaller. To quantitatively test their dissipation law in our setup, one should further reduce the Ekman number while giving particular attention to the measure of the uniform vorticity flow that enters the calculation of $D_ν^{BL,lam}$. However, reaching $Re_{bl} > 400$ will be numerically challenging.
# 4 Precession driven dynamos

We now add a magnetic field and solve the whole system of equations (1)-(4). In all our simulations, the initial magnetic and velocity fields are random. We have produced two databases totaling more than 900 simulations. The first set of more than 750 simulations consists in a broad search in the 4-dimensional parameter space (E, η, Po, Pm) which resulted in only few dynamos. In this set, when an inner-core is present it is conducting with the same conductivity as the fluid. The non-dynamo runs of this set were extensively used in the previous sections. The second set of runs uses an insulating inner-core and is focused on several 1-dimensional paths across the parameter space, yielding a few tens of self-sustained dynamos. The explored parameters are summarized in figure 13, and the full databases are made freely available at https://doi.org/10.6084/m9.figshare.7017137.

## 4.1 Beyond Tilgner (2005)

Tilgner (2005) has shown that precessing spheres can generate dynamos, either driven by the Ekman pumping of the forced basic laminar flow (at relatively large values of E), or driven by the anti-symmetric flow associated with instabilities (at smaller E). More recently, Lin et al. (2016) found that large scale vortices are sometimes generated by these instabilities, and that they contribute to magnetic field generation.

We first focus on the precession rate Po = 0.3 and precession angle $\alpha = 120^{\circ}$ for which the onset for dynamo action has been determined by Tilgner (2005) in his figure 4. Our results are summarized in figure 14, showing that our non-linear dynamos are correctly separated by his critical magnetic Prandtl number $Pm_c$ curve (the solid line in figure 14). Thanks to today's computing facilities and to the highly efficient XSHELLS code, we were able to further decrease the viscosity by a factor 10. Since the forcing (Po = 0.3) is kept constant here, this leads to turbulent flows, which require high resolutions and a high degree of parallelization to simulate. At large and small Ekman numbers, we observe the apparition of two local minima for the critical magnetic Prandtl number $Pm_c$ of the dynamo onset: one for $E \approx 1.5 \times 10^{-3}$ and one for $E \approx 10^{-4}$. If the former can be explained by the transition from base flow driven dynamos to instability driven dynamos, the latter is more difficult to interpret. Indeed, Tilgner (2005) explains the decrease of the critical magnetic Prandtl number $Pm_c$ in the range $10^{-4} \le E \le 10^{-3}$ by assuming that dynamo action takes place above a critical magnetic Reynolds number $Rm_c \propto E_a^{1/2} Pm_c / E$ based on the anti-symmetric energy $E_a$. For $E \le 10^{-4}$, this law is no longer valid. The turbulent fluctuations seem to have a negative effect on dynamo action, both in terms of onset ($Pm_c$ increases when decreasing E) and field intensity (as shown by the diminishing circle area in figure 14). What happens when the viscosity is lowered is illustrated in figures 15 and 16. Although the large scale flow presents similar strength and shape at $E = 1.25 \times 10^{-4}$ and at $E = 1.87 \times 10^{-5}$, small-scale instabilities develop near the outer shell in the latter case. This results in a shedding of the surface magnetic field to small scales. In figure 16, comparing $E = 1.87 \times 10^{-5}$ to $E = 1.25 \times 10^{-4}$, the large scale magnetic field is reduced by a factor 100, while the magnetic energy peaks at scales that will be strongly attenuated with distance, and likely to be undetectable at the planet's surface.

While small-scale turbulent fluctuations might in certain cases sustain a large-scale magnetic field (by a so-called mean-field dynamo, e.g. Moffatt, 1970), our results suggest the opposite in a precessing sphere. In addition, a destructive effect of small-scale fluctuations on the dynamo has also been reported for precessing cubes (Goepfert & Tilgner, 2018) and cylinders (Nore et al., 2014). However, a constructive effect of small-scale fluctuations cannot be excluded in the very low Pm regime relevant for planetary cores but out of reach with direct numerical simulations. Further dedicated studies are needed to address this issue.

However, this path with fixed Po = 0.3 is not appropriate for typical planetary regimes, for which Po ≪ 1, and we now vary the precession rate.

## 4.2 Varying the precession rate

Surprisingly, when setting Po = 0.2 or Po = 0.4 instead of Po = 0.3, the critical magnetic Prandtl number $Pm_c$ required for dynamo action increases (see Fig. 13a). This means that the flow generated by a precession rate Po = 0.2 or Po = 0.4 is less efficient than the one obtained at precession rate Po = 0.3. This already highlights the complicated landscape in which we are trying to find dynamos. Specific values
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13: Saturated dynamos (stars), self-killing dynamos (triangles) and non-dynamo (circles) represented in two parameter-space planes. Symbols outlined in pink are from Lin et al. (2016); a grey outline indicates a conducting inner-core; a blue outline stands for a small insulating inner-core ($\eta = 0.1$) or no inner-core (full sphere). The x-axis is the standard Ekman number based on the gap width.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 14: Successful dynamo (circle) and non-dynamo (cross) simulations for the precession parameters of Tilgner (2005): $\eta = 0.1$, $Po = 0.3$, $\alpha = 120^{\circ}$, with a stress-free insulating inner core. The solid blue line is the critical $Pm$ for dynamo action found by Tilgner (2005), and the dashed blue line is the law $Pm_c = 300E_a^{-1/2}E(1+Po)/(1-\eta)^2$ that he proposed (cast to our definition of $E$). The area of the dots is proportional to the magnetic energy.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 15: Radial velocity in the equatorial plane (top) and radial magnetic field at the core surface (bottom) for two different values of the viscosity, near the onset of dynamo action. Lower viscosity (right: $E = 1.87 \times 10^{-5}$, $Pm = 1$) results in much smaller scales in both velocity and magnetic fields than the larger viscosity (left: $E = 1.25 \times 10^{-4}$, $Pm = 0.75$). Both cases have $Po = 0.3$, $\alpha = 120^{\circ}$ and a small insulating inner-core ($\eta = 0.1$).</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 16: Magnetic and kinetic energy spectra as a function of spherical harmonic degree $ℓ$ at the fluid surface for the magnetic field and below the Ekman layer for the velocity field. The black dashed lines indicate slopes of -3 for the kinetic energy spectra. The two cases have $Po = 0.3$, $α = 120°$, a small insulating inner-core ($η = 0.1$) and $Pm = 1$ and differ by their Ekman number $E = 1.25 × 10^{-4}$ and $E = 1.87 × 10^{-5}$, the latter displaying more energy at smaller scales.</figcaption></figure>

of the precession rate will lead to dynamos whereas neighboring values will not. This is apparent in figure 13, where circles (decaying magnetic field) and stars (dynamos) are entangled. We have not been able to find simple parameter combinations that allowed to disentangle them. For instance, we introduce a precession-based magnetic Reynolds number $Rm^* = |Po \sin \alpha|Pm/E$. Figure 13b shows that a stable dynamo is found for a low $Rm^* = 60$ (at $E = 10^{-4}$, $Po = 0.02$, $Pm = 0.3$), while several cases at $Rm^* \ge 4000$ do not produce a magnetic field. Furthermore, the power-based scaling laws that govern convective dynamos (e.g Christensen et al., 2009; Oruba & Dormy, 2014) do not work here. A possible reason being that the power is injected by viscous coupling and that laminar viscous dissipation at the boundaries remains dominant in the accessible parameter range (see Fig. 12).

## 4.3 Low viscosity dynamos and large-scale vortices

### 4.3.1 Stable dynamos

Decreasing $Po$ and $E$ together, we find a few stable dynamos in full spheres and spherical shells with a small inner-core ($η = 0.1$). The dynamos obtained at the lowest viscosities and forcing (low $Pm$, $E$, $Po$) are all associated with large-scale vortices (LSV, see figure 10 and 20 for examples). The importance of LSV for dynamo action has been already highlighted by Lin et al. (2016), and our study confirms that
<figure><img src="image_20.png" /><figcaption>Figure 17: Time-evolution of the magnetic energy in dynamos at $E = 7.94 \times 10^{-5}$, $Po = 0.02$ and $\alpha = 90^\circ$, with a small insulating inner-core $\eta = 0.1$. While $Pm = 0.3$ is a stationary dynamo, the stronger Lorentz force at $Pm = 0.5$ leads to the loss of the magnetic field. One magnetic diffusion time is $R^2(1 - \eta)^2\mu\gamma$. The link with large-scale cyclones is shown in supplementary animations https://doi.org/10.6084/m9.figshare.7063652.</figcaption></figure>

they play an important role for dynamo action (see also Guervilly et al., 2015, in the context of rotating convection). With stable, persistent LSV, the magnetic energy is rather stable (case $Pm = 0.3$ in figure 17), allowing to obtain dynamos at low viscosity ($Pm < 1$, $E \le 10^{-4}$), seemingly relevant for planetary cores.

As an example, at $E = 7.94 \times 10^{-5}$, $Po = 0.02$, $\alpha = 90^\circ$, we found a stable saturated dynamo at $Pm = 0.3$ (see figure 17). Three stable LSV are seen during the whole simulation, unaffected by the Lorentz force. Furthermore, the fluid rotation vector does not change significantly from the corresponding hydrodynamic case.

### 4.3.2 Self-killing dynamos

However, when increasing the electrical conductivity to $Pm = 0.5$, after the exponential growth of the magnetic field, the Lorentz force becomes strong enough to alter the flow so that the magnetic field decays and never recovers, even after the magnetic field has decayed to very low intensity (case $Pm = 0.5$ in figure 17). While the three LSV are present in the growing phase, they wither away when the magnetic field reaches saturation value.

We found other self-killing precessing dynamos as indicated by the triangles in figure 13. At $E = 10^{-5}$, $Po = 0.005$, $\alpha = 90^\circ$, $\eta = 0$, large-scale vortices are observed together with a growing magnetic field for $0.2 \le Pm \le 1$ until the Lorentz force kills the vortices and the magnetic field immediately decays. This is illustrated in figure 18a. For the limited time we could run these self-killing dynamos, the LSV and hence also the magnetic field have not been able to recover, even though the field has reached levels where the Lorentz force is negligible. Figure 18b shows that the mean rotation axis of the fluid changes slightly but permanently when the magnetic field reaches its peak intensity. We hypothesize that the slight change
<figure><img src="image_21.png" /><figcaption>Figure 18: Top: ratio of magnetic energy $E_b$ over time-averaged antisymmetric kinetic energy $\bar{E}_a$. Bottom: relative variations (in percent) of the projection of the fluid rotation axis on the planet spin axis. Both cases have $E = 10^{-5}$, $\alpha = 90^\circ$, $Pm = 0.3$ and no solid inner-core. They differ only by their precession rate $Po$. The saturation of the magnetic energy at $E_b \approx 0.1 \bar{E}_a$ leads to a change in the rotation axis of the fluid, as hinted by the vertical dashed lines. In the case $Po = 0.005$, this change seems permanent despite the loss of the magnetic field. The link with large-scale cyclones is shown in supplementary animations https://doi.org/10.6084/m9.figshare.7063652.</figcaption></figure>
in rotational state of the fluid is enough to prevent the reformation of the LSV. Furthermore, the system stays around the second rotational state even though the magnetic field has vanished, which suggests a hydrodynamic bistability.

Self-killing dynamos have already been reported in simple laminar dynamo models (Fuchs et al., 1999) or turbulent experiments (Miralles et al., 2015). To our knowledge, it is the first time such self-killing dynamos are reported in a self-consistent, turbulent setup. Hence, it highlights that kinematic precession dynamos (a growing magnetic field without the Lorentz force) do not imply that strong magnetic fields can be sustained once the Lorentz force is taken into account.

### 4.3.3 Intermittent dynamos

In addition, many other dynamos show large fluctuations of their magnetic energy of about a factor 10 to 100, suggesting that the Lorentz force is often pushing the flow to a different attractor before quickly recovering. When the magnetic energy is low, the LSV develop and the magnetic field can grow. When the magnetic energy saturates at a high enough level, the Lorentz force sometimes kills the LSV and thus the magnetic energy decays to a lower level. This behavior is rather common and illustrated by case $Pm = 1$ in figure 17, and case $Po = 0.007$ in figure 18. At the planet's surface, this may appear as an intermittent dynamo, with stronger magnetic field alternating with undetectable magnetic field.

### 4.3.4 Small-scale surface magnetic field at low viscosity

Keeping $E = 10^{-5}$, $\alpha = 90^{\circ}$, $\eta = 0$ and increasing the precession rate from $Po = 0.005$ – a self-killing dynamo – to $Po = 0.007$, the LSV are now able to withstand the Lorentz force, and at $Pm = 0.3$ the dynamo saturates with the magnetic energy fluctuating within a factor 10. This is the stable dynamo we obtained with parameters closest to planetary values, and we double-checked by also computing it in the mantle frame with lower time steps (see §2.2). Time-evolution of the magnetic energy is shown in figure 18a. A snapshot of the corresponding flow and field is shown in figure 19, while the LSV are highlighted in figure 20. Two large-scale cyclones are seen in the bulk, together with small-scale vorticity fluctuations. Near the outer shell, a thick layer, much thicker than the Ekman layer, of intense small-scale vorticity is seen. This leads to a small-scale magnetic field at the surface, while the larger-scale, stronger field does not escape the bulk. This shift towards small-scale surface magnetic field seems robust as the viscosity is decreased toward planetary values.

Note however that, at even lower $Pm$, we cannot exclude the emergence of a large-scale magnetic field produced by the small-scale turbulence (see e.g. Moffatt, 1970).

# 5 Application to the Moon

## 5.1 Time evolution of the lunar precession

Precession has been suggested for driving turbulence and dynamo magnetic fields in the past Moon liquid core (Dwyer et al., 2011). In the light of the present results, we propose to revisit the time evolution of precession driven flows in the lunar core. In the following, we will consider a lunar metallic liquid core of radius $R = 350$ km.

During the Moon history, the variation of the precession angle $\alpha$ can be related to the variation of the semi-major axis $a$ of the lunar orbit by (Eq. 5 of Dwyer et al., 2011)

$$180^{\circ} - \alpha = 0.1075^{\circ}\tilde{a}^{10} - 0.0332^{\circ}\tilde{a}^{9} - 1.0008^{\circ}\tilde{a}^{8} + 0.6110^{\circ}\tilde{a}^{7} + 2.7016^{\circ}\tilde{a}^{6} - 1.7281^{\circ}\tilde{a}^{5} \\ -2.3280^{\circ}\tilde{a}^{4} - 1.4509^{\circ}\tilde{a}^{3} + 6.9951^{\circ}\tilde{a}^{2} - 6.6208^{\circ}\tilde{a} + 5.5828^{\circ}, \quad (27)$$

where $\alpha$ is given in degrees, $\tilde{a} = (a/R_E - 46.6308)/7.7288$ with the Earth radius $R_E$ (this formula gives $\alpha$ for $34.2 \le a/R_E \le 60.2$). Then, $a$ can be related to time using the so-called nominal model of Dwyer et al. (2011), shown in their figure S2 and reproduced in Fig 21(a).

Assuming that the Moon remains synchronized during its history, the lunar spin rate $\Omega_s$ is given by its orbital rate. Thus, using the Kepler law $\Omega_s \propto a^{-3/2}$, one can calculate $E$. Then, we obtain $Po$ by
<figure><img src="image_22.png" /><figcaption>Figure 19: Snapshots (corotating with the fluid) at $E = 10^{-5}$, $Pm = 0.3$, $Po = 0.007$, $\alpha = 90^{\circ}$ and no solid inner-core. The radial magnetic field is shown on the surface. In the bulk: the axial vorticity along fluid rotation axis (top) and the magnetic intensity (bottom – in logarithmic scale).</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 20: Vorticity along the fluid rotation axis in the frame rotating with the fluid, averaged along the fluid rotation axis, for $E = 10^{-5}$, $Po = 0.007$, $\alpha = 90^{\circ}$, $Pm = 0.3$. Turbulent boundary layers have been excluded from the average. Two intense cyclones (red) can be seen around a large central anti-cyclonic (blue) region.</figcaption></figure>

extracting the lunar precession rate $\Omega_p$ from the Fig. 19 of Touma & Wisdom (1994). The time evolution of $Po$ and $E$ over the lunar history are presented in Fig. 21(b).

## 5.2 Flow stability at the lunar CMB during its history

One can calculate the stability of the lunar liquid core for the parametric instability (CSI-CMB) and the boundary layer instability (BL-CMB). To do so, we define a general parameter $\zeta$ as an estimate for the onset distance, given by $\zeta = \epsilon/(K_{CMB}E^{3/10})$ for the CSI-CMB and $\zeta = \epsilon/(K_{BL}E^{1/2})$ for the BL-CMB. An instability is thus expected in both cases when $\zeta > 1$. The results are shown in Fig. 22. It confirms the fact that a CSI-CMB can be currently expected in the Moon, as already proposed by Lin et al. (2015). Beyond this confirmation, this figure furthermore shows that both instabilities are clearly expected during the whole lunar history, with a BL-CMB significantly more unstable.

## 5.3 Turbulent torque and dissipation in the lunar liquid core

According to Fig. 22, the local Reynolds number $Re = \epsilon/E^{1/2}$ at the lunar CMB has decreased during the lunar evolution, from $Re = 4.10^5$, 4 Ga ago, to its current value of $Re = 10^4$. These values are well in the regime, in which Sous et al. (2013) observe turbulent Ekman layers. A turbulent friction can thus be expected at the lunar CMB, as previously suggested for the current lunar core (Williams et al., 2001). Naturally, a different friction would lead to a different differential rotation strength $\epsilon$, and we should thus calculate $\epsilon$ in a self-consistent way in presence of a modified (turbulent) friction. To do so, we simply replace the viscous (laminar) term $\mathcal{L}\Gamma_{\nu}$ in equations (39)-(41) by the following turbulent damping term:

$$ \mathcal{L}\Gamma_{\nu} = \lambda_t ||\Omega - \tilde{\Omega}_s|| (\Omega - \tilde{\Omega}_s), \qquad (28) $$

with a coefficient $\lambda_t$. The dimensionless dissipation can thus generally be written as

$$ D_{\nu} = -I_{c} \lambda_{t} \epsilon^{3}, \qquad (29) $$
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

Figure 21: Time evolution of lunar values. (a) $\alpha$ (solid line), $a$ (dashed line). Note that Dwyer et al. (2011) consider negative values of $Po$ with $\alpha \in [0, 90^{\circ}]$ whereas, here, our convention is to consider positive $Po$ and $\alpha > 90^{\circ}$ for retrograde precession. (b) $Po$ (solid line) and $E$ (dashed line), where $E$ is calculated for the liquid core (with $R = 350$ km, $\nu = 10^{-6} \text{ m}^2\cdot\text{s}^{-1}$).
<figure><img src="image_26.png" /><figcaption>Figure 22: Stability of the lunar liquid core considering the BL-CMB (solid line), with $K_{BL} = 60$, and the CSI-CMB (dashed line), with $K_{CMB} = 8$. Each instability is expected for $\zeta > 1$. Since $\zeta = \epsilon/(K_{BL}E^{1/2})$ for the BL-CMB, the local Reynolds number $Re = \epsilon/E^{1/2}$ at the lunar CMB is also given by the solid line, with a factor $K_{BL} = 60$.</figcaption></figure>

where $\lambda_t$ remains to be obtained. Obtaining expressions of $\lambda_t$ requires to describe the (local) turbulent stress associated to the shear velocity $v_{sh} = (\Omega - \tilde{\Omega}_s) \times r$ generated by the differential rotation $\Omega - \tilde{\Omega}_s$ at the inner and outer boundaries. Noting $\tau_v$ the (local) surface stress per unit of mass, one can write $\tau_v = -\kappa|v_{sh}|v_{sh}$, where $\kappa$ can be seen as a (local) drag coefficient. To close the equations, we need to specify $\kappa$, where the physics of the friction coupling is hidden.

Focusing first on laminar flows, the model of Sous et al. (2013) predicts such flows for $Re \lesssim 150$ and prescribes $\kappa = \xi\epsilon\sqrt{E}/v_{sh}$, where $\xi$ is a constant of order unity ($\xi = 1$ in Sous et al. (2013)). One can then calculate the associated viscous torque $\Gamma_v = \int_S r \times \tau_v dS$ on the surface $S$ of the fluid boundary. Noting the colatitude $\vartheta$, we have $|r \times \tau_v| \propto r^3 \sin^3 \vartheta$, which gives

$$ \lambda_t = -5 \xi \frac{1 + \eta^4 \vartheta}{1 - \eta^5} \frac{E^{1/2}}{\epsilon}, \qquad (30) $$

with $\vartheta = 1$ (resp. $\vartheta = 0$) for a no-slip (resp. stress-free) inner boundary. For $Ek \ll 1$ and no-slip boundaries, our equation (22) is exactly recovered with $\xi = 2.62/5 \approx 0.52$, including the correct dependency in $\eta$ (see equation 8). In the model of Sous et al. (2013), using this value of $\xi$ allows thus to switch naturally from the validated laminar dissipation (22) to a turbulent dissipation.

In the turbulent regime, it is usually assumed that $\kappa$ does not vary in space. Under this hypothesis, the associated viscous torque is then given by

$$ \Gamma_v = \frac{3\pi^2}{4} (1 + \eta^5) \kappa |v_{sh}| v_{sh}, \qquad (31) $$

which recovers equation (55) of Williams et al. (2001), obtained in the particular case $\eta = 0$. Using equation (28), equation (31) leads to

$$ \lambda_t = -\frac{45\pi}{32} \frac{1 + \eta^5 \vartheta}{1 - \eta^5} \kappa, \qquad (32) $$

in the turbulent regime (with $\vartheta$ defined as above). Note the different dependency in $\eta$ compared to equation (30).
<figure><img src="image_27.png" /><figcaption>Figure 23: Dimensionless dissipation of the model based on Sous et al. (2013), for an (arbitrary) illustrating case, in presence of a smooth CMB (dashed red line) or a rough CMB with $z_r = 10^{-5}$ (solid black line), using the parameters $Po = 0.02$, $\alpha = \pi/3$, $\eta = 0$, $(A, B) = (3.3, 3)$, $\xi = 0.52$. We also show, as a horizontal dash-dotted green line, the turbulent model of Yoder (1981), and, as a dotted blue line, the laminar reduced model (22) given by equation (9), using $\lambda_r = -2.62$, $\lambda_i = 0$.</figcaption></figure>

To close the equations in the turbulent regime, a value has to be chosen for $\kappa$ (assumed to be uniform in space). As a first approach, Yoder (1981) simply considered a constant $\kappa = 0.002$ based on Bowden (1953). Later, refined models based on turbulent non-rotating boundary-layer theory have been proposed (Yoder, 1995; Williams et al., 2001). By contrast, for rotating turbulent flows ($Re > 150$), Sous et al. (2013) proposed a self-consistent approach where $\kappa = C_d \cos \alpha_0$ depends on $\epsilon$ through

$$C_d = \frac{u_{0\star}^2}{\epsilon^2} = \frac{k^2}{\left(\ln\left|\frac{u_{0\star}}{2z_0}\right| - A\right)^2 + B^2}, \quad (33)$$

and $\sin \alpha_0 = (Bu_{0\star})/(k\epsilon)$. Here, $k = 0.4$ is the von Karman constant, $(A, B)$ are constants obtained from measurements, $u_{0\star}$ is the unknown friction velocity and $\alpha_0$ is the so-called cross-isobar angle from the geostrophic flow due to the tilt of the velocity vector in the boundary layer. From atmospheric measurements, we typically have $(A, B) = (1.3, 4.4)$, whereas laboratory experiments rather give $(A, B) = (3.3, 3)$ (Sous et al., 2013). Noting $z_r$ the dimensionless root mean square roughness height of the boundary, $z_0 = 0.11E/u_{0\star}$ for a smooth CMB (i.e. $z_r u_{0\star}/E < 60$), whereas $z_0 = z_r/30$ for a rough CMB (i.e. $z_r u_{0\star}/E > 60$). From equation (32) and (33) we obtain $\lambda_t$, which we substitute in (31) to obtain the turbulent viscous torque. We finally self-consistently solve for the torque balance (39)-(41). This model is implemented in the updated FLIPPER program (initially introduced in Cébron, 2015), a MATLAB script calculating the theoretical uniform vorticity flow in precessing ellipsoids (https://www.mathworks.com/matlabcentral/fileexchange/50612-flipper).

Let's consider a generic case to illustrate a possible scenario for the viscous dissipation as we lower
<figure><img src="image_28.png" /><figcaption>Figure 24: Dissipated (dimensional) power in the lunar core given by models, compared to the current dissipation (horizontal dash-dotted line) of ~ 76 MW. Laminar model as a dashed line, model of Yoder (1981) as a dotted line, and the model of Sous et al. (2013) as solid lines (smooth CMB and $z_r = 10^{-5}$ for the lowermost and the uppermost ones at -1 Ga).</figcaption></figure>

the Ekman number from numerically accessible values to planetary settings. In this example show in Fig. 23, (Po, α, η) are fixed to arbitrary values and the transition from smooth to rough boundary, i.e. $z_r u_{0*}/E \sim 60$, is taken to arise at $E \simeq 10^{-8}$. At moderate Ekman number $E \gtrsim 10^{-5}$ the flows remain weakly non-linear such that the dissipation is dominated by the laminar processes in the boundary layer leading to $D_\nu \propto I_c \epsilon^2 \sqrt{E}$. At low enough Ekman number $E \lesssim 10^{-5}$, the boundary layer becomes turbulent for the precession parameters considered in this example, and the dissipation is then estimated using the friction model derived from Sous et al. (2013). Two regimes must be distinguished, in the range $10^{-8} \lesssim E \lesssim 10^{-5}$ the roughness is buried in the boundary layer and the dissipation remains weakly dependent on the Ekman number, decreasing with E. When the boundary layer thickness becomes small compared to the roughness of the boundary, $z_r$ becomes the relevant length scale for the dissipation, leading to $D_\nu \propto I_c \epsilon^3$, independent of E, as proposed by Yoder (1981) for the lunar core.

Considering now the lunar core, we take into account a polar flattening of $2.5 \times 10^{-5}$ by solving the torque balance (39)-(41) with the lunar parameters. The dimensional dissipation are shown in Fig. 24 for the different models of viscous torques discussed above. The lunar fluid dissipation currently observed in the Lunar Laser Ranging (LLR) data is 60 MW according to Williams et al. (2001). Since the dissipation is proportional to the so-called fluid core coupling parameter K/C (see e.g. Williams et al., 2001), this value can be updated to 88 MW and 76 MW using the more recent values of K/C respectively given in Williams et al. (2014) and Williams & Boggs (2015). Recovering that a laminar dissipation is not consistent with the currently observed dissipation (Williams et al., 2001), we also show here that the dissipation is expected to be turbulent during the whole lunar history. Note that the turbulent upper bound (26) of Kerswell (1996) gives much larger dissipation ($> 10^{14}$ W in Fig. 24), which shows that our turbulent model based on Sous et al. (2013) is consistent with this theoretical upper bound. Note finally that taking an inner core into account, even with $\eta = 0.7$ (Weber et al., 2011), does not strongly modify the dissipation (modification by a less than 56% for a smooth CMB).
# 6 Conclusion

In this study, we have characterized the precession driven instabilities necessary to sustain magnetic fields in planetary cores. At larger viscosities, the boundary layers remain stable, while parametric resonances (CSI) lead to bulk instabilities. At lower viscosities, the Ekman boundary layer becomes highly turbulent while the CSI is still stirring the bulk. Varying the inner-core size, we characterize both the evolution of the onset and the dissipation with the shell aspect ratio $\eta$. We find that it influences only weakly the onset, which is compatible with experimental findings in a librating ellipsoid (Lemasquerier et al., 2017). While our simulations are still dominated by laminar dissipation, the dissipation in planetary cores may be governed by turbulence in the boundary layers. Based on our numerical results and the experimental work of Sous et al. (2013), we have derived a self-consistent model of the dissipation in precession spherical shells, including both turbulent friction, inner-core and rotation effects. Extrapolating our model to the lunar core, we predict dissipation compatible with the observed LLR data, with little sensitivity to the size of an inner-core.

Adding the magnetic field, we examine the dynamo action over a wide range of parameters towards the planetary regime ($E, Pm, Po \ll 1$). At the lowest investigated viscosity $E = 10^{-5}$ we found a self-sustained magnetic field at $Pm = 0.3$ and $Po = 0.007$. Besides the laminar dynamos at high viscosity which are not relevant for planets, three types of dynamo behaviours emerge at low viscosity ($Pm < 1, E \le 10^{-4}$): (i) stable dynamos, where the magnetic energy reaches a statistically steady state with low to moderate fluctuations; (ii) intermittent dynamos, where the system oscillates between two states with different mean magnetic energies corresponding to two slightly different directions of the fluid rotation axis; (iii) self-killing dynamos, where the magnetic field grows exponentially, saturates but finally decays. At low viscosity, two or three large-scale cyclonic vortices (LSV) are observed during the initial exponential growth of the magnetic field in all three cases. We suspect that far from their onset, LSV can withstand the Lorentz force leading to stable dynamos. Furthermore, in the regime of $Po < 0.1$, stable dynamos without LSV have been seen only with $Pm \ge 1$. For the parameters investigated here, the presence of LSV allows low $Pm$ dynamos, while small-scale turbulent fluctuations are detrimental to dynamo action. Our results suggest that LSV play a key role in magnetic field generation in (spherical) planetary cores.

Despite the large number of simulations, predictive scaling laws remain elusive for the kinetic energy stored in the instabilities, for the onset of dynamo action, and for the magnetic field strength. These last two quantities have a non-monotonic behaviour when varying the key control parameters $E$ and $Po$, and an asymptotic regime has yet to be reached. Nore et al. (2014) and Goepfert & Tilgner (2018) also draw similar conclusions for dynamos in a precessing cylinder and cube, respectively. This suggests that the non-monotonic behaviour is not linked to the spherical shape, but rather to the precession itself. For $Pm \le 2$, the magnetic energy seems capped by the turbulent kinetic energy, in contrast to convective dynamos where magnetic energy overcomes kinetic energy as viscosity is lowered (e.g. Schaeffer et al., 2017). Predicting the turbulent fluctuation level is an outstanding issue, as only a few of our simulations reach the regime where departures from laminar dissipation are observed. Exploring this challenging regime relevant for planetary cores is however a necessary first step towards extracting useful scaling laws – if such laws exist for this system. We release our simulation database to enable further investigations and contributions.

None of our dynamos produce a predominantly dipolar field. Furthermore, at small Ekman number, while the magnetic field is generated in the bulk at the large scale of the LSV, the surface magnetic field is more and more dominated by small scales. This shift towards small-scale surface fields as the viscosity is lowered was observed for other spherical dynamos driven by the boundaries (Monteux et al., 2012). Indeed, intense small-scale turbulence develops in the boundary layers, shredding the magnetic field to small scales. This contrasts with convective dynamos for which a wide range of parameters lead to dipolar fields (e.g. Kutzner & Christensen, 2002).

For a precessing spheroid, a topography driven instability may occur (Vidal & Cébron, 2017) while the boundary layers remain stable. In this case, no small-scale turbulence would shred the magnetic field in the boundary layer, permitting a surface field on the same scale as in the bulk.
Acknowledgments

We thank Yufeng Lin for sharing his data with us. We are also grateful to Mathieu Dumberry and two other anonymous reviewers. Numerical simulations were performed using the XSHELLS code which is freely available at https://bitbucket.org/nschaeff/xshells. Our dataset built from our simulations is available at https://doi.org/10.6084/m9 figshare.7017137.

This work was performed using HPC resource Occigen from CINES under allocations x2016047382, A0020407382 and A0040407382 made by GENCI. Part of the computations were also performed on the Froggy platform of CIMENT (https://ciment.ujf-grenoble.fr), supported by the Rhône-Alpes region (CPER07_13 CIRA), OSUG@2020 LabEx (ANR10 LABX56) and Equip@Meso (ANR10 EQPX-29-01). NS and DC were supported by the French *Agence Nationale de la Recherche* under grant ANR-14-CE33-0012 (MagLune). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). RL would like to acknowledge support from the European Research Council (ERC Advanced grant 670874 ROTANUT) and the HPC resources of the Royal Observatory of Belgium.

References

Boisson, J., Cébron, D., Moisy, F., & Cortet, P.-P., 2012. Earth rotation prevents exact solid-body rotation of fluids in the laboratory, *EPL (Europhysics Letters)*, **98**(5), 59002.

Bowden, K., 1953. Note on wind drift in a channel in the presence of tidal currents, *Proc. R. Soc. Lond. A*, **219**(1139), 426-446.

Busse, F. H., 1968. Steady fluid flow in a precessing spheroidal shell, *Journal of Fluid Mechanics*, **33**(04), 739-751.

Cappanera, L., Guermond, J.-L., Léorat, J., & Nore, C., 2016. Two spinning ways for precession dynamo, *Physical Review E*, **93**(4), 043113.

Cébron, D., 2015. Bistable flows in precessing spheroids, *Fluid Dynamics Research*, **47**(2), 025504.

Cébron, D. & Hollerbach, R., 2014. Tidally driven dynamos in a rotating sphere, *The Astrophysical Journal Letters*, **789**(1), L25.

Cébron, D., Le Bars, M., & Meunier, P., 2010. Tilt-over mode in a precessing triaxial ellipsoid, *Physics of Fluids*, **22**(11), 116601.

Christensen, U. R., Holzwarth, V., & Reiners, A., 2009. Energy flux determines magnetic field strength of planets and stars, *Nature*, **457**(7226), 167.

Dwyer, C. A., Stevenson, D. J., & Nimmo, F., 2011. A long-lived lunar dynamo driven by continuous mechanical stirring, *Nature*, **479**(7372), 212-214.

Faller, A. J., 1991. Instability and transition of disturbed flow over a rotating disk, *Journal of Fluid Mechanics*, **230**, 245-269.

Fuchs, H., Rädler, K.-H., & Rheinhard, M., 1999. On self-killing and self-creating dynamos, *Astronomische Nachrichten: News in Astronomy and Astrophysics*, **320**(3), 129-133.

Giesecke, A., Vogt, T., Gundrum, T., & Stefani, F., 2018. Nonlinear large scale flow in a precessing cylinder and its ability to drive dynamo action, *Physical review letters*, **120**(2), 024502.

Goepfert, O. & Tilgner, A., 2016. Dynamos in precessing cubes, *New Journal of Physics*, **18**, 103019.

Goepfert, O. & Tilgner, A., 2018. Mechanisms for magnetic field generation in precessing cubes, *Geo-physical & Astrophysical Fluid Dynamics*, pp. 1-13.
Goto, S., Matsunaga, A., Fujiwara, M., Nishioka, M., Kida, S., Yamato, M., & Tsuda, S., 2014. Turbulence driven by precession in spherical and slightly elongated spheroidal cavities, *Physics of Fluids*, **26**(5), 055107.

Greenspan, H. P. G. et al., 1968. *The theory of rotating fluids*, CUP Archive.

Guervilly, C., Hughes, D. W., & Jones, C. A., 2015. Generation of magnetic fields by large-scale vortices in rotating convection, *Physical Review E*, **91**(4), 041001.

Hollerbach, R. & Kerswell, R., 1995. Oscillatory internal shear layers in rotating and precessing flows, *Journal of Fluid Mechanics*, **298**, 327–339.

Hollerbach, R., Nore, C., Marti, P., Vantieghem, S., Luddens, F., & Léorat, J., 2013. Parity-breaking flows in precessing spherical containers, *Physical Review E*, **87**(5), 053020.

Hough, S. S., 1895. The oscillations of a rotating ellipsoidal shell containing fluid, *Philosophical Transactions of the Royal Society of London A*, pp. 469–506.

Kerswell, R. R., 1995. On the internal shear layers spawned by the critical regions in oscillatory Ekman boundary layers, *Journal of Fluid Mechanics*, **298**, 311–325.

Kerswell, R. R., 1996. Upper bounds on the energy dissipation in turbulent precession, *Journal of Fluid Mechanics*, **321**, 335–370.

Kerswell, R. R., 2002. Elliptical instability, *Annual Review of Fluid Mechanics*, **34**(1), 83–113.

Kida, S., 2011. Steady flow in a rapidly rotating sphere with weak precession, *Journal of Fluid Mechanics*, **680**, 150–193.

Kida, S., 2018. Steady flow in a rotating sphere with strong precession, *Fluid Dynamics Research*, **50**(2), 021401.

Kutzner, C. & Christensen, U., 2002. From stable dipolar towards reversing numerical dynamos, *Physics of the Earth and Planetary Interiors*, **131**(1), 29–45.

Le Bars, M., Cébron, D., & Le Gal, P., 2015. Flows driven by libration, precession, and tides, *Annual Review of Fluid Mechanics*, **47**, 163–193.

Lemasquerier, D., Grannan, A. M., Vidal, J., Cébron, D., Favier, B., Le Bars, M., & Aurnou, J. M., 2017. Libration-driven flows in ellipsoidal shells, *Journal of Geophysical Research: Planets*, **122**(9), 1926–1950.

Lin, Y., Noir, J., & Jackson, A., 2014. Experimental study of fluid flows in a precessing cylindrical annulus, *Physics of Fluids*, **26**(4), 046604.

Lin, Y., Marti, P., & Noir, J., 2015. Shear-driven parametric instability in a precessing sphere, *Physics of Fluids*, **27**(4), 046601.

Lin, Y., Marti, P., Noir, J., & Jackson, A., 2016. Precession-driven dynamos in a full sphere and the role of large scale cyclonic vortices, *Physics of Fluids*, **28**(6), 066601.

Lingwood, R. J., 1997. Absolute instability of the ekman layer and related rotating flows, *Journal of Fluid Mechanics*, **331**, 405–428.

Lorenzani, S., 2001. *Fluid instabilities in precessing ellipsoidal shells*, Ph.D. thesis, Niedersächsische Staats-und Universitätsbibliothek Göttingen.

Malkus, W. V. R., 1968. Precession of the Earth as the cause of geomagnetism, *Science*, **160**(3825), 259–264.
Marti, P., Schaeffer, N., Hollerbach, R., Cébron, D., Nore, C., Luddens, F., Guermond, J.-L., Aubert, J., Takehiro, S., Sasaki, Y., Hayashi, Y.-Y., Simitev, R., Busse, F., Vantieghem, S., & Jackson, A., 2014. Full sphere hydrodynamic and dynamo benchmarks, *Geophysical Journal International*, **197**(1), 119–134.

Matsui, H., Heien, E., Aubert, J., Aurnou, J. M., Avery, M., Brown, B., Buffett, B. A., Busse, F., Christensen, U. R., Davies, C. J., et al., 2016. Performance benchmarks for a next generation numerical dynamo model, *Geochemistry, Geophysics, Geosystems*, **17**(5), 1586–1607.

Miralles, S., Plihon, N., & Pinton, J.-F., 2015. Lorentz force effects in the bullard–von kármán dynamo: saturation, energy balance and subcriticality, *Journal of Fluid Mechanics*, **775**, 501–523.

Moffatt, H., 1970. Turbulent dynamo action at low magnetic reynolds number, *Journal of Fluid Mechanics*, **41**(2), 435–452.

Monteux, J., Schaeffer, N., Amit, H., & Cardin, P., 2012. Can a sinking metallic diapir generate a dynamo?, *Journal of Geophysical Research: Planets*, **117**(E10).

Noir, J. & Cébron, D., 2013. Precession-driven flows in non-axisymmetric ellipsoids, *Journal of Fluid Mechanics*, **737**, 412–439.

Noir, J., Jault, D., & Cardin, P., 2001. Numerical study of the motions within a slowly precessing sphere at low Ekman number, *Journal of Fluid Mechanics*, **437**, 283–299.

Noir, J., Cardin, P., Jault, D., & Masson, J.-P., 2003. Experimental evidence of non-linear resonance effects between retrograde precession and the tilt-over mode within a spheroid, *Geophysical Journal International*, **154**(2), 407–416.

Nore, C., Léorat, J., Guermond, J., & Luddens, F., 2011. Nonlinear dynamo action in a precessing cylindrical container, *Physical Review E*, **84**(1), 016317.

Nore, C., Lorat, J., Guermond, J. L., Cappanera, L., & Luddens, F., 2014. Dynamo action in precessing cylinders, in *9th PAMIR International Conference Fundamental and Applied MHD*.

Oruba, L. & Dormy, E., 2014. Predictive scaling laws for spherical rotating dynamos, *Geophysical Journal International*, **198**(2), 828–847.

Poincaré, H., 1910. Sur la précession des corps déformables, *Bulletin Astronomique, Serie I*, **27**, 321–356.

Pozzo, M., Davies, C., Gubbins, D., & Alfe, D., 2012. Thermal and electrical conductivity of iron at Earth’s core conditions, *Nature*, **485**(7398), 355–361.

Rieutord, M., 2001. Ekman layers and the damping of inertial r-modes in a spherical shell: application to neutron stars, *The Astrophysical Journal*, **550**(1), 443.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations, *Geochemistry, Geophysics, Geosystems*, **14**(3), 751–758.

Schaeffer, N., Jault, D., Nataf, H.-C., & Fournier, A., 2017. Turbulent geodynamo simulations: a leap towards Earths core, *Geophysical Journal International*, **211**(1), 1–29.

Sloudsky, T., 1895. De la rotation de la Terre supposée fluide à son intérieur, *Bull. Soc. Imp. Natur. Mosc.*, **IX**, 285–318.

Sous, D., Sommeria, J., & Boyer, D., 2013. Friction law and turbulent properties in a laboratory ekman boundary layer, *Physics of Fluids*, **25**(4), 046602.

Springel, V., 2010. E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh, *Monthly Notices of the Royal Astronomical Society*, **401**(2), 791–851.

Stegman, D., Jellinek, A., Zatman, S., Baumgardner, J., & Richards, M., 2003. An early lunar core dynamo driven by thermochemical mantle convection, *Nature*, **421**(6919), 143–146.
Stewartson, K. & Roberts, P. H., 1963. On the motion of liquid in a spheroidal cavity of a precessing rigid body, *Journal of Fluid Mechanics*, **17**(01), 1–20.

Tilgner, A., 2005. Precession driven dynamos, *Physics of Fluids*, **17**(3), 034104.

Tilgner, A., 2007. Kinematic dynamos with precession driven flow in a sphere, *Geophysical & Astrophysical Fluid Dynamics*, **101**(1), 1–9.

Tilgner, A. & Busse, F. H., 2001. Fluid flows in precessing spherical shells, *Journal of Fluid Mechanics*, **426**, 387–396.

Touma, J. & Wisdom, J., 1994. Evolution of the earth-moon system, *The Astronomical Journal*, **108**, 1943–1961.

Vidal, J. & Cébron, D., 2017. Inviscid instabilities in rotating ellipsoids on eccentric Kepler orbits, *Journal of Fluid Mechanics*, **833**, 469–511.

Vidal, J., Cébron, D., Schaeffer, N., & Hollerbach, R., 2018. Magnetic fields driven by tidal mixing in radiative stars, *Monthly Notices of the Royal Astronomical Society*, **475**(4), 4579–4594.

Vormann, J. & Hansen, U., 2018. Numerical simulations of bistable flows in precessing spheroidal shells, *Geophysical Journal International*, **213**(2), 786–797.

Weber, R. C., Lin, P.-Y., Garnero, E. J., Williams, Q., & Lognonne, P., 2011. Seismic detection of the lunar core, *science*, **331**(6015), 309–312.

Williams, J. G. & Boggs, D. H., 2015. Tides on the moon: Theory and determination of dissipation, *Journal of Geophysical Research: Planets*, **120**(4), 689–724.

Williams, J. G., Boggs, D. H., Yoder, C. F., Ratcliff, J. T., & D., J. O., 2001. Lunar rotational dissipation in solid body and molten core, *Journal of Geophysical Research: Planets*, **106**(E11), 27933–27968.

Williams, J. G., Konopliv, A. S., Boggs, D. H., Park, R. S., Yuan, D.-N., Lemoine, F. G., Goossens, S., Mazarico, E., Nimmo, F., Weber, R. C., et al., 2014. Lunar interior properties from the grail mission, *Journal of Geophysical Research: Planets*, **119**(7), 1546–1578.

Wu, C.-C. & Roberts, P. H., 2013. On a dynamo driven topographically by longitudinal libration, *Geophysical & Astrophysical Fluid Dynamics*, **107**(1-2), 20–44.

Yoder, C. F., 1981. The free librations of a dissipative moon, *Phil. Trans. R. Soc. Lond. A*, **303**(1477), 327–338.

Yoder, C. F., 1995. Venus' free obliquity, *Icarus*, **117**(2), 250–286.

Zhang, K., Chan, K. H., & Liao, X., 2010. On fluid flows in precessing spheres in the mantle frame of reference, *Physics of Fluids*, **22**(11), 116604.

## A Theoretical solutions for the forced base flow

In this appendix, we consider the usual planetary relevant limit $Po \ll 1$ considered in the literature. In this limit, the unit of time $\Omega_o^{-1}$ corresponds to the usual unit of time $\Omega_s^{-1}$.
## A.1 The model of Busse (1968)

In the frame of precession, the three Cartesian components ($\Omega_x, \Omega_y, \Omega_z$) of the dimensionless fluid rotation vector $\Omega$ are governed by the three following theoretical equations (e.g. Noir et al. (2003); Cébron et al. (2010))

$$ \Omega_z = \Omega_x^2 + \Omega_y^2 + \Omega_z^2, \qquad (34) $$

$$ -P_z \Omega_y = (\lambda_r \Omega_x \Omega_z^{1/4} + \lambda_i \Omega_y \Omega_z^{-1/4}) \sqrt{E}, \qquad (35) $$

$$ P_x \Omega_y = -\lambda_r \Omega_z^{1/4} (1 - \Omega_z) \sqrt{E}, \qquad (36) $$

with $P_x = P_0 \sin \alpha$ and $P_z = P_0 \cos \alpha$ the two dimensionless components of the precession vector $\Omega_p$ along x and z (the y-component is zero). Equations (34)-(36) are exactly the equations (20)-(22) of Noir et al. (2003), or equations (21)-(23) of Cébron et al. (2010) in the particular case of a sphere (no deformation, and $P_y = 0$ in their equations). As shown by Noir et al. (2003), this system of equations is equivalent to the well-known implicit expression (3.19) of Busse (1968). Equation (34) is the so-called no spin-up condition (solvability condition 3.14 of Busse (1968), or equation (12) of Noir et al. (2003)) given that it forbids any differential rotation along $\Omega$. Equations (35)-(36) are simply obtained from a torque balance (see Noir et al. (2003); Cébron et al. (2010) for details).

In equations (34)-(36), we have noted the spin-over damping factor $\lambda = \lambda_r + i\lambda_i$, given by

$$ \lambda_{inv}^{sphere} = - \frac{3[19(1-i) + 9\sqrt{3}(1+i)]}{28\sqrt{2}} \approx -2.62 + 0.258i \qquad (37) $$

for a spherical container ($\eta = 0$), in the inviscid limit $E \ll 1$. For finite values of $E$, Noir et al. (2001) has obtained empirically $\lambda_r \approx -2.62 - 1.36E^{0.27}$, and a fit of the results of Hollerbach & Kerswell (1995) gives $\lambda_i \approx 0.258 + 1.25E^{0.21}$.

Even if equations (34)-(36) are obtained without any inner core, corrections have been proposed for the case of the sphere to take a spherical inner core into account. Using the dimensionless inner radius $r_i$, it has been proposed to simply modify $\lambda$ by the factor $(1 + \eta^4)/(1 - \eta^5)$ for a no-slip inner core Hollerbach & Kerswell (1995), and by $1/(1 - \eta^5)$ for a stress-free inner core Tilgner & Busse (2001). Using the results of Hollerbach & Kerswell (1995) for the spherical shell, the corrections for viscous and aspect ratio effects can be combined following

$$ \lambda \approx \left[ \lambda_{inv}^{sphere} - 1.36E^{0.27} + i1.25E^{0.21} \right] \frac{1 + \eta^4}{1 - \eta^5}, \qquad (38) $$

with $\lambda = \lambda_r + i\lambda_i$.

## A.2 Approximate explicit solution

As shown by Noir & Cébron (2013), (34)-(36) can be obtained as fixed points of a dynamical model for $\Omega$, given by (see equations A14-A 16 of Noir & Cébron (2013) for a spheroid)

$$ \frac{\partial \Omega_x}{\partial t} = P_z \Omega_y - (1 - \gamma) [P_z \Omega_y + \Omega_y \Omega_z] + \mathcal{L} \Gamma_{\nu} \cdot \hat{x}, \qquad (39) $$

$$ \frac{\partial \Omega_y}{\partial t} = P_x \Omega_z - P_z \Omega_x + (1 - \gamma) [P_z \Omega_x + \Omega_x \Omega_z] + \mathcal{L} \Gamma_{\nu} \cdot \hat{y}, \qquad (40) $$

$$ \frac{\partial \Omega_z}{\partial t} = -P_x \Omega_y - (1 - \gamma) P_x \Omega_y + \mathcal{L} \Gamma_{\nu} \cdot \hat{z}, \qquad (41) $$

where $\gamma = (2a^2)/(a^2 + c^2)$ represents the ratio of the polar to equatorial moment of inertia, where $\Gamma_\nu$ is the viscous torque, and where $\mathcal{L}$ is a 3x3 matrix given by equation (A3) of Noir & Cébron (2013). For the spherical shell, $\mathcal{L}$ reduces to $\mathcal{L} = 15\delta_{ij}/(8\pi(1 - \eta^5))$, using the Kronecker delta $\delta_{ij}$. As detailed by Cébron (2015), equations (34)-(36) are then recovered with

$$ \mathcal{L} \Gamma_{\nu} = \sqrt{\Omega E} \left[ \lambda_r \begin{pmatrix} \Omega_x \\ \Omega_y \\ \Omega_z - 1 \end{pmatrix} + \frac{\lambda_i}{\Omega} \begin{pmatrix} \Omega_y \\ -\Omega_x \\ 0 \end{pmatrix} \right]. \qquad (42) $$
To obtain tractable analytical solutions, we need a simpler set of equations. We thus follow Noir & Cébron (2013) who linearize $\mathcal{L}\Gamma_ν$ by assuming that the fluid rotates at the same rate than the boundaries, i.e. $\Omega = 1$, which gives

$$ \mathcal{L}\Gamma_ν = \sqrt{E} \left[ \lambda_r \begin{pmatrix} \Omega_x \\ \Omega_y \\ \Omega_z - 1 \end{pmatrix} + \lambda_i \begin{pmatrix} \Omega_y \\ -\Omega_x \\ 0 \end{pmatrix} \right]. \quad (43) $$

Focusing on stationary solutions of equations (39)-(41) in the sphere ($\gamma = 1$) with the viscous term (43), the explicit solution for $\Omega$ is then

$$ \Omega_x = \frac{[\lambda_i + \chi \cos \alpha] \chi \sin \alpha}{\chi(\chi + 2\lambda_i \cos \alpha) + |\Delta|^2}, \quad (44) $$

$$ \Omega_y = -\frac{\chi \lambda_r \sin(\alpha)}{\chi(\chi + 2\lambda_i \cos \alpha) + |\Delta|^2}, \quad (45) $$

$$ \Omega_z = \frac{\chi(\chi \cos^2 \alpha + 2\lambda_i \cos \alpha) + |\Delta|^2}{\chi(\chi + 2\lambda_i \cos \alpha) + |\Delta|^2}, \quad (46) $$

where $\chi = Po/\sqrt{E}$, and $|\Delta|^2 = \lambda_r^2 + \lambda_i^2$. The differential rotation $\epsilon$ of the fluid with the boundary is thus:

$$ \epsilon = \frac{|\chi \sin \alpha|}{\sqrt{\chi(\chi + 2\lambda_i \cos \alpha) + |\Delta|^2}}. \quad (47) $$

Note that these explicit expressions for $\Omega$ allows to clarify quantitatively various observations, previously noticed in the literature. For instance, Noir & Cébron (2013) have considered the resonance of $\epsilon$ for a fixed Rossby number $Ro = Po \sin \alpha$, i.e. the value of $Po$ where $\epsilon$ is maximum when $Ro$ is maintained constant. They have noticed that the role of $\lambda_i$ is mainly to shift the resonance peak, that the reduced model (43) always gives at $Po = 0$ for the sphere when $\lambda_i = 0$ is assumed. We can thus use expressions (44)-47 to clarify this observation. For a given $Ro$, we obtain that the resonance is indeed reached for $Po = 0$ when $\lambda_i = 0$, but the resonance is shifted to

$$ \chi = Po/\sqrt{E} = -\lambda_i \sqrt{1 + Ro^2/\lambda_i^2} \quad (48) $$

when $\lambda_i \neq 0$. It is straightforward to show that, at the resonance, $\Omega_x$ is zero whereas $\Omega_y$ and $\Omega_z$ are respectively maximum and minimum.

One can finally notice that the additional hypothesis $\lambda_i = 0$, considered e.g. by Noir & Cébron (2013) and Cébron (2015), allows to simplify equations (44)-(46) into

$$ \Omega = \left( \frac{1}{2} \frac{\sin 2\alpha}{1+x^2}, -\frac{x \sin \alpha}{1+x^2}, \frac{x^2 + \cos^2 \alpha}{1+x^2} \right), \quad (49) $$

which gives

$$ \epsilon = \frac{|\sin \alpha|}{\sqrt{1+x^2}}, \quad (50) $$

where $x = \lambda_r/\chi$. Note that $\Omega_x$, $\Omega_z$ and $\epsilon$ are maximum for $x = 0$, whereas $\Omega_y$ is maximum for $x = \pm 1$. It is interesting to note that equation (50) is exactly equation (53) of Williams et al. (2001), but we manage here to obtain an explicit expression for $x$.