2.1 Model equations, input parameters and methods

We consider a convecting, electrically conducting, incompressible fluid in a self-gravitating spherical shell between radii $r_{i}$ and $r_{o}$ , with $r_{i}/r_{o}=0.35$ as in the Earth’s core at present. The shell is rotating about an axis $\boldsymbol{e}_{z}$ with an angular velocity $\unicode[STIX]{x1D6FA}$ . We solve for Boussinesq convection for a velocity field $\boldsymbol{u}$ and a density anomaly field $C$ , magnetic induction in the magnetohydrodynamic approximation for a magnetic field $\boldsymbol{B}$ , with the following set of dimensionless equations:

Here $\boldsymbol{r}$ is the radius vector. Time is scaled with the inverse of the rotation rate $\unicode[STIX]{x1D6FA}^{-1}$ . Length is scaled with the shell gap $D=r_{o}-r_{i}$ . Velocity is scaled with $\unicode[STIX]{x1D6FA}D$ . Magnetic induction is scaled by $(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D707})^{1/2}\unicode[STIX]{x1D6FA}D$ , where $\unicode[STIX]{x1D70C}$ is the fluid density and $\unicode[STIX]{x1D707}$ the magnetic permeability of the fluid. Tables 1 and 2 summarise the values of the input parameters for our survey, together with their Earth estimates. We recall the definitions for the Ekman number $E$ , magnetic Prandtl number $Pm$ and introduce the Prandtl number $Pr$  as:

Here $\unicode[STIX]{x1D708},\unicode[STIX]{x1D702}$ have already been introduced as the fluid viscosity and magnetic diffusivity, and $\unicode[STIX]{x1D705}$ is the diffusivity of the density anomaly field $C$ . Two different sets of boundary conditions are investigated. The standard type (ST) refers to the standard boundary conditions often used in the numerical dynamo modelling community (e.g. Christensen & Aubert Reference Christensen and Aubert2006): rigid (no-slip) and electrically insulating boundary conditions, and fixed density anomaly at both boundaries, with an imposed difference $\unicode[STIX]{x0394}C$ between inner and outer boundary (which also serves as the scale for the dimensionless density anomaly). In this case the density anomaly source term $S$ vanishes in (2.3), and the modified Rayleigh number $Ra^{\ast }$ appearing in (2.1) reads:

Here $g_{o}$ is the gravity at the outer boundary. The modified Rayleigh number relates to the canonical Rayleigh number $Ra=g_{o}\unicode[STIX]{x0394}CD^{3}/\unicode[STIX]{x1D70C}\unicode[STIX]{x1D705}\unicode[STIX]{x1D708}$ through

Note that the square root of the modified Rayleigh number is classically referred to as the convective Rossby number (Gilman Reference Gilman1977). The second set of boundary conditions is derived from the coupled Earth (CE)-type dynamo set-up (Aubert et al. Reference Aubert, Finlay and Fournier2013). Here the idea is to minimise the residual influence of boundary layers, and provide a geophysically more relevant description of global coupling by using stress-free, electrically conducting, and fixed mass anomaly flux conditions at both boundaries. We will show in § 4 that results obtained within set-ups ST and CE are largely similar, meaning that the nature of the asymptotic state is not influenced by the choice of boundary conditions. In type CE, the imposed mass anomaly flux crossing the shell is defined as $F$ , and the density anomaly scale for non-dimensionalisation is defined as $F/4\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FA}D^{3}$ . The dimensionless, spherical parts of the mass anomaly flux at the inner and outer boundaries are fixed to $\unicode[STIX]{x2202}C/\unicode[STIX]{x2202}r(r_{i})=-Pr/\,Er_{i}^{2}$ and $\unicode[STIX]{x2202}C/\unicode[STIX]{x2202}r(r_{o})=0$ . The modified flux Rayleigh number appearing in (2.1) then becomes

and mass flux conservation implies that the source term in (2.3) for type CE is $S=-3/(r_{o}^{3}-r_{i}^{3})$ (Aubert, Labrosse & Poitou Reference Aubert, Labrosse and Poitou2009). Throughout this work, the notation $Ra_{(F)}^{\ast }$ introduced in (2.1) refers either to $Ra^{\ast }$ or $Ra_{F}^{\ast }$ depending on whether we consider the ST or the CE set-up, respectively.

Table 2.

Set of numerical models with input and output parameters (see § 2 for definitions), together with estimates for Earth’s core (see § 3). Type CE refers to numerical models performed in the coupled Earth set-up (fixed mass anomaly flux, stress-free, electrically conducting boundary conditions). Type ST refers to models performed in the standard set-up (fixed density anomaly, no-slip, electrically insulating boundary conditions). Types D and L respectively refer to direct and large-eddy simulations (see § 3.4). Model 0 is similar to the original CE dynamo (Aubert et al. Reference Aubert, Finlay and Fournier2013, parameters with subscript 0 in the text). Model 1 is an equivalent standard dynamo (parameters with subscript 1 in the text). All models have $r_{i}/r_{o}=0.35$ and $Pr=1$ . All large-eddy simulations have $l_{h}=30$ . Note that the magnetic Reynolds number $Rm$ is constant to within $\pm 10\,\%$ , and that the Elsasser number $\unicode[STIX]{x1D6EC}$ for a given boundary condition type (CE or ST) is also constant to within $\pm 10\,\%$ .

In addition to this basic set-up, type CE also features a number of refinements, listed below, which help the simulation output to match the statics and kinematics of the geomagnetic field, particularly regarding the geomagnetic westward drift and the geographical localisation of the temporal magnetic field variations. Although these are mostly irrelevant to the fundamental force balance and scaling results to be presented in the following sections, they have been kept here for future exploitation of the results in a geophysical context. The fluid shell is magnetically coupled to a solid inner sphere of radius $r_{i}$ (the inner core) and with the same electrical conductivity, axially rotating at a rotation rate $\unicode[STIX]{x1D6FA}_{ic}$ . The inner sphere couples to a solid spherical shell between radii $r_{o}$ and $1.83r_{o}$ (the mantle), also axially rotating at a rate $\unicode[STIX]{x1D6FA}_{m}$ , through a restoring remote torque $\unicode[STIX]{x1D6E4}=\unicode[STIX]{x1D709}(\unicode[STIX]{x1D6FA}_{ic}-\unicode[STIX]{x1D6FA}_{m})$ meant to arise from gravitational coupling, with $\unicode[STIX]{x1D709}=0.75$ . The outer solid shell has an dimensionless electrical conductance $\unicode[STIX]{x1D6F4}=10^{-4}$ also enabling magnetic coupling with the fluid shell (see Pichon, Aubert & Fournier Reference Pichon, Aubert and Fournier2016, equations (28) and (29), for complete definitions of $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D6F4}$ ). Angular momentum of the coupled system is preserved, with moments of inertia for the inner sphere, fluid shell, and outer solid shell in Earth-like ratios (see also Pichon et al. Reference Pichon, Aubert and Fournier2016). Finally, a longitudinal hemispheric modulation of the mass anomaly flux is imposed at the inner core boundary on top of the homogeneous part, with a peak-to-peak amplitude $0.8Pr/\,Er_{i}^{2}$ . A spatial modulation of the mass anomaly flux at the outer boundary is also imposed, with a peak-to-peak amplitude $0.1Pr/\,E$ . The patterns and relative orientations of these mass anomaly flux heterogeneities are those prescribed in Aubert et al. (Reference Aubert, Finlay and Fournier2013).

We have carried out 10 numerical cases of type CE and 7 cases of type ST (table 2). Type ST cases were integrated using the MagIC numerical implementation (Wicht Reference Wicht2002; Gastine et al. Reference Gastine, Wicht, Barik, Putigny and Duarte2016, freely available at www.github.com/magic-sph/magic), which uses a Chebyshev decomposition in the radial direction with $NR$ nodal points, and a lateral spherical harmonic decomposition up to degree and order $l_{max}$ . Type CE cases were integrated using the PARODY-JA numerical implementation (Dormy, Cardin & Jault Reference Dormy, Cardin and Jault1998; Aubert, Aurnou & Wicht Reference Aubert, Aurnou and Wicht2008) using a second-order finite differencing scheme in the radial direction with $NR$ grid points in the fluid, 50 grid points in the inner sphere and the same spherical harmonic decomposition as MagIC. Both implementations use a toroidal–poloidal decomposition of the solenoidal vector fields $\boldsymbol{u}$ and $\boldsymbol{B}$ (2.4), (2.5), the same Crank–Nicolson/Adams–Bashforth semi-implicit adaptive time stepping scheme of second order in time and are benchmarked against each other (Christensen et al. Reference Christensen, Aubert, Busse, Cardin, Dormy, Gibbons, Glatzmaier, Honkura, Jones and Kono2001). Both implementations resort to the library SHTns (Schaeffer Reference Schaeffer2013, freely available at https://bitbucket.org/nschaeff/shtns), for efficient handling of the spherical harmonic transforms. All calculations have been performed in the full domain and no azimuthal symmetry was assumed. More extreme parameter regimes could be reached in the CE set-up, owing to the absence of viscous and density anomaly boundary layers. The least forced reference models 0 and 1 in the ST and CE set-ups were already largely supercritical, with Rayleigh numbers respectively 40 and 400 times their critical values. Note however that these figures are different because supercriticality is defined relative to the flux in CE case, and the dynamical regimes of models 0 and 1 are otherwise equivalent. All models produced a self-sustained magnetic field with a dominant axial dipole that did not reverse polarity, and with a morphology exhibiting detailed resemblance to that of the geomagnetic field (see § 4.2). No bistable states were found, though such states were not specifically searched for. Our initial conditions for models 0 and 1 were indeed an axially dipolar magnetic field, zero velocity field and a random density anomaly perturbation. We then progressed by half-decades towards lower Ekman numbers by iteratively initialising each new model with the final state of the previous calculation. Large-eddy simulations (§ 3.4) were computed first, and then expanded into fully resolved numerical simulations by removing the hyperdiffusive treatment, a strategy causing only weak transients and ensuring significant savings in numerical cost.

2.2 Output parameters

In the following, we will analyse the following time-averaged, integral outputs, which are also summarised in table 1 and detailed in table 2. The magnetic field amplitude in the shell is characterised either by the Lehnert number

or the Elsasser number

where $B$ is the dimensional, root-mean-squared magnetic field amplitude in the fluid shell. Table 2 lists the Elsasser numbers in our numerical cases, and we recall that $\unicode[STIX]{x1D706}=\sqrt{\unicode[STIX]{x1D6EC}E/Pm}$ . The velocity field amplitude is characterised either by the Rossby number

the hydrodynamic Reynolds number

or the magnetic Reynolds number

where $U$ is the dimensional, root-mean-squared velocity field amplitude in the fluid shell. Table 2 lists the magnetic Reynolds numbers, and we recall that $Re=Rm/Pm$ and $Ro=ERm/Pm$ . The squared Alfvén number

measuring the ratio of kinetic to magnetic energy can also be derived from table 2 as $A^{2}=ERm^{2}/Pm\unicode[STIX]{x1D6EC}$ . The efficiency of mass anomaly transfer in the ST cases is measured by the Nusselt number $Nu=FD/4\unicode[STIX]{x03C0}r_{i}r_{o}\unicode[STIX]{x1D705}\unicode[STIX]{x0394}C$ .

The volumetric convective power is defined as

where $V=4\unicode[STIX]{x03C0}(r_{o}^{3}-r_{i}^{3})/3$ is the volume of the fluid shell. On time average, the convective power equates to the rate of gravitational potential energy release, which itself relates to the mass anomaly flux. In the CE set-up this leads to a relationship between $p$ and $Ra_{F}^{\ast }$ (Aubert et al. Reference Aubert, Labrosse and Poitou2009) while in the ST set-up this may be more precisely assessed by relating $p$ to an advected flux $Ra^{\ast }(Nu-1)E/Pr$ (Christensen & Aubert Reference Christensen and Aubert2006). Determining the proportionality factor $\unicode[STIX]{x1D6FE}=p/Ra_{F}^{\ast }$ (CE set-up) or $\unicode[STIX]{x1D6FE}=p/(Ra^{\ast }(Nu-1)E/Pr)$ (ST set-up) requires knowledge of the gravitational potential difference between the radii for density anomaly injection and mixing. The constant $\unicode[STIX]{x1D6FE}$ can be exactly determined in a configuration of condensed mass central gravity $g\sim 1/r^{2}$ (Gastine, Wicht & Aurnou Reference Gastine, Wicht and Aurnou2015). Approximate values of $\unicode[STIX]{x1D6FE}$ depending only on geometry can also be analytically determined in the present situation of sufficiently supercritical convection and linear radial gravity, through an assumption on the density anomaly profile (ST set-up, Christensen & Aubert Reference Christensen and Aubert2006) or on the average gravitational potential relevant to density anomaly mixing (CE set-up, Aubert et al. Reference Aubert, Labrosse and Poitou2009). Figure 1 shows that the deviations of our simulations from these approximate theories remain weak throughout the investigated range and should also remain negligible at Earth’s core conditions, as expected for strongly supercritical convection in the rotationally dominated regime (Oruba Reference Oruba2016). In the following, we will thus consider that there is proportionality between $p$ and $Ra_{F}^{\ast }$ in the CE set-up, such that power is an input parameter. In the ST set-up we will also consider that $p$ is proportional to $Ra^{\ast }(Nu-1)E/Pr$ , meaning that power is an output parameter.

Figure 1.

Proportionality (a) between the convective power $p$ and the dimensionless mass anomaly flux $Ra_{F}^{\ast }$ (CE cases) or $Ra^{\ast }(Nu-1)E/Pr$ (ST cases), presented together with its compensated scaling (b). The green line is the analytic prediction relevant to the CE set-up from the perfect-mixing theory (Aubert et al. Reference Aubert, Labrosse and Poitou2009), with proportionality constant $\unicode[STIX]{x1D6FE}=0.332$ . The orange line is the analytic prediction $\unicode[STIX]{x1D6FE}=0.480$ obtained in the ST set-up by Christensen & Aubert (Reference Christensen and Aubert2006). In (b), shaded regions represent the $\pm 1$ standard deviation of temporal fluctuations relative to the time average. In the CE set-up, the large-eddy simulations (see § 3.4) exhibit a weak residual trend for $\unicode[STIX]{x1D6FE}$ which should amount to approximately $20\,\%$ at Earth’s core value $Ra_{F}^{\ast }=O(10^{-12})$ , a negligible error compared to the geophysical uncertainty in the determination of $Ra_{F}^{\ast }$  (3.14).

The volumetric rate of ohmic dissipation is defined as

and is expressed in table 2 through the ohmic dissipation fraction $f_{ohm}=D_{\unicode[STIX]{x1D702}}/p$ . The equivalent ohmic dissipation length as defined in Christensen & Tilgner (Reference Christensen and Tilgner2004) is also reported in table 2 and reads

We finally use a standard measure (Wicht & Christensen Reference Wicht and Christensen2010; Teed et al. Reference Teed, Jones and Tobias2015) for the enforcement of the Taylor constraint (Taylor Reference Taylor1963), by assessing the cancellation level of the Lorentz force acting on axial cylinders:

Here $s,\unicode[STIX]{x1D711},z$ are cylindrical coordinates, $\boldsymbol{e}_{\unicode[STIX]{x1D711}}$ is the unit vector in the azimuthal direction, and the overbar represents an average taken in the azimuthal direction. The vertical integrals are evaluated between the lower and upper heights $z_{-,+}$ at which the axial cylinder of radius $s$ intersects the spherical shell boundaries. In table 2 we report ${\mathcal{T}}$ as an average of $|{\mathcal{T}}(s)|$ over axial cylinders between $s=0$ and $s=r_{o}$ , hence also encompassing fluid regions above and below the inner sphere inside the axial cylinder tangent to the inner core (the tangent cylinder).

3.1 Outline of the MAC theory

The D13 theory introduces three length scales to describe the asymptotic solutions. First, $d_{/\!/}$ describes the typical size of convection structure along the rotation axis. The second scale $d_{\bot }$ is the typical size of convection structures in a plane perpendicular to the rotation axis. The third scale is finally the magnetic dissipation length $d_{min}$ which we introduced above in its dimensionless form. The columnar structure of convection rolls is a robust feature of spherical convective dynamos (Christensen & Aubert Reference Christensen and Aubert2006; Yadav et al. Reference Yadav, Gastine, Christensen and Duarte2013b ), prompting

The MAC balance evaluated from the curled, dimensional version of (2.1) leads to

where $C$ denotes a typical density anomaly, complementing the already introduced typical velocity and magnetic fields $U$ and $B$ . Likewise, the dimensional balance between the rate of convective energy input and magnetic dissipation (Christensen & Tilgner Reference Christensen and Tilgner2004; Christensen & Aubert Reference Christensen and Aubert2006) may now be written

where we have used the equivalence between convective power and mass anomaly flux (figure 1). Using $F\sim UCD^{2}$ , and assuming a context of dominant magnetic dissipation i.e. $f_{ohm}\approx 1$ , one sees that satisfying both the balance between buoyancy and magnetic forces in (3.2) and the energy budget (3.3) implies that

This implies that the large-scale vorticity and the vorticity at the scale of magnetic dissipation are equivalent. This is an important result as it indicates that the turbulent energy transfers which are essential to the dynamics occur between the large scale $d_{\bot }$ at which the magnetic field sustains itself by drawing convective power (balance between buoyancy and Lorentz forces), and the small scale $d_{\min }$ at which this power is ohmically dissipated. If one further requires the large- and small-scale vorticities to be independent of the system rotation rate and diffusivities, then from (3.3) it follows that the magnetic field itself is independent of the rotation rate and diffusivities. Dimensional analysis finally yields the following scaling, corresponding to the initial proposal of Christensen & Aubert (Reference Christensen and Aubert2006):

One finally obtains the dimensionless power-driven, diffusivity-free scalings

together with the following relationships for the length scales

In (3.8) $d_{min}$ thus primarily scales with the magnetic Reynolds number, in accordance with the initial result of Christensen & Tilgner (Reference Christensen and Tilgner2004). In the range between model 0 (table 2) where $Ro=10^{-2}$ , and the conditions of Earth’s core where $Ro\approx 3\times 10^{-6}$ (see § 3.2), the secondary dependence with the Rossby number $Ro^{1/8}$ should remain small, within a factor 3. Numerical dynamo surveys (Christensen Reference Christensen2010; Stelzer & Jackson Reference Stelzer and Jackson2013; Yadav et al. Reference Yadav, Gastine and Christensen2013a ) generally suggest an even weaker dependence of $d_{min}$ in $Ro$ with power-law exponents between 1/20 and 1/10. The discrepancy with D13 may be rooted in the variations of $f_{ohm}$ that are still significant in the numerical data set but not considered in the asymptotic MAC scaling theory. A reasonable approximation is thus to neglect the Rossby number dependence in (3.8). Consistency with the vorticity equivalence (3.4) then also demands the stronger assumption to neglect the Rossby number dependence in (3.9).

3.2 Path theory in the CE set-up

The calculations listed in table 2 do not represent a systematic sampling of the parameter space, but are rather chosen to follow a path connecting the classical numerical models such as the original CE dynamo (Aubert et al. Reference Aubert, Finlay and Fournier2013) to Earth’s core conditions. To introduce this concept it is best to use a reasoning based on time scale ratios, with our ultimate goal being to bring five relevant time scales in Earth-like proportions. These are the inverse rotation rate $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FA}}=1/\unicode[STIX]{x1D6FA}$ , the convective overturn time $\unicode[STIX]{x1D70F}_{U}=D/U$ , the Alfvén time $\unicode[STIX]{x1D70F}_{A}=\sqrt{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D707}}D/B$ , the time scale for convective power input $\unicode[STIX]{x1D70F}_{F}=(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70C}D^{4}/g_{o}F)^{1/3}$ (analogous to the Kelvin–Helmholtz time scale), and the magnetic diffusion time $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}=D^{2}/\unicode[STIX]{x1D702}$ . Less relevant to us are the viscous diffusion time $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}=D^{2}/\unicode[STIX]{x1D708}$ and the density anomaly diffusion time $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}}=D^{2}/\unicode[STIX]{x1D705}$ . For these latter times we do not target Earth-like ratios relative to the other times, but simply conditions corresponding to the magnetic diffusion time $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ being much smaller than $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}$ and $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}}$ , so that ohmic losses are the dominant source of dissipation.

We first illustrate the strategy with CE cases. Our starting point is a model similar to the original coupled Earth dynamo (Aubert et al. Reference Aubert, Finlay and Fournier2013), which we label model 0 (see table 2). In the following, subscript 0 denotes the parameters relative to this model. Four relevant time scale ratios for this model are

Our goal corresponds to conditions in the Earth’s core. To estimate these we use the values $B=4~\text{mT}$ (Gillet et al. Reference Gillet, Jault, Canet and Fournier2010), $U=5\times 10^{-4}~\text{m}~\text{s}^{-1}$ (e.g. Aubert Reference Aubert2014), $\unicode[STIX]{x1D70C}=11000~\text{kg}~\text{m}^{3}$ , $D=2260~\text{km}$ , $\unicode[STIX]{x1D6FA}=7.29\times 10^{-5}~\text{s}^{-1}$ , which are known with some certainty. We further use the less certain values $\unicode[STIX]{x1D702}=0.3{-}3~\text{m}^{2}~\text{s}^{-1}$ for the magnetic diffusivity, and $Q_{ad}=4{-}15~\text{TW}$ for the core–mantle boundary adiabatic heat flow, the uncertainty being due to an ongoing debate on core thermal and electrical conductivities (e.g. Pozzo et al. Reference Pozzo, Davies, Gubbins and Alfè2012; Konôpková et al. Reference Konôpková, McWilliams, Gómez–Pérez and Goncharov2016; Ohta et al. Reference Ohta, Kuwayama, Hirose, Shimizu and Ohishi2016). We assume that heat flow at the core–mantle boundary is exactly adiabatic, such that the dynamo is entirely bottom driven, as is the case in our CE models. Taking into account the thermodynamic efficiencies of inner core crystallisation and latent heat release (as computed in Lister Reference Lister2003), the dimensional, volumetric power is then $p\approx 0.2Q_{ad}/V$ . The modified, flux-based Rayleigh number can then be deduced using the proportionality exhibited in figure 1 as $Ra_{F}^{\ast }=p/\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FA}^{3}D^{2}$ , with $\unicode[STIX]{x1D6FE}=0.33$ . The resulting time scale ratios for Earth’s core are then

We see that along a sensible path connecting models and Earth, we would need to preserve at least the magnetic Reynolds number and the Elsasser number, our starting values representing reasonable midpoints in the range inferred for Earth’s core. We will see below that these two conditions in fact reduce to only one once the MAC balance is considered. We thus first decide to keep $Rm$ constant. From the discussion in § 3.1 this implies that $d_{min}$ should be close to invariant along the path. If we assume invariance for $d_{min}$ , then the vorticity equivalence (3.4) also implies that $d_{\bot }$ is invariant. The first part of the balance (3.2) then reduces to $\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FA}U\sim g_{o}C$ , and we also have $F\sim UCD^{2}$ . This yields the classical thermal-wind balance (e.g. Starchenko & Jones Reference Starchenko and Jones2002; Aurnou et al. Reference Aurnou, Andreadis, Zhu and Olson2003; Aubert Reference Aubert2005; Pichon et al. Reference Pichon, Aubert and Fournier2016), which is written here using the time scales introduced above:

With a constant $d_{min}$ and $f_{ohm}\approx 1$ , the energy budget (3.3) now becomes

Note that (3.18) and (3.19), together with a constant $Rm=\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D70F}_{U}$ , indeed yield a constant Elsasser number $\unicode[STIX]{x1D6EC}=\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FA}}/\unicode[STIX]{x1D70F}_{A}^{2}$ along the path. With these relationships at hand, we may now mathematically define the rules for input parameters along the path. Any model may be characterised by a path parameter $\unicode[STIX]{x1D716}$ such that

In the CE set-up, $\unicode[STIX]{x1D716}$ is then an input parameter controlling the mass anomaly flux and convective power, or alternatively the rotation rate since $Ra_{F}^{\ast }=(\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FA}}/\unicode[STIX]{x1D70F}_{F})^{3}$ . The rapidly rotating asymptotic regime is obtained for $\unicode[STIX]{x1D716}\rightarrow 0$ . With this definition, $\unicode[STIX]{x1D716}=1$ corresponds to the conditions of model 0, and $\unicode[STIX]{x1D716}=10^{-7}$ is appropriate to describe Earth’s core conditions given the Earth values reported above. Using the definition of $Ra_{F}^{\ast }$ in terms of time scales, together with (3.18) and a constant $Rm=\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}/\unicode[STIX]{x1D70F}_{U}$ we then obtain

We see here again that $\unicode[STIX]{x1D716}=1$ still describes model 0, while $\unicode[STIX]{x1D716}=10^{-7}$ indeed yields a correct Earth value $E/Pm=3.8\times 10^{-9}$ for the magnetic Ekman number. As mentioned above, we also wish to sufficiently increase $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}$ and $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}}$ relative to $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ in order to ensure a dominantly ohmic dissipation. We adopt

such that

We see that setting $\unicode[STIX]{x1D716}=10^{-7}$ implies $Pm=7.9\times 10^{-4}$ and $E=3\times 10^{-12}$ , meaning that although at the end of the path we have indeed rendered $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ and $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FA}}$ much shorter than $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D708}}$ and $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D705}}$ , we have not reached the expected Earth values for their ratios. However, we argue that $Pm$ and $E$ are small enough for thermal and viscous diffusivities to effectively become irrelevant, as is the case in Earth’s core. In contrast, it is a far more important result that the path achieves a correct value of the Earth magnetic Ekman number $E/Pm$ , since this number involves $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D6FA}}$ and $\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$ , two time scales which have a major impact on the MAC balance and on energy dissipation. Other dependences on $\unicode[STIX]{x1D716}$ for $E$ and $Pm$ than those in (3.22), (3.24) may besides be chosen, as long as $E/Pm$ still varies like $\sqrt{\unicode[STIX]{x1D716}}$ , and as long as $Pm$ decreases more slowly than $E^{3/4}$ along the path, an empirical condition (Christensen & Aubert Reference Christensen and Aubert2006) for maintaining self-sustained dynamo action that is also satisfied by (3.22), (3.24). For instance, perhaps a more elegant (but more computationally demanding) approach could be to simply define the path from the sole requirement to match the Earth’s core values of $Ra_{F}^{\ast }$ , $E$ and $Pm$ at $\unicode[STIX]{x1D716}=10^{-7}$ , i.e. $Ra_{F}^{\ast }=\unicode[STIX]{x1D716}(Ra_{F}^{\ast })_{0}$ , $E=\unicode[STIX]{x1D716}^{1.4}E_{0}$ and $Pm=\unicode[STIX]{x1D716}^{0.9}Pm_{0}$ . Along this alternative path (which satisfies all the above requirements) the constancy of $Rm$ , $\unicode[STIX]{x1D6EC}$ , together with the gradual enforcement of the MAC balance should then presumably emerge as results rather than prescriptions.

With the path input parameters now fully defined, we turn to the expected scalings for the outputs along the path. We have already seen that our path implies

The CE-type models in table 2 provide a first consistency check of our theoretical approach, as $Rm$ and $\unicode[STIX]{x1D6EC}$ are indeed constant to within $\pm 10\,\%$ for models chosen along the path down to $\unicode[STIX]{x1D716}=3\times 10^{-4}$ , and, by construction, obviously match the expected Earth values for $\unicode[STIX]{x1D716}=10^{-7}$ . This also applies to the ST models, though with a slightly different baseline for $\unicode[STIX]{x1D6EC}$ . Diffusivity-free scalings for other outputs can be derived using our path definitions (3.20), (3.22), (3.24) together with the relationships (3.18), (3.19), (3.25):

Here $Ro_{0}=1.11\times 10^{-2}$ , $\unicode[STIX]{x1D706}_{0}=1.60\times 10^{-2}$ and $A_{0}=0.69$ are respectively the Rossby, Lehnert and Alfvén numbers of model 0 (tables 1 and 2). Extrapolations of the scalings (3.27)–(3.29) to Earth’s core conditions ( $\unicode[STIX]{x1D716}=10^{-7}$ ) yield $Ro=3.5\times 10^{-6}$ , $\unicode[STIX]{x1D706}=2.8\times 10^{-4}$ and $A=1.2\times 10^{-2}$ . These are strikingly close to the estimates $Ro=3\times 10^{-6}$ , $\unicode[STIX]{x1D706}=2.1\times 10^{-4}$ and $A=1.5\times 10^{-2}$ obtained using the Earth values introduced above for $U,B,\unicode[STIX]{x1D70C},\unicode[STIX]{x1D6FA},D$ . Stating the result in another way, these scalings together with the output of model 0 provide independent predictions $U=6\times 10^{-4}~\text{m}~\text{s}^{-1}$ and $B=5.3~\text{mT}$ which compare very favourably with current estimates obtained through geophysical methods (Buffett Reference Buffett2010; Gillet et al. Reference Gillet, Jault, Canet and Fournier2010; Aubert Reference Aubert2014). While there is a rough equipartition of the kinetic and magnetic energy in model 0 ( $A_{0}^{2}=0.48$ ), the predicted separation increases with decreasing $\unicode[STIX]{x1D716}$ and is large at Earth’s core conditions, i.e. $A^{2}=O(10^{-4})$ .

In summary, here we have used a spatially invariant approximation of the D13 theory to establish a smooth (i.e. devoid of abrupt physical transitions) theoretical connection between classical numerical dynamos such as model 0 and the conditions of Earth’s core, along a uni-dimensional parameter space path constrained by the MAC balance and the need to maintain a constant value of the magnetic Reynolds number. While the extrapolations performed above are certainly supportive of this connection, the spatial invariance assumed in the path theory implies scaling exponents for the main outputs (3.25)–(3.29) differing from the values predicted by the (in principle) asymptotically correct D13 theory along the path. From equations (3.6)–(3.9), (3.20)–(3.24), we indeed get the D13 scalings

We note that in the case of $d_{min}$ and $d_{\bot }$ , the dependences in $\unicode[STIX]{x1D716}$ are indeed marginal and support spatial invariance to the extent that there is less than an order of magnitude variation between $\unicode[STIX]{x1D716}=1$ and $\unicode[STIX]{x1D716}=10^{-7}$ . The same also holds for $Rm$ , but interestingly not for $\unicode[STIX]{x1D6EC}$ . Both the D13 and spatially invariant path theories will be checked in § 4.3 against numerical data covering a wide portion of the path, leading to an assessment of the quality of the spatially invariant approximation. As can be expected, we will show that the direct numerical simulation results support the D13 set of exponents (3.30)–(3.35), while the results of large-eddy simulations to be introduced in § 3.4 are best described by the set of exponents (3.25)–(3.29), since these simulations assume spatial invariance to some extent.

3.3 Path theory in the ST set-up

The formulation of the path theory is straightforward in the CE set-up, because the convective power is in fact an input as it relates to the imposed mass anomaly flux $Ra_{F}^{\ast }$ (figure 1). The situation changes in the more widely studied ST set-up, where the power is an output. Starting from model 1 (table 2, model parameters subscripted with 1 in the following), which operates with the ST set-up in the same physical regime as model 0 with the CE set-up, we now wish to define the path parameter $\unicode[STIX]{x1D716}$ in terms of the power $p$ (which is then fully equivalent to what we did in the CE set-up, equation (3.20)):

We wish to also keep the same dependences on $\unicode[STIX]{x1D716}$ for $E$ , $Pm$ and $Pr$ as those introduced in (3.22)–(3.24), i.e. $E=\unicode[STIX]{x1D716}E_{1}$ , $Pm=\sqrt{\unicode[STIX]{x1D716}}~Pm_{1}$ and $Pr=1$ . To obtain the path rule for the input parameter $Ra^{\ast }$ , we recall the proportionality between the convective power and the mass anomaly flux (figure 1), which in the ST set-up is (Christensen & Aubert Reference Christensen and Aubert2006):

We thus need a relationship between the Nusselt and Rayleigh numbers. Such a relationship usually involves the canonical Rayleigh number $Ra$ and its critical value for convection onset $Ra_{c}\sim E^{-4/3}$ (Busse Reference Busse1970). It can generically be expressed as $Nu-1\sim (Ra/Ra_{c})^{\unicode[STIX]{x1D6FD}}\sim Ra^{\unicode[STIX]{x1D6FD}}E^{4\unicode[STIX]{x1D6FD}/3}$ , and with equation (2.10) and the fact that $Pr=1$ this yields $Nu-1\sim (Ra^{\ast })^{\unicode[STIX]{x1D6FD}}E^{-2\unicode[STIX]{x1D6FD}/3}$ . The convective power then scales with $p\sim (Ra^{\ast })^{\unicode[STIX]{x1D6FD}+1}E^{-2\unicode[STIX]{x1D6FD}/3+1}$ , and the requirements $p\sim \unicode[STIX]{x1D716}$ and $E\sim \unicode[STIX]{x1D716}$ thus imply $Ra^{\ast }=\unicode[STIX]{x1D716}^{2\unicode[STIX]{x1D6FD}/3(\unicode[STIX]{x1D6FD}+1)}~Ra_{1}^{\ast }$ . Rotating convection and dynamo studies suggest values ranging from $\unicode[STIX]{x1D6FD}=6/5$ (Christensen & Aubert Reference Christensen and Aubert2006; Aurnou Reference Aurnou2007; King et al. Reference King, Soderlund, Christensen, Wicht and Aurnou2010) to $\unicode[STIX]{x1D6FD}=3/2$ (the diffusivity-free prediction, Gillet & Jones Reference Gillet and Jones2006; Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). These yield very similar dependences $\unicode[STIX]{x1D716}^{4/11}$ to $\unicode[STIX]{x1D716}^{2/5}$ for $Ra^{\ast }$ . In the following we choose to retain the simpler expression

for completing the definition of the ST path. We note that the canonical Rayleigh number $Ra$ then increases with $\unicode[STIX]{x1D716}^{-8/5}$ , which is steeper than the increase of the critical value $Ra_{c}\sim \unicode[STIX]{x1D716}^{-4/3}$ , meaning that supercriticality indeed increases as we progress along the path.

Direct numerical simulations (table 2, ST/D cases) show that the obtained $\unicode[STIX]{x1D716}$ , reported in table 2 as the ratio $p/p_{1}$ , very closely matches the intended $\unicode[STIX]{x1D716}$ defined as $E/E_{1}$ . This provides a posterior check for our analysis. For the large-eddy simulations (ST/L cases) to be introduced in § 3.4, however, the obtained value is lower than the intended value, and the difference increases with decreasing $\unicode[STIX]{x1D716}$ . The hyperdiffusive treatment of the velocity and density anomaly fields indeed leads to a loss of the convective power carried through the interaction of small-scale velocity and density anomaly. Note though that this problem is specific to the ST set-up, since in the CE set-up the relationships $p/p_{0}\approx Ra_{F}^{\ast }/(Ra_{F}^{\ast })_{0}=E/E_{0}=\unicode[STIX]{x1D716}$ hold with markedly better accuracy (table 2).

3.4 Large-eddy simulations

Exhibiting numerical solutions for low values of $\unicode[STIX]{x1D716}$ is difficult because hydrodynamic turbulence sets up as one progresses along the path, as witnessed by the Reynolds number scaling that can be derived from (3.27):

with $Re_{0}=372$ in the CE case. As a result, the exploration of the path with a fully resolved direct numerical simulation (DNS) is currently possible only down to $\unicode[STIX]{x1D716}=0.1$ (cases CE/D and ST/D in table 2). Achieving lower values of $\unicode[STIX]{x1D716}$ requires to define relevant large-eddy simulations (LES, cases CE/L and ST/L in table 2). We have seen in § 3.1 that the essential part of the energy transfers should occur between the scales $d_{\bot }$ and $d_{min}$ , with $d_{\bot }>d_{min}$ and the separation between the two being mainly controlled by the level of magnetic turbulence $Rm$ (see (3.8) and (3.9)). We have also seen that we could approximate $d_{min}$ by a constant along a path of constant $Rm$ , implying that $d_{\bot }$ is also a constant. These points suggest that the hydrodynamic turbulence gradually setting up at scales below $d_{min}$ as we progress along this path is irrelevant to the determination of the large-scale structure. This latter structure is indeed controlled by a large-scale MAC balance at scale $d_{\bot }$ and only a small part of the convective power does cascade further below the scale $d_{min}$ , the major part being dissipated there through ohmic losses.

We may thus parameterise the hydrodynamic turbulence below scale $d_{min}$ without too much loss in physical relevance and accuracy, and our numerical computation may be restricted to the sole length scale range $d_{min}\leqslant d\leqslant D$ , thus alleviating the common scale-separation problems encountered at extreme control parameter values. Using $d_{min}\approx \sqrt{1/\,2Rm}$ (Christensen & Tilgner (Reference Christensen and Tilgner2004) a relationship also well verified in our data set, see table 2) and $Rm=1000$ , this corresponds to a numerical expansion up to spherical harmonic degree $l_{max}=140$ . In CE/L and ST/L large-eddy cases (table 2) we have used either $l_{max}=133$ or $l_{max}=170$ for reasons of numerical efficiency. We adopt a hyperdiffusive treatment recently proposed by Nataf & Schaeffer (Reference Nataf and Schaeffer2015) that we apply to the diffusion of momentum and density anomaly, but not to that of the magnetic field, which remains fully resolved. In the numerical implementation, the principle is to use effective diffusivities $\unicode[STIX]{x1D708}_{eff},\unicode[STIX]{x1D705}_{eff}$ that depend on the harmonic degree $l$ and the molecular diffusivities $\unicode[STIX]{x1D708},\unicode[STIX]{x1D705}$ according to

Here $l_{h}$ is the cutoff degree below which the hyperdiffusive treatment is not applied. For this treatment to remain physically accurate, it should satisfy the following requirements: (i) it should not perturb the large-scale MAC balance, i.e. $l_{h}$ should be sufficiently larger than $l_{\bot }=\unicode[STIX]{x03C0}D/d_{\bot }$ , and (ii) the additional viscous energy losses that it implies should remain small relatively to ohmic losses. These requirements are satisfied here by prescribing $l_{h}=30$ , a value which, as we will see in § 4.1, is three times larger than our estimate $l_{\bot }\approx 10$ , and values of $q_{h}$ close to 1 (table 2), enabling a very smooth increase of hyperdiffusivity with harmonic degree. In § 4 we will validate this treatment by comparing DNS and LES simulations in the range where both are feasible, and by demonstrating that the additional viscous losses are indeed negligible in the asymptotic limit. There is some leeway in the choices of $q_{h}$ and $l_{h}$ . For the CE/L case performed at $\unicode[STIX]{x1D716}=0.1$ we have indeed checked that varying $q_{h}$ in the range $1.02\leqslant q_{h}\leqslant 1.07$ yields similar results to within 5 %. The same holds when varying $l_{h}$ in the range $15\leqslant l_{h}\leqslant 60$ in the CE/L case performed at $\unicode[STIX]{x1D716}=0.33$ . Given the numerical schemes used to solve the equations, the hyperdiffusive treatment only applies to the lateral directions of the numerical calculation, making it still necessary to increase the radial resolution as $\unicode[STIX]{x1D716}$ decreases. This is however a far more tractable numerical problem than that of expanding the grid in all three spatial directions. Combined with the fact that we avoid treating viscous boundary layers, this permitted to calculate CE cases down to $\unicode[STIX]{x1D716}=3.33\times 10^{-4}$ and an Ekman number $E=10^{-8}$ . This is the lowest Ekman number reached to date in self-sustained spherical convective dynamos, though it should be acknowledged that this value applies to the degree range $l\leqslant 30$ only, after which the effective Ekman number increases.

4.1 Spatial structure of the asymptotic MAC force balance

In figure 2 we first analyse the force balance operating in the ST simulations (similar results hold for the CE simulations). To this end, for each force $\boldsymbol{f}$ we follow Soderlund et al. (Reference Soderlund, King and Aurnou2012, Reference Soderlund, Sheyko, King and Aurnou2015) and compute

Figure 2.

Force balance spectrum in the ST numerical simulations (viscous boundary layers excluded), presented as a function of the spherical harmonic degree $l$ and normalised relative to the peak amplitude of the pressure gradient, for two DNS (a,b) and two LES (c,d) cases performed at varying $\unicode[STIX]{x1D716}$ along the path. (b) and (d) present a comparison of DNS and LES at a constant $\unicode[STIX]{x1D716}$ value. Note that the maximal spherical harmonic degree $l_{max}$ of the calculations differs in the four cases (see table 2). The solid lines represent the time-averaged force balance, and the shaded regions represent the $\pm 1$ standard deviation in the instantaneous force balance. The harmonic degree $l_{\bot }$ is defined as the intersection between the buoyancy and Lorentz forces, such that these two forces equally contribute to balance the part of the Coriolis force not balanced by the pressure gradient. For comparison purposes, the predictions $\unicode[STIX]{x03C0}/d_{\bot }$ obtained from equation (4.2) using the values in table 2 are also reported in (a–d).

There are a few novelties of our calculation relative to Soderlund et al. (Reference Soderlund, King and Aurnou2012, Reference Soderlund, Sheyko, King and Aurnou2015). First, we have access to the part of the Coriolis force not balanced by the pressure gradient (the ageostrophic part). Second, we reveal the actual force balance in the bulk of the fluid by excluding viscous boundary layers from the calculation. Note that this methodology is similar to that used in Yadav et al. (Reference Yadav, Gastine, Christensen, Wolk and Poppenhaeger2016b ) and the results are hence directly comparable. A third point which is original to the present work is that we decompose $f_{rms}^{2}$ into a sum of contributions $f_{l}^{2}$ along spherical harmonic degrees and present $f_{l}$ as a function of $l$ . At the order following that of the diagnostic balance between Coriolis and pressure forces, a MAC balance is robustly observed in all simulations between the ageostrophic Coriolis, buoyancy and Lorentz forces. This force balance is structurally similar among all simulations (compare figure 2 a–d). The harmonic degree $l_{\bot }\approx 10$ corresponding to an optimal MAC balance may be identified as the intersection between the contributions of the Lorentz and buoyancy forces, and is largely invariant, both in the DNS and LES simulations. For harmonic degrees smaller than $l_{\bot }$ (length scales larger than $d_{\bot }$ ) the main balance is between the buoyancy and ageostrophic Coriolis forces, and the Lorentz force substitutes to the buoyancy force for harmonic degrees larger than $l_{\bot }$ (length scales smaller than $d_{\bot }$ ). The temporal variability of this basic structure is quite low, implying that the MAC balance holds at all instants in the numerical simulations. Independent predictions $l_{\bot }=\unicode[STIX]{x03C0}/d_{\bot }$ (reported in figure 2) can be obtained from equation (3.4) and the values in table 2, where the dimensionless D13 length scale $d_{\bot }$ is

Note that (4.2) reintroduces the contribution of $f_{ohm}$ that should arise from (3.3) but has been neglected in (3.4). These predictions are essentially invariant and closely match the measure obtained by the crossings of forces in figure 2. This confirms that $l_{\bot }$ indeed corresponds to the D13 scale $d_{\bot }$ . This also demonstrates the conjecture that the large scale at which the convective dynamo is organised is indeed controlled by the vorticity equivalence (3.4) that follows from the MAC balance. It is interesting to note that the starting point of the ST path (model 1, figure 2 a) already exhibits a well-respected MAC balance, though the contributions of inertial and viscous forces are still sizeable. Decreasing $\unicode[STIX]{x1D716}$ by a factor 10 in a DNS simulation (figure 2 b) decreases the contribution of inertia by a factor 3 approximately, as expected from the predicted Rossby number scaling (3.27) along the path. The gap between inertial and viscous forces also enlarges, though at a pace somewhat slower than that predicted by Reynolds number scaling (3.39), hinting at a non-trivial length scale for viscous dissipation which we leave outside the scope of the present study.

The agreement between DNS and LES performed at a similar value of $\unicode[STIX]{x1D716}$ (figure 2 b,d) is excellent regarding the harmonic structure and relative amplitudes of pressure, Coriolis, buoyancy, Lorentz and inertial forces up to degree $l\leqslant 70$ . The LES obviously shows an increase in viscous forces after degree $l_{h}=30$ but it is worth noting that up to $l=l_{max}$ these do not reach a first-order dominance, such that the nature of the small-scale force balance is also preserved. An interesting point is that for degrees $l\leqslant l_{h}$ , the viscous force also slightly increases in the LES relatively to the DNS. Finally, comparing the DNS case with $\unicode[STIX]{x1D716}=1$ to the LES case with $\unicode[STIX]{x1D716}=6.42\times 10^{-3}$ (figure 2 a,c), we see that despite the use of hyperdiffusivity, this latter simulation has reached a nearly inviscid regime in the harmonic degree range $l\leqslant 50$ , with the viscous force at least two orders of magnitude lower than the MAC forces. In this harmonic degree range, inertia is also rendered largely subdominant, with its typical amplitude one to two orders of magnitude below the MAC forces. All these points underline the success of the path and LES approaches in reaching the MAC balance. This regime has been observed down to the lowest value of $\unicode[STIX]{x1D716}$ that we could achieve, halfway in logarithmic distance between the classical models and the Earth conditions along the parameter path. These results provide strong support to the claims that the MAC balance is indeed the planetary, asymptotic regime, and that our numerical models are close to this regime.

Figure 3.

Time-averaged magnetic (a,c) and kinetic (b,d) energy spectra in the CE (top row) and ST (bottom row) set-ups, for DNS (a,b) and LES (c,d) simulations performed at varying $\unicode[STIX]{x1D716}$ . To obtain the normalised magnetic energy spectra in (a,c), the magnetic field amplitude is normalised by its scaling prediction $\unicode[STIX]{x1D706}_{0,1}\unicode[STIX]{x1D716}^{1/4}$ (3.28), where $\unicode[STIX]{x1D706}_{0,1}$ stands for the Lehnert number of model 0 (CE) or 1 (ST). Similarly, in (b,d) the amplitude of the velocity field is also normalised by $Ro_{0,1}\unicode[STIX]{x1D716}^{1/2}$ (3.27). The total normalised energy is close to 1 in all cases. For reference, the dashed vertical lines mark the harmonic degree $l_{\bot }$ corresponding to $d_{\bot }$ as identified in figure 2.

Figure 4.

Comparison of time-averaged energy spectra obtained in DNS and LES at a given $\unicode[STIX]{x1D716}$ value corresponding to an Ekman number $E=3\times 10^{-6}$ , for CE (a) and ST (b) cases. The amplitude of magnetic and velocity fields are both normalised by $\unicode[STIX]{x1D706}_{0,1}\unicode[STIX]{x1D716}^{1/4}$ . The total normalised magnetic energy is then close to 1 while the total normalised kinetic energy is close to the separation factor $A^{2}=0.2$ (CE) and $A^{2}=0.14$ (ST).

4.2 Spatial invariance of the solutions along the path

We now turn to the analysis of the spatial structure of the fields along the path. The large-scale invariance can be first demonstrated by examining energy spectra from direct numerical simulations (figure 3 a,b) performed in the range $0.1\leqslant \unicode[STIX]{x1D716}\leqslant 1$ . The magnetic spectral energy density profile (figure 3 a) is structurally invariant in the range $l\leqslant l_{diss}$ , where $l_{diss}\approx 133$ is the harmonic degree corresponding to the magnetic dissipation length scale $d_{min}$ . The spectral profile of kinetic energy density (figure 3 b) is also invariant in the large-scale range $l\leqslant 30$ but shows an enrichment in smaller scales as $\unicode[STIX]{x1D716}$ decreases. Large-eddy simulations (figure 3 c,d) present even higher levels of structural similarity, both for the velocity and magnetic fields, throughout the resolved spectral range. CE-type models are generally more invariant than ST-type models, owing to their better control on $\unicode[STIX]{x1D716}$ . Comparing a DNS and a LES performed at a similar value of $\unicode[STIX]{x1D716}$ (figure 4) it is clear that despite using a significantly lower value of $l_{max}$ , the LES accurately captures the spectral distribution of magnetic energy almost up to degree $l_{max}=133$ , and the kinetic energy distribution at least up to degree 30. These results are found to hold irrespectively of the boundary condition choice, with CE and ST yielding comparable energy distributions, at the exception of velocity for $l=1$ which is enhanced in the CE set-up by an inhomogeneous mass anomaly forcing at the inner boundary (Aubert et al. Reference Aubert, Finlay and Fournier2013). In all cases, the kinetic energy spectra feature a peak around the value $l_{\bot }\approx 10$ already identified in the force balance (figure 2), implying a predominant organisation of convective structures around this scale.

Figure 5.

Evolution of velocity planforms (snapshots, units of $\unicode[STIX]{x1D702}/D$ or magnetic Reynolds number units) in CE-type calculation as $\unicode[STIX]{x1D716}$ decreases (a–c), and comparisons at fixed $\unicode[STIX]{x1D716}$ between DNS and LES (b) and between CE and ST cases (d). In (a) and (c), the left and right panels respectively present snapshots of the radial velocity in a half equatorial plane (red is outwards), and a meridional cut of the azimuthal velocity (blue is westwards). In (b) and (d) all panels represent snapshots of the equatorial radial velocity.

To check how these spectral results translate into the physical space, we further examine equatorial and meridional velocity planforms in the CE set-up (figure 5 a–c). From the kinetic energy spectra observed in figure 3, we expect the equatorial planforms of radial velocity to remain structurally similar as $\unicode[STIX]{x1D716}$ varies, particularly in the large-scale range where an azimuthal wavenumber $m\approx l_{\bot }\approx 10$ should emerge. Comparing the CE-type DNS simulations at $\unicode[STIX]{x1D716}=1$ and $\unicode[STIX]{x1D716}=0.1$ (figure 5 a,b, left panels), we indeed observe the preservation of structures at this large scale, while a decrease in $\unicode[STIX]{x1D716}$ leads to the refinement of small-scale details. Compared to its LES counterpart, the DNS at $\unicode[STIX]{x1D716}=0.1$ (figure 5 b) also refines the small-scale details while not changing the large-scale structure, thus visually confirming the relevance of the LES. Comparing equatorial planforms obtained in the CE and ST DNS set-ups (figure 5 d), we find structural similarity close to the inner boundary and a richer small-scale content in the ST set-up close to the outer boundary (as previously shown by Sakuraba & Roberts Reference Sakuraba and Roberts2009). This reflects a fundamental change in how convective power is transported. In the ST set-up, the part of the convective power carried by these small-scale velocity structures together with their corresponding density anomalies (see figure 6 e) is not accounted for in the LES (compare the values of $p$ in table 2), leading to deviations from the ideal path (see § 3.3) and to a decrease of the LES accuracy. In contrast, CE-type LES simulations preserve the DNS convective power remarkably well (compare again $p$ values in table 2), save for a slight change of mixing properties (see trend for $\unicode[STIX]{x1D6FE}=p/Ra_{F}^{\ast }$ in figure 1 b) that should remain negligible at Earth’s core conditions. This confirms that convective power is mostly transferred at large scales in the CE set-up.

Turning now to meridional planforms, the axially columnar structure of the azimuthal velocity field (figure 5 a,c, right panels) outside the tangent cylinder is also preserved as $\unicode[STIX]{x1D716}$ decreases, together with retrograde (westward) polar vortices at the upper and lower ends of the tangent cylinder, confirming the ansatz $d_{/\!/}\sim D$ made in § 3. The power input in the system indeed increases as $\unicode[STIX]{x1D716}$ decreases, but so does also the rotational constraint, leading to a preservation of the columnar structure. Types CE and ST have zonal flows of similar amplitude despite the change in boundary conditions from stress free to rigid (not shown, previously already documented by Yadav et al. Reference Yadav, Gastine and Christensen2013a ). Livermore, Bailey & Hollerbach (Reference Livermore, Bailey and Hollerbach2016) have suggested that the amplitude of zonal flows may asymptotically scale differently depending on mechanical boundary conditions if these are limited by the residual viscosity. This does not apply here because zonal flows in spherical, convective, dipole-dominated dynamos such as those discussed here are thermal-wind limited (Aubert Reference Aubert2005).

Figure 6.

Evolution of equatorial density anomaly planforms (half-plane snapshots, orange denotes lighter fluid) as $\unicode[STIX]{x1D716}$ decreases (a–d) and comparisons at fixed $\unicode[STIX]{x1D716}$ between DNS and LES (b,c) and between CE and ST cases (b,e). In CE cases (a–d) the dimensionless density anomaly $C$ is normalised by $C_{0}\unicode[STIX]{x1D716}^{-1/2}$ , with $C_{0}=6000$ . In the ST case (e) $C$ , is normalised by $Pr/E$ with $Pr=1$ and $E=3\times 10^{-6}$ .

Equatorial planforms of the density anomaly (figure 6) also lead to similar observations, with broadly structurally similar patterns at large scales (figure 6 a–d) and the DNS simply refining the picture obtained with the LES (figure 6 b,c). Both CE and ST set-ups produce dominant bottom-originated convective plumes (figure 6 b,e), as theoretically expected, but the ST set-up additionally produces small-scale, secondary return plumes originated at the outer boundary. In the CE set-up, the dimensionless density anomaly $C$ is found to scale with $C_{0}\unicode[STIX]{x1D716}^{-1/2}$ , as can be predicted from the expression of convective power (2.18), together with the velocity scaling (3.27) and the fact that the ratio $p/Ra_{F}^{\ast }$ is constant (figure 1), and we find $C_{0}=6000$ . The typical dimensional density anomaly $\unicode[STIX]{x1D6FF}C$ in Earth’s core may then be estimated through the following diffusivity-free, power-driven scaling obtained from equations (2.11) and (3.20):

with $C_{0}(Ra_{F}^{\ast })_{0}=0.162$ . Using $\unicode[STIX]{x1D716}=10^{-7}$ together with the previously introduced values $D=2260~\text{km}$ , $\unicode[STIX]{x1D6FA}=7.29\times 10^{-5}~\text{s}^{-1}$ and $g_{0}=10~\text{m}~\text{s}^{-2}$ , this yields $\unicode[STIX]{x1D6FF}C/\unicode[STIX]{x1D70C}=6\times 10^{-8}$ . Converting this anomaly into a temperature deviation $\unicode[STIX]{x1D6FF}T$ from the core adiabat, we obtain $\unicode[STIX]{x1D6FF}T=\unicode[STIX]{x1D6FF}C/\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D70C}=6\times 10^{-3}~\text{K}$ , where $\unicode[STIX]{x1D6FC}=10^{-5}~\text{K}^{-1}$ is the thermal expansion coefficient. Either $\unicode[STIX]{x1D6FF}C$ or $\unicode[STIX]{x1D6FF}T$ are in good agreement with geophysical estimates (e.g. Aurnou et al. Reference Aurnou, Andreadis, Zhu and Olson2003; Christensen & Aubert Reference Christensen and Aubert2006), which again underlines the relevance of the path chosen to connect model and Earth’s core conditions.

Figure 7.

Evolution of the radial magnetic field pattern at the model outer boundary (snapshots, units of $\sqrt{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D707}\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FA}}$ or Elsasser units, orange is outwards) in the CE set-up as $\unicode[STIX]{x1D716}$ decreases (a–c), and comparison at fixed $\unicode[STIX]{x1D716}$ between DNS and LES (b,d).

Figure 8.

Radial magnetic field pattern at the model outer boundary (snapshots, units of $\sqrt{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D707}\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FA}}$ or Elsasser units, orange is outwards), filtered below spherical harmonic degree and order 13, in the CE set-up for a DNS at $\unicode[STIX]{x1D716}=1$ (a), and a LES at $\unicode[STIX]{x1D716}=10^{-3}$ (b).

Figures 7 and 8 present maps of the radial magnetic field at the model outer boundary for CE-type calculations. The full resolution magnetic field patterns (figure 7) are again remarkably invariant, with an axial dipole-dominated morphology, high-latitude flux patches of normal polarity at the intersection of the outer boundary with the tangent cylinder, westward-drifting, low-latitude patches of normal polarity, a magnetic flux deficit inside the tangent cylinder due to the action of polar vortices and localised inverse flux patches. The good results obtained in the original CE dynamo model (similar to model 0, figure 7(a), see also Aubert et al. Reference Aubert, Finlay and Fournier2013) as regards the static and kinematic similarity of the model output to the geomagnetic field are thus preserved, and demonstrated here to pertain to the asymptotic regime. Figure 8 further shows that when filtered to the typical resolution $l\leqslant 13$ of modern geomagnetic observations, the model output preserves the Earth-like morphological properties (as objectively defined in Christensen et al. Reference Christensen, Aubert and Hulot2010) of the initial model 0. We also note in figures 7 and 8 that the field amplitude presented in Elsasser units remains constant as $\unicode[STIX]{x1D716}$ decreases, in agreement with (3.26). Finally, considering the ratio of the root-mean-squared magnetic field inside the shell to that at the outer boundary, the CE models yield a constant value of approximately 5 if the full outer boundary magnetic field is considered, or 7 if that field is truncated to spherical harmonic degree 13. This latter value is in agreement with geophysical estimates (e.g. Christensen & Aubert Reference Christensen and Aubert2006; Buffett Reference Buffett2010; Gillet et al. Reference Gillet, Jault, Canet and Fournier2010).

We finish our structural analysis by examining radial profiles of the kinetic and magnetic energies (figure 9) at a fixed $\unicode[STIX]{x1D716}$ value, for DNS and LES simulations carried out in the CE and ST set-ups. These profiles show that the choice of boundary conditions has a weak influence on the distribution of energy in the shell. Indeed, the ST and CE kinetic energy profiles (figure 9 a) are strikingly similar, with the only noticeable differences residing in the ST viscous boundary layers that are absent in the CE set-up. Likewise, the ST and CE magnetic energy profiles (figure 9 b) only differ in their amplitude, as already noted in table 2 (see also figure 11 b). Turning now to the agreement between DNS and LES, the best results are obtained in the CE set-up, while in the ST set-up the discrepancies are most visible near the outer boundary. There, the inhibition of the small-scale radial motion yields a quieter zone where the magnetic energy can concentrate.

Figure 9.

Spherical- and time-averaged radial distributions of the kinetic and magnetic energies $r^{2}u^{2}$ (a) and $r^{2}B^{2}$ (b), for CE models (black) computed at $\unicode[STIX]{x1D716}=0.1$ , and ST models (grey) computed at $\unicode[STIX]{x1D716}=0.1$ (DNS) and $\unicode[STIX]{x1D716}=0.0891$ (LES). DNS (solid lines) and LES (dashed lines) calculations are presented. Velocity and magnetic fields are presented in units of $\sqrt{\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FA}}$ and $\sqrt{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D707}\unicode[STIX]{x1D702}\unicode[STIX]{x1D6FA}}$ (Elsasser units), respectively, such that the radial integral of the magnetic energy density is the Elsasser number $\unicode[STIX]{x1D6EC}$ and the radial integral of the kinetic energy density is $A^{2}\unicode[STIX]{x1D6EC}$ , revealing the energy separation factor $A^{2}$ .

4.3 Evolution of diagnostics and asymptotic scalings along the path

Figure 10.

Evolution of the Rossby number $Ro$ (a), Reynolds $Re$ and magnetic Reynolds $Rm$ numbers (b) with $\unicode[STIX]{x1D716}$ .

Figure 11.

Evolution of the Lehnert number $\unicode[STIX]{x1D706}$ (a) and Elsasser number $\unicode[STIX]{x1D6EC}$ (b) with $\unicode[STIX]{x1D716}$ . In (a) the Lehnert number is corrected by the factor $f_{ohm}^{1/2}$ in order to account for the variations of the ohmic dissipation fraction within the scaling theory (e.g. Christensen & Aubert Reference Christensen and Aubert2006; Davidson Reference Davidson2013). These variations are significant in the DNS cases but unimportant in the LES cases (see figure 12 a). In the CE-type LES case indeed, the corrected best fit presented in (a) is closely similar to the uncorrected fit, which yields $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D716}^{0.25}$ . In (b), shaded regions represent the $\pm 1$ standard deviation of temporal fluctuations relative to the time average.

Figure 10 presents the evolution of the velocity diagnostics $Ro,Re,Rm$ with  $\unicode[STIX]{x1D716}$ . DNS and LES diagnostics are here again in broad agreement. Results obtained with type CE and ST collapse on each other once presented as functions of $\unicode[STIX]{x1D716}$ , illustrating the convergence between the two types of boundary conditions. CE and ST-type DNS cases respectively follow $Ro=Ro_{0}\unicode[STIX]{x1D716}^{0.44\pm 0.02}$ , $Ro=Ro_{1}\unicode[STIX]{x1D716}^{0.45\pm 0.01}$ and CE-type LES cases follow $Ro=Ro_{0}\unicode[STIX]{x1D716}^{0.49\pm 0.01}$ (figure 10 a). These best fit exponents are in excellent agreement with the D13 prediction 4/9 from equation (3.30) and the path theory prediction 1/2 from scaling (3.27), respectively. This is fully expected since the DNS leaves the spatial structure of the solution free while the LES constrains it in the lateral directions. The slope difference between the two theories is resolvable here thanks to the extremely low scatter of the numerical data and to the 3 decades available in LES cases. The path approach is thus demonstrably advantageous when compared to systematic samplings of the parameter space (e.g. Christensen & Aubert Reference Christensen and Aubert2006). The best fit exponent 0.46 obtained for $Ro$ in the ST-type LES simulations is significantly below the predicted value 1/2, because the $\unicode[STIX]{x1D716}$ values of these simulations are misplaced relative to the values that they should have along the ideal path, as we have seen in § 3.3 and table 2. It is thus generally not advisable to attempt an interpretation of ST-type LES scaling exponents. Finally, the magnetic Reynolds number (figure 10 b) is confirmed to be roughly constant, as already seen in table 2.

Figure 11 presents the evolution of the magnetic field diagnostics $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D6EC}$ with $\unicode[STIX]{x1D716}$ . CE and ST-type calculations show more difference than in figure 10, but still produce very comparable results. DNS and LES diagnostics are again in good agreement. Within the standard error 0.01 on exponent determination, CE and ST-type DNS cases (figure 11 a) are in close agreement with the ohmic fraction-corrected scaling  $\unicode[STIX]{x1D706}/f_{ohm}^{1/2}\sim \unicode[STIX]{x1D716}^{1/3}$ (see 3.31) expected from D13 and Christensen & Aubert (Reference Christensen and Aubert2006). LES calculations performed in the CE set-up yield the corrected Lehnert number scaling $\unicode[STIX]{x1D706}/f_{ohm}^{1/2}\sim \unicode[STIX]{x1D716}^{0.27\pm 0.01}$ , and the uncorrected scaling $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{0}\unicode[STIX]{x1D716}^{0.25\pm 0.01}$ , both also in close agreement with the prediction (3.28) from the path theory. We do not attempt to scale the LES simulations performed in the ST set-up, for the same reason as above. The Elsasser number $\unicode[STIX]{x1D6EC}$ is confirmed to be roughly constant (figure 11 b) in CE simulations. The CE-type DNS cases thus appear to invalidate the D13 scaling (3.33), but this is the consequence of the persistent increase in $f_{ohm}$ not being accounted for in (3.33). ST simulations have a higher $\unicode[STIX]{x1D6EC}$ baseline than CE simulations, and the LES cases show a residual decrease with decreasing $\unicode[STIX]{x1D716}$ , which can be explained by noting that their rotation rate increases too rapidly with respect to the injected power, as witnessed by the misplaced values of $\unicode[STIX]{x1D716}$ relative to the ideal path (table 2).

We next turn to the analysis of the dissipation diagnostics $f_{ohm}$ and $d_{min}$ in our simulations (figure 12). The ohmic dissipation fraction $f_{ohm}$ (figure 12 a) is largely similar in CE and ST cases. LES simulations produce $f_{ohm}$ values somewhat below those obtained from DNS simulations performed at similar values of $\unicode[STIX]{x1D716}$ , an expected consequence of increased viscous losses due to hyperdiffusivity. However, as $\unicode[STIX]{x1D716}$ decreases these additional losses are dominated by the increasing ohmic losses, such that LES simulations gradually approach the expected asymptotic behaviour $f_{ohm}=1$ . This shows that along the chosen path, the additional viscous losses indeed become asymptotically negligible, thus validating our hyperdiffusive treatment of turbulence. The magnetic dissipation length scale $d_{min}$ is confirmed to be weakly variable within the investigated $\unicode[STIX]{x1D716}$ range (note the narrow linear scale used for the ordinate axis of figure 12 b). A closer analysis shows that the DNS simulations in ST and CE set-ups support a decrease $d_{min}\sim \unicode[STIX]{x1D716}^{0.09\pm 0.01}$ , matching the dependence $\unicode[STIX]{x1D716}^{1/12}$ predicted by the D13 theory (3.34). LES simulations in the ST set-up do not show a systematic evolution of $d_{min}$ with $\unicode[STIX]{x1D716}$ . LES simulations in the CE set-up show a weak variation $d_{min}\sim \unicode[STIX]{x1D716}^{0.022\pm 0.002}$ . This can be expected because the LES constrain the length scale to be constant in the two lateral directions, leaving only the radial length scale free to evolve. The discrepancy between DNS and LES, or between the D13 and path theories should remain quite low at Earth’s core asymptotic conditions. For $\unicode[STIX]{x1D716}=10^{-7}$ we indeed predict $d_{min}\approx 0.006$ (corresponding to approximately 14 km) if we follow the DNS trend, and $d_{min}\approx 0.015$ (or 34 km) following the LES trend.

Figure 12.

Evolution of the ohmic dissipation fraction $f_{ohm}$ (a) and the magnetic dissipation length $d_{min}$ (b) with $\unicode[STIX]{x1D716}$ . Shaded regions represent the $\pm 1$ standard deviation of temporal fluctuations relative to the time average.

Figure 13.

Evolution of the Alfvén number $A$ (a), and of the level of Taylor constraint enforcement ${\mathcal{T}}$ (b) with $\unicode[STIX]{x1D716}$ . The Alfvén number in (a) is the square root of the kinetic to magnetic energy ratio. In (b), the Taylor constraint enforcement level on cylindrical surfaces (large symbols) is as defined in (2.21). For reference, the enforcement level on spherical surfaces (where ${\mathcal{T}}$ should be unconstrained) is also presented (small symbols). In (b), shaded regions represent the $\pm 1$ standard deviation of temporal fluctuations relative to the time average.

Finally, we investigate in figure 13 the level of enforcement of strong-field dynamo action and of the Taylor constraint in our simulations. Figure 13(a) shows that at low values of $\unicode[STIX]{x1D716}$ , the solutions have a kinetic energy much smaller than the magnetic energy, as witnessed by the low values of the squared Alfvén number $A^{2}$ measuring their ratio. Together with values of the Elsasser number $\unicode[STIX]{x1D6EC}$ in excess of 1 (figure 11 b), this clearly characterises a regime of strong-field dynamo action. The power-law exponents obtained for $A$ in the DNS simulations range between 0.17 and 0.19, in consistency with the results obtained in figures 10 and 11. LES calculations performed in the CE set-up yield $A=A_{0}\unicode[STIX]{x1D716}^{0.24\pm 0.01}$ , in close agreement with the prediction (3.29) from the path theory. Analysing the variations of ${\mathcal{T}}$ in CE-type simulations (figure 13 b), we find that the level of Taylor constraint enforcement increases (i.e. ${\mathcal{T}}$ decreases) steadily with decreasing $\unicode[STIX]{x1D716}$ , both in DNS and LES simulations. This again confirms that the hyperdiffusive treatment performed in the LES simulations introduces dynamically negligible additional amounts of the viscous force in asymptotic conditions. Respectively to the time-averaged values of ${\mathcal{T}}$ , we find low instantaneous deviations of $\pm 20\,\%$ , meaning that the Taylor constraint is enforced at all times in addition to holding everywhere in the fluid domain. As theoretically expected, cylindrical reference surfaces are particular for the Taylor constraint in the sense that the integral Lorentz force remains unconstrained on other surfaces such as spheres. Qualitatively, the enforcement of the Taylor constraint links with strong-field dynamo action and the low value of the Alfvén number $A=\unicode[STIX]{x1D70F}_{A}/\unicode[STIX]{x1D70F}_{U}$ . This ratio indeed also determines how many torsional oscillation periods $\unicode[STIX]{x1D70F}_{A}$ can occur within a convective overturn time $\unicode[STIX]{x1D70F}_{U}$ to minimise (by virtue of the Lenz law) the Taylor state deviations that trigger these oscillations. More quantitatively, the level of Taylor constraint enforcement should scale like the ratio of inertia at the Alfvén time scale and magnetic forces at scale $d_{\bot }$ , leading to (in dimensional form)

This suggests that ${\mathcal{T}}$ should primarily scale with $A$ , as evidenced by the similarity between the power-law exponents obtained in figures 13(a,b). We thus propose a scaling law of the form ${\mathcal{T}}={\mathcal{T}}_{0}\,Ad_{\bot }/A_{0}(d_{\bot })_{0}$ . Using $A_{0}=0.69$ , ${\mathcal{T}}_{0}\approx 0.2$ and $d_{\bot }/(d_{\bot })_{0}=\unicode[STIX]{x1D716}^{1/9}$ (3.35), at Earth’s core conditions i.e. $\unicode[STIX]{x1D716}=10^{-7}$ and $A=1.2\times 10^{-2}$ we predict

such that the Taylor constraint should be strongly but not completely enforced, to leave the possibility for torsional oscillations to exist.