2. Governing equations

Thermal convection and magnetic field generation in the Earth’s core are modelled in the present study using the most classical set of equations. The rotating incompressible MHD equations are coupled to the heat equation under the Boussinesq approximation. Convection is driven by an imposed temperature difference across a spherical shell (of inner radius $r_{i}$ and outer radius $r_{o}$ ). Magnetic field generation by dynamo action requires a flow with an appropriate geometry and sufficient amplitude, which can be achieved if the control parameter (measuring the efficiency of the thermal driving) is increased away from the onset of convection. The parameter space for such dynamos has been extensively studied by Christensen and collaborators (e.g. Christensen et al. Reference Christensen, Olson and Glatzmaier1999), providing a detailed description of the ‘phase diagram’ for dynamo action in this set-up (i.e. the region in the parameter space for which different dynamo solutions are produced). The governing equations are solved in a spherical shell ( $r_{i}/r_{o}=0.35$ ) and in a rotating reference frame. The reference frame is such that the velocity vanishes on both spheres (no-slip boundaries), a temperature difference is maintained for all time across the shell, and both the inner and the outer domains are assumed electrically insulating. The equations governing the velocity $\boldsymbol{u}$ , magnetic field $\boldsymbol{B}$ and temperature field $T$ are then, in non-dimensional form,

In the above equations $L=r_{o}-r_{i}$ has been used as length scale, ${\it\tau}_{{\it\eta}}=L^{2}/{\it\eta}$ as time scale, and $({\it\Omega}{\it\mu}_{0}{\it\rho}_{0}{\it\eta})^{1/2}$ as magnetic field scale. This non-dimensional form is well suited for strong-field dynamos (e.g. Fearn Reference Fearn1998). The following non-dimensional numbers have been introduced: the Ekman number $\mathit{E}={\it\nu}/{\it\Omega}L^{2}$ , the magnetic Prandtl number $\mathit{Pm}={\it\nu}/{\it\eta}$ , the Roberts number $q={\it\kappa}/{\it\eta}$ , and the Rayleigh number $Ra={\it\alpha}g{\rm\Delta}TL/{\it\kappa}{\it\Omega}$ , where $g=g_{o}/r_{o}$ with $g_{o}$ the gravity at $r=r_{o}$ (note that the Rayleigh number is here modified from its most standard definition to account for the stabilising effect of rotation). It is also useful to define the Prandtl number $Pr={\it\nu}/{\it\kappa}\equiv Pm/q$ . In the present work, $\mathit{E}$ is set to $3\times 10^{-4}$ and $Pr$ to unity; thus it follows that $q=\mathit{Pm}$ in the sequel.

These equations are numerically integrated using the Parody code, originally developed by the author and improved with several collaborators (see Dormy Reference Dormy1997; Schrinner et al. Reference Schrinner, Petitdemange, Dormy and Schrinner2012). The numerical resolution in the simulations reported here is 132 grid points in radius, with spherical harmonic decomposition of degrees up to $\ell _{max}=256$ and modes up to $m_{max}=64$ . The models were integrated for up to 10 magnetic diffusion times. In order to ensure the validity of the new solutions presented here, these simulations were also kindly reproduced by V. Morin using the Magic code, developed by Glatzmaier and modified by Christensen and J. Wicht. Both codes have been validated through an international benchmark (see Christensen et al. Reference Christensen, Aubert, Cardin, Dormy, Gibbons, Glatzmaier, Honkura, Jones and Kono2001).

3. Weak- and strong-field dynamos

Following the same approach as Morin & Dormy (Reference Morin and Dormy2009), we study the bifurcation from the purely hydrodynamic solution to the dynamo state using the Rayleigh number as control parameter. At fixed $E$ and $Pm$ , the Rayleigh number needs to exceed a given value for a dynamo solution (non-vanishing field) to exist. We report here direct numerical simulations performed at large values of the magnetic Prandtl number $Pm$ . One may object that such a parameter regime is irrelevant to dynamo action in liquid metals (characterised by a small magnetic Prandtl number). We will however argue that considering large values of $Pm$ can compensate for the excessive role of inertial terms in numerical dynamo models, and is a necessary consequence of the large values assumed by the Ekman number.

Figure 1 presents the bifurcation diagrams obtained for $Pm=q=14$ and $Pm=q=18$ . The magnetic field strength, as measured by the classical Elsasser number ${\it\Lambda}=B^{2}/(2{\it\Omega}{\it\rho}{\it\mu}{\it\eta})$ is represented versus the Rayleigh number, normalised by its value at the onset of thermal convection $\mathit{Ra}_{c}$ (here $\mathit{Ra}_{c}=60.8$ ). Each point on this figure corresponds to a time-averaged fully three-dimensional simulation. The time variability of the dynamo mode is reported using the standard deviation.

Figure 1.

Weak- and strong-field branches for $\mathit{E}=3\times 10^{-4}$ and for (a) $\mathit{Pm}=14$ and (b)  $\mathit{Pm}=18$ . Symbols indicate the time-averaged Elsasser number (●, stable; ○, unstable); the time variability of the dynamo mode is reported using the standard deviation (indicated with error bars). The insets present the same graphs in lin–log scale, so that the weak branch is more clearly visible.

Figure 1 is characterised by a supercritical bifurcation (as reported in Morin & Dormy (Reference Morin and Dormy2009) for their ‘large’ values of the magnetic Prandtl number, i.e. $\mathit{Pm}=6$ ). However, this first dynamo branch rapidly destabilises to a second branch of much stronger amplitude. This strong-field branch can be maintained for decreasing values of the Rayleigh number. The magnetic field is dominated by the axial dipole on both branches (see figure 2). The strong-field branch on figure 1(a) is hysteretic to the onset of the dynamo itself: once on this branch, the control parameter can be decreased below the critical value for dynamo bifurcation, while maintaining a dipolar magnetic field.

Figure 2.

Radial component of the magnetic field produced at the surface of the model for $\mathit{Pm}=18$ and $\mathit{Ra}/\mathit{Ra}_{c}=1.73$ on the weak-field (a) and strong-field (b) branches. The field is characterised by a strong axial dipolar component on both branches.

This new branch completes the sequence of bifurcation diagrams introduced in Morin & Dormy (Reference Morin and Dormy2009), and the complete three-dimensional bifurcation diagram (including the results of Morin & Dormy (Reference Morin and Dormy2009)) for $\mathit{E}=3\times 10^{-4}$ and $\mathit{Pr}=1$ is reported on figure 3 versus $\mathit{Ra}/\mathit{Ra}_{c}$ and $\mathit{Pm}$ . The corresponding two-dimensional bifurcation diagrams, for lower values of $\mathit{Pm}$ , are available in Morin & Dormy (Reference Morin and Dormy2009). On such a three-dimensional diagram, the dynamo bifurcation can also be envisaged at fixed value of $\mathit{Ra}/\mathit{Ra}_{c}$ and varying $\mathit{Pm}$ .

Figure 3.

Three-dimensional bifurcation diagram for a fixed Ekman number $\mathit{E}=3\times 10^{-4}$ . Solid lines mark linear interpolation between realised direct numerical simulations. Dashed lines offer a plausible interpretation of the unstable branches.

Figure 3 demonstrates how the transition between the different types of bifurcation takes place for different values of $Pm$ . The study of Morin & Dormy (Reference Morin and Dormy2009) indicates that, as $E$ is decreased (in the moderate range numerically achievable), the overall bifurcation diagram remains largely unaltered but shifted towards lower values of $Pm$ and larger $\mathit{Ra}/\mathit{Ra}_{c}$ .

Transitions between these two branches of dynamo solutions are obtained by varying only slightly the control parameter at the edge of a given branch. This can produce either a runaway field growth (figure 4 a), or a catastrophic collapse (figure 4 b) of the magnetic field. The time at which the forcing (as measured by the Rayleigh number) has been modified (by less than 3 % in each case) is indicated by an arrow on each graph.

Figure 4.

Transition between the weak- and strong-field branches for $\mathit{E}=3\times 10^{-4}$ and $\mathit{Pm}=18$ (same as figure 1 b). Runaway growth from the weak- to the strong-field branch as $\mathit{Ra}/\mathit{Ra}_{c}$ is increased (arrow) from 1.78 to 1.83 (a) and catastrophic decay to the weak-field branch as it is decreased (arrow) from 1.73 to 1.68 (b). The lower branch (left part of (a) and right part of (b) corresponds to small, but non-vanishing, magnetic fields.

No significant changes on the typical length scale of the flow can be reported by comparing the weak- and strong-field branches. This is probably due to the fact that the viscous length scale is not very small at the value of the Ekman number considered here ( $3\times 10^{-4}$ ). Smaller values of the Ekman number are undoubtedly needed if one is to appreciate a change in the typical length scale of the flow. One can note however that the Nusselt number is 80 % larger on the strong-field branch than on the weak-field branch.

4. Force balance

The bifurcation diagram presented on figure 1 is reminiscent of a longstanding theoretical expectation originally introduced as the ‘weak’- and ‘strong’-field branches (see Roberts Reference Roberts, Roberts and Soward1978; Roberts & Soward Reference Roberts and Soward1992). The existence of these two branches in the limit of vanishing viscous forces was introduced through the investigation of magneto-convection studies (for reviews, see Fearn, Roberts & Soward Reference Fearn, Roberts, Soward, Galdi and Straughan1986; Proctor Reference Proctor, Proctor and Gilberts1994).

Soward (Reference Soward1979) investigated the onset of magneto-convection in the cylindrical annulus configuration with sloping boundaries. He found that in most cases the critical Rayleigh number first starts to increase with the Elsasser number, until ${\it\Lambda}\sim O(\mathit{E}^{1/3})$ , before decreasing. This pointed to the probable existence of a weak-field branch, and the occurrence of a turning point marking the end of the weak-field branch when ${\it\Lambda}\sim O(\mathit{E}^{1/3})$ . Simultaneously, Fearn (Reference Fearn1979a ,Reference Fearn b ) performed a similar study in the spherical geometry. There again, the Rayleigh number for the thermal Rossby mode may first increase with increasing Elsasser number, yet it eventually decreases to reach a minimum for ${\it\Lambda}\sim O(1)$ . A more recent study of magneto-convection (Jones et al. Reference Jones, Mussa and Worland2003) focused on the ‘weak-field’ regime and confirmed its existence.

The above asymptotic scenario assumes a small value of both $E$ and $Pm$ , whereas direct numerical simulations in the self-excited dynamo regime require overestimated values of both numbers.

In order to test the above ideas in the numerical simulations, we need to investigate the dominant force balance relevant to these dynamo modes. Figure 5 presents an instantaneous cross-section of the zonal velocity on both branches for a given parameter set ( $\mathit{Pm}=18$ , $\mathit{Ra}/\mathit{Ra}_{c}=1.73$ ). The contour intervals are equally spaced between the minimum and the maximum value for each figure. The zonal flow is nearly three times larger on figure 5(a), so that the contour intervals are not identical on the two plots. On the one hand, the weak-field branch saturates while the zonal flow remains essentially geostrophic; the flow is characterised by quasi-geostrophic convection columns. The zonal flow on the strong-field branch, on the other hand, strongly departs from bidimensionality, demonstrating that the rapid rotation constraint has been relaxed. On figure 5(b), a localised jet appears near the equator, which marks a clear departure from geostrophy. The flow is in general less anisotropic along the direction of the axis of rotation. If this corresponds to the weak-field versus strong-field branches as introduced by Roberts, it implies that the Lorentz forced has achieved a balance with the Coriolis term and thus relaxed the rapid rotation constraint. This dominant balance is usually referred to as the ‘magnetostrophic balance’.

Figure 5.

Azimuthal velocity in a meridional cross-section for an arbitrary time on the weak-field branch (a) and on the strong-field branch (b) for the same parameter set ( $\mathit{Pm}=18$ , $\mathit{Ra}/\mathit{Ra}_{c}=1.73$ ).

The Elsasser number ${\it\Lambda}$ was introduced to measure an order of magnitude of the relative strength of the Lorentz force with respect to the Coriolis force. It achieves this aim remarkably well in asymptotic studies, for the huge distinction between the strong-field balance, characterised by ${\it\Lambda}\sim O(1)$ , and the weak-field branch, characterised by ${\it\Lambda}\sim O(\mathit{E}^{1/3})$ . In numerical works, however, as small parameters (such as the Ekman number) are not asymptotically small, the measure provided by this non-dimensional number is then not accurate enough. Finer estimates of this force balance can then be constructed. Introducing $\{\cdot \}$ as an ‘order-of-magnitude’ operator, we can write

where $U$ is a typical, say root mean square (r.m.s.), value for the velocity field, $B$ a typical value for the magnetic field, and $\ell _{B}$ the typical magnetic dissipation length scale (see also Oruba & Dormy Reference Oruba and Dormy2014a ).

The classical definition of the Elsasser number is obtained by assuming $\ell _{B}\sim L$ and a statistical balance between induction and diffusion of the magnetic field,

which yields $U\sim {\it\eta}/L$ . Then (4.1) provides the standard expression for the Elsasser number ${\it\Lambda}=B^{2}/(2{\it\Omega}{\it\rho}{\it\mu}{\it\eta})$ . This expression provides a sensible description of the force balance for asymptotic studies, yet finer estimates appear to be needed for numerical studies.

One can note, for example, that $U\sim {\it\eta}/L$ amounts to assuming $Rm\sim O(1)$ . A finer description of the force balance (4.1) can be obtained by estimating $U$ via $\mathit{Rm}\,{\it\eta}/L$ . Inserting this definition in (4.1) yields

Table 1 presents a comparison of the classical Elsasser number ${\it\Lambda}$ and the modified Elsasser number ${\it\Lambda}^{\prime }$ on both branches. The magnetic Reynolds number $Rm$ is here defined on the r.m.s. velocity, and the typical magnetic dissipation length scale $\ell _{B}$ is defined, as in Oruba & Dormy (Reference Oruba and Dormy2014a ), as

Figure 6 presents the time variation of the modified Elsasser number on the weak-field branch (dashed) and on the strong-field branch (solid line) for the same parameter set, $\mathit{Pm}=18$ , $\mathit{Ra}/\mathit{Ra}_{c}=1.73$ .

Figure 6.

Fluctuation of the modified Elsasser number on both branches for the same parameter set ( $\mathit{Pm}=18$ , $\mathit{Ra}/\mathit{Ra}_{c}=1.73$ ).

The modified Elsasser number, offering a finer description of the force balance, reveals that the Lorentz force is significantly weaker than the Coriolis force on the weak-field branch and that the two terms are indeed of comparable amplitude on the strong-field branch.

Table 1.

Typical estimates of the Elsasser number and the modified Elsasser number ( ${\it\Lambda}^{\prime }={\it\Lambda}L/(\mathit{Rm}\ell _{B})$ , see text) on the weak- and strong-field branches. The modified Elsasser number offers a finer measurement of the force balance.

The orders of magnitude derived above indicate that the anticipated balance between the Coriolis and Lorentz forces is plausible. To achieve a finer validation than simple orders of magnitude, we can assess whether the two terms tend to balance each other locally in space. To this aim, one can consider the curl of the momentum equation (2.1), neglecting both the inertial term and the viscous term,

If we now consider the radial component of the above equation, i.e. its toroidal component, the first term on the right-hand side disappears. The remaining two terms were computed numerically at a given instant in time and on a cross-section in an arbitrary meridional plane. These quantities are presented on figure 7.

Figure 7.

Meridional cross-sections on the strong-field branch at a given time for the same parameters and same phase as in figure 5(b). On the left side is presented $-2\partial \boldsymbol{u}/\partial z\boldsymbol{\cdot }\boldsymbol{e}_{r}$ , and on the right side $(\boldsymbol{{\rm\nabla}}\times ((\boldsymbol{{\rm\nabla}}\times \boldsymbol{B})\times \boldsymbol{B}))\boldsymbol{\cdot }\boldsymbol{e}_{r}$ using the same colour range.

Deviations between the two cross-sections can imply only non-vanishing inertial and/or viscous effects. Estimations of these terms reveals that the viscous term accounts for the differences visible on the figures (inertia being one order of magnitude smaller). The comparison reveals such effects (in particular, in viscous boundary layers), but otherwise clearly demonstrates that the radial component of the curl of the Lorentz force balances that of the Coriolis force, as expected in the strong-field limit.

The viscous force will of course not always be negligible in the parameter regime considered here. It can be more important at some places or times. To illustrate this, the quantities represented on figure 7 are represented at a later time in the form of three-dimensional isosurfaces on figure 8. The blue and red isosurfaces, respectively, correspond to $\pm 90\,\%$ of the peak values. While deviations from magnetostrophy are obvious, in particular comparing the centre of each figure, the dominant magnetostrophic balance is highlighted. Deviations are here primarily due to viscous forces. Boundary layers have not been represented in these figures.

Figure 8.

Three-dimensional isosurfaces of the radial component of the curl of the Coriolis force (a) and the Lorentz force (b).

Figures 7 and 8 highlight the balance between the non-gradient part of the Lorentz and Coriolis terms. Differences are primarily due to viscous effects and remain small (but non-vanishing) on average. The r.m.s. value of the sum of the two quantities plotted in these figures, averaged in time, exceeds by a factor 5 that of their difference.