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

HAL
open science

# Toward scale-separated weak-field spherical dynamos

R. Teed, E. Dormy

► To cite this version:

R. Teed, E. Dormy. Toward scale-separated weak-field spherical dynamos. Physical Review Fluids, 2025, 10 (10), pp.103702. 10.1103/9r2f-t4wk. hal-05407951

HAL Id: hal-05407951

https://hal.science/hal-05407951v1

Submitted on 9 Dec 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
Toward scale-separated weak-field spherical dynamos

R. J. Teed*

School of Mathematics and Statistics, University of Glasgow,
University Place, Glasgow G12 8SQ, United Kingdom

E. Dormy†

Département de Mathématiques et Applications, UMR-8553, École Normale Supérieure,
CNRS, PSL University, 75005 Paris, France

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

(Received 13 June 2025; accepted 4 September 2025; published 15 October 2025)

Recent numerical experiments of dynamo action relevant to the generation of the geomagnetic field have produced different regime branches identified within bifurcation diagrams. Notable are separate branches in which the resultant magnetic field is either weak or strong. Weak-field solutions can be identified by the prominent role of viscosity on the motion whereas the magnetic field has a leading order effect on the flow in strong-field solutions. For a given Ekman number, $\mathcal{E}$ (measuring the ratio of viscosity to rotational effects), the existence of these branches and bistability between them is reliant on a small enough magnetic Ekman number, $\mathcal{E}_m$ (measuring the ratio of magnetic diffusion to rotational effects, so $\mathcal{E}/\mathcal{E}_m = P_m$, the magnetic Prandtl number). Both branches are known to produce large scale dipolar magnetic fields but do not exhibit an expected scale separation between the flow and magnetic field. In this work, by reducing $\mathcal{E}_m$, we identify a variety of dynamo states on the weak-field branch beyond the known dipolar solutions. Specifically, hemispherical and nondipolar dynamos were found, in addition to the usual dipolar solutions. Some solutions exhibit clear scale separation between small-scale flow and large-scale magnetic field, despite the large ratio of viscosity to magnetic diffusion. Numerical solutions in this regime have not been observed before and they offer a first connection with earlier theoretical work based on mean-field theory.

DOI: 10.1103/9r2f-t4wk

# I. INTRODUCTION

Magnetic fields generated by dynamo action are a common feature of astronomical bodies throughout the solar system and the universe beyond. The dynamo process involves the conversion of kinetic to magnetic energy with buoyancy typically providing the energy source. For instance, Earth's magnetic field is constantly being replenished through complex motions within its electrically conducting molten iron outer core [1]. Such motions are a result of thermocompositional convection driven by the release of heat or light elements at the boundary with the solid inner core.

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*Contact author: Robert.Teed@glasgow.ac.uk</footnote>

<footnote>†Contact author: Emmanuel.Dormy@ens.fr</footnote>
A longstanding problem in dynamo theory is an explanation of how small scale turbulent convective motions can generate large scale ordered magnetic field (e.g., Earth's dipolar magnetic field). In the case of Earth, unmagnetized convective motions are expected to be less than a few kilometers in size [2] whereas the magnetic field is able to operate on the scale of the planet itself. This observation opens the way for an analytical treatment of the dynamo problem based on scale separation. This approach, known as "mean-field" theory was pioneered by Ref. [3] and introduces two scales thereby splitting the flow and field into "mean" and "fluctuating" parts. For a review of mean-field theory, see Refs. [4–6].

As computation power has expanded, direct numerical solution (DNS) of the governing equations has become an increasingly vital tool in the understanding of the problem. Through these endeavors, a variety of dynamo regimes—dipolar, multipolar, hemispherical, weak-field, strong-field, etc.—have been identified by exploring the parameter space [7–10]. Much focus has been given to regimes that are supposedly “Earth-like,” e.g., Refs. [11,12] as well as to models tailored to specific other planets, e.g., Ref. [13]. Some studies have aimed to connect direct numerical simulations with the mean-field theory [14,15] but have not shown the necessary scale-separation. Recent years have seen a focus on analyzing the magnitude and hierarchy of the individual forces (the “force balance”) that appear in the governing equations [16–19]. This can be used a measure of a given dynamo regime’s appropriateness for modeling the Geodynamo (or other planetary dynamos).

Weak- and strong-field branches of dynamo action (identified by the size of the ratio of rotational to magnetic forces), based on the effect of varying the magnetic field strength in magnetoconvection studies [20–24], were hypothesized by Ref. [25]. In this bifurcation configuration, viscously dominated weak-field solutions appear first, as the onset of dynamo action, and this is likely how the Geodynamo first manifested. At this weak level of convective driving, the viscous force allows for some violation of the Taylor-Proudman constraint (which otherwise demands motion is invariant in the direction of the rotation axis leading to no convection and no dynamo action). These weak solutions become unstable at larger driving leading to runaway growth of the magnetic field onto a branch of strong-field solutions; for a review of these ideas see Refs. [1,5]. The separate weak- and strong-field dipolar branches and bistability between them can be found using DNS [10,26], provided the magnetic Ekman number, $E_m$, is small enough. Equivalently, this can instead be posed as requiring the ratio of the viscosity and magnetic diffusivity (i.e., the magnetic Prandtl number, $P_m$) to be large enough. These parameters will be defined in the following section. Alternatively, at larger $E_m$ and higher driving, bistability can also be found between weak-field dipolar solutions and a multipolar branch [27,28] as inertia gradually replaces viscosity as the primary force enabling violation of the Taylor-Proudman constraint. The multipolar branch is a distinct weak-field solution in which inertia, but not magnetic forces, are important. The focus of our current study will be on the viscous weak-field branch and its transformation to the strong-field branch. For this reason, in what follows it should be understood that “weak-field” refers to the viscous weak-field branch.

Since weak-field solutions are located close to the onset of (nonmagnetic) convection, they are often dominated by the tall, thin, columnar flows of linear/weakly nonlinear nonmagnetic convection; see, e.g., Ref. [29]. In keeping with those solutions, the forces of the weak-field branch therefore exhibit a dominant viscous-Archimedean-Coriolis (VAC) balance. The existence of a generated weak magnetic field does not change the generic columnar structure of the flow pattern but it is strong enough to affect the preferred mode and drift rate of convection, compared to the nonmagnetic case [30]. This is in stark contrast to the strong-field regime where flows become more complicated and three-dimensional with a dominant magnetic-Archimedean-Coriolis (MAC) balance.

The Geodynamo is expected to lie within the strong-field regime although the weak-field regime remains of interest both from a theoretical point-of-view and as a means to address the onset of the Geodynamo. The aim of this work is to investigate the weak-field branch at moderate rotation rates where its dependence on the magnetic Ekman number can be explored at relatively low computational expense. In addition to known dipolar solutions, we find a range of nondipolar solutions exist on the weak-field branch, including hemispherical and quadrupole-dominated models.
Weak-field solutions give way to the strong-field branch at large enough driving. Interestingly, several simulations demonstrate a clear separation between a small lengthscale in the flow and a large lengthscale in the magnetic field. This is a new kind of numerical dynamo that establishes a first connection between simulations and classical mean-field dynamo theory.

Our methods are discussed in Sec. II, followed by the results of several simulations in Sec. III and a discussion in Sec. IV.

## II. METHODS

We aim to numerically solve equations governing the evolution of the flow and magnetic field in a geometrical setup relevant to the present-day Earth's core. As such, we consider a spherical shell rotating about the vertical coordinate with rotation vector $\Omega = \Omega e_z$. Using spherical coordinates, $(r, \theta, \phi)$, the shell is confined between an inner sphere of radius, $r_i$, and outer sphere of radius, $r_o > r_i$ such that the aspect ratio $\chi = r_i/r_o = 0.35$, as appropriate for the present-day Earth. Within the shell resides a Boussinesq, electrically conducting fluid of constant density, $\rho_0$, with intrinsic diffusivity parameters of the fluid (the kinematic viscosity, $\nu$, thermal diffusivity, $\kappa$, and magnetic diffusivity, $\eta$), each assumed to be constant. Gravity acts radially inwards such that $\mathbf{g} = -g_0 \mathbf{r}$ where $\mathbf{r} = (1 - \chi) r e_\nu$. A temperature difference of $\Delta T$ is maintained between the boundaries, allowing differential heating to drive convection.

The coupled set of partial differential equations governing the evolution of velocity, $\mathbf{u}$, pressure, $p$, temperature, $T$, and magnetic induction, $\mathbf{B}$, are

$$ \mathcal{E}_m \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = -\nabla p - 2 \mathbf{e}_z \times \mathbf{u} + (\nabla \times \mathbf{B}) \times \mathbf{B} + \tilde{R} a T \mathbf{r} + \mathcal{E} \nabla^2 \mathbf{u}, \quad (1a) $$

$$ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T = q \nabla^2 T, \quad (1b) $$

$$ \frac{\partial \mathbf{B}}{\partial t} - \nabla \times (\mathbf{u} \times \mathbf{B}) = \nabla^2 \mathbf{B}, \quad (1c) $$

$$ \nabla \cdot \mathbf{u} = 0, \quad (1d) $$

$$ \nabla \cdot \mathbf{B} = 0, \quad (1e) $$

where we have nondimensionalized using lengthscale, $d = r_o - r_i$, timescale, $d^2/\eta$, temperature scale, $\Delta T$, and magnetic scale $\sqrt{\rho_0\mu_0\eta\Omega}$. The nondimensional parameters appearing in our set of equations are the magnetic Ekman number, $\mathcal{E}_m$, modified Rayleigh number, $\tilde{R}a$, Ekman number, $\mathcal{E}$, Roberts number, $q$, defined as

$$ \mathcal{E}_m = \frac{\eta}{\Omega d^2}, \quad (2a) $$

$$ \tilde{R}a = \frac{\alpha g_0 \Delta T d}{\Omega \eta}, \quad (2b) $$

$$ \mathcal{E} = \frac{\nu}{\Omega d^2}, \quad (2c) $$

$$ q = \frac{\kappa}{\eta}, \quad (2d) $$

where $\alpha$ is the coefficient of thermal expansion. The Rayleigh number we use is a modified version of the (rotational) Rayleigh number, $Ra = \alpha g_0 \Delta T d / \Omega \kappa$, which itself is an updated version of the classical (nonrotational) Rayleigh number $R = \alpha g_0 \Delta T d^3 / \nu \kappa$. Hence the various definitions are related by $\tilde{Ra} = qRa = q\mathcal{E}R$. For convenience, we also introduce $\tilde{Ra}' = Ra / Ra_c = \tilde{Ra} / qRa_c$ where $Ra_c$ is the value of $Ra$ for the onset of nonmagnetic convection (for a given $\mathcal{E}$). This allows us to measure the level of supercriticality (above convective onset) of each simulation by $\tilde{Ra}'$. We note that equations (1a)–(1e) are commonly posed alternatively, using a different set of input
parameters. These could typically involve the (thermal) Prandtl number, $P_r \equiv \mathcal{E}/q\mathcal{E}_m = \nu/\kappa$ and the magnetic Prandtl number, $P_m \equiv \mathcal{E}/\mathcal{E}_m = \nu/\eta$, as inputs, rather than $\mathcal{E}_m$ and $q$. We will thus also make reference to $P_m$ in our results section. We use impenetrable, rigid, electrically insulating, isothermal boundary conditions. We numerically solve the governing equations (1a)–(1e) (subject to our boundary conditions), using the Leeds spherical dynamo code [31,32]. This code has been used for nearly 20 years to perform a wide range of studies. Data produced as part of this study are available [33].

Output parameters we discuss in our results include the total kinetic and magnetic energy, which under our nondimensionalization take the form

$$E_K = \frac{1}{2} \int \mathbf{u}^2 dV \quad \text{and} \quad E_M = \frac{1}{2\mathcal{E}_m} \int \mathbf{B}^2 dV.$$

By decomposing the energies into spherical harmonic degree, $l$, and order, $m$ (denoted $E_K^{(l,m)}$ and $E_M^{(l,m)}$), we form power spectra of the energy by summing over either $l$ or $m$: $(E_K)_l$, $(E_M)_l$, $(E_K)_m$, and $(E_M)_m$. As a measure of magnetic field strength we use the modified Elsasser number:

$$\Lambda' = \frac{B^2 d}{2\mu_0 \rho_0 \Omega \eta \mathcal{R}_m \ell_B},$$

where $\ell_B$ is the magnetic dissipation lengthscale and the magnetic Reynolds number, $\mathcal{R}_m = Ud/\eta$, is given by $\sqrt{2E_K}$ under our choice of nondimensionalization. We additionally define a scale dependent magnetic Reynolds number, $(\mathcal{R}_m)_l = \sqrt{2(E_K)_l \pi l}$. To measure the morphology of the magnetic field at the outer boundary we define the relative dipole moment by

$$f_{\text{dip}} = \left( \frac{E_M^{(1,0)}(r_o)}{\sum_{l=0}^{12} \sum_{m=0}^{l} E_M^{(l,m)}(r_o)} \right)^{1/2},$$

and form the Lowes-Mauersberger spectrum [34,35] using the Gauss coefficients ($g_l^m$ and $h_l^m$):

$$g_l = \sum_m (g_l^m)^2 + (h_l^m)^2.$$

For each given force, $\mathbf{F}$, appearing in equation (1a), we report solenoidal forces [19] and, in particular, the power spectrum in $l$ of compensated curls of forces, $\hat{C}_l = (l(l+1))^{-1/2}C_l$, where

$$C_l^2 = 2 \sum_{m=0}^{l} \int_{r_1+b}^{r_0-b} [(\nabla \times \mathbf{F})^2]_{lm} r^2 dr,$$

$b$ is the boundary layer thickness, and the prime denotes a halving of the $m=0$ term of the summation. See Ref. [19] for details of the idea of solenoidal forces and, in particular, Eq. (S.12) of its Supplemental Material for more detail of the definition of $\hat{C}_l$.

It has long been known [29,36,37] that thermal convection onsets in the form of columns of convection with a preferred azimuthal wavenumber and drifting as a thermal Rossby wave. Since several of the simulations we present are close to this onset regime, we report the azimuthal wavenumber, $\hat{m}$, and the (prograde) drift rate, $\omega$, of the columns. The former is simply the value of $m$ at the peak value of $(E_K)_m$, and the latter is calculated from the complex growth rate. The growth rate is itself found from the time-derivative of the leading spectral coefficient, $f_{\hat{m}\hat{m}}$ of a suitable measured quantity; see, e.g., Ref. [38]. We choose to use the radial velocity sampled at $r \approx \hat{r} = 0.55$ (i.e., just above the bottom boundary) and the drift rate is calculated from this. Additionally, we introduce a new parameter to measure the relative importance of the Rossby wave to the flow in the system. This is

$$W = \frac{\tau_{\omega}}{\tau_{U}} \sim \frac{2\pi/\omega}{\ell/U} \sim \frac{2\pi/\omega}{2\pi\hat{r}/mU} = \frac{mU}{\hat{r}\omega},$$
<figure><img src="image_4.png" /><figcaption>FIG. 1. (a) Modified Elsasser number, $\Lambda'$, and (b) the dipolarity, $f_{\text{dip}}$, as a function of supercriticality, $\tilde{Ra}'$. Each point represents a single simulation where data has been volume averaged and time-averaged. The following symbols are used for the given dynamo regime: filled/unfilled circles for strong/weak-field dipolar; triangles for weak-field hemispherical; squares for weak-field nondipolar; diamonds for fluctuating multipolar; crosses for failed dynamos (with data moved to the lower horizontal axis for visibility in $\tilde{Ra}'$ space).</figcaption></figure>

where $\tau_\omega$ is the time for a column of fluid to drift to its next neighbor's position and $\tau_U$ is the turnover time of a flow particle around a column of fluid. Additionally, $\ell$ is the horizontal lengthscale of a column of fluid and $U$ is the rms velocity. For systems dominated by waves (flows) we expect $W < 1$ ($W > 1$).

### III. RESULTS

We begin this section with the presentation of a suite of simulations discussed in this work. All simulations are performed with $\mathcal{E} = 10^{-4}$ and $q = P_m = \mathcal{E}/\mathcal{E}_m$ (or, equivalently, $P_t = 1$). We focus on three values of the magnetic Ekman number ($\mathcal{E}_m \in \{10^{-4}, 8.3 \times 10^{-6}, 4 \times 10^{-6}\}$) or, equivalently, $P_m \in \{1, 12, 25\}$) and vary the Rayleigh number. Since our focus is on solutions that emerge on distinct weak- and strong-field branches of dynamo action, and given computational limitations on the value of $\mathcal{E}$, the values of $P_m$ studied must be sufficiently large [1, 5, 10, 19] (also see discussion in Sec. IV). We focus on a single value of the Ekman number to allow for adequate coverage of the solution space in $\mathcal{E}_m$ space and $\tilde{Ra}'$ space. Note that simulations performed at the low values of $\mathcal{E}_m$ used here are numerically demanding.

Weak- and strong-field (or other) solutions can be tracked through the modified Elsasser number [Fig. 1(a)] and the (non)dipolar nature of the solution by $f_{\text{dip}}$ [Fig. 1(b)]. Weak-field and strong-field solutions exist for small enough $\mathcal{E}_m$ as evidenced by the large jump in the magnitude of $\Lambda'$ at $\tilde{Ra}' \approx 2.2$ for $\mathcal{E}_m = 8.3 \times 10^{-6}$ or $\tilde{Ra}' \approx 1.8$ for $\mathcal{E}_m = 4 \times 10^{-6}$. There exists bistability between the weak- and strong-field branches for $\mathcal{E}_m = 8.3 \times 10^{-6}$ (and very likely also for $\mathcal{E}_m = 4 \times 10^{-6}$, though we did not confirm this due to computational expense). For $\mathcal{E}_m = 10^{-4}$, separate weak- and strong-field branches are not apparent, and the solutions do not exhibit either weak- or strong-field behavior. Indeed, the Elsasser number never attains a strong-field value of $\Lambda' \sim \mathcal{O}(1)$.

The majority of our solutions are dipolar, indicated by the large values of $f_{\text{dip}}$. In particular, solutions on the strong-field branch for both $\mathcal{E}_m = 8.3 \times 10^{-6}$ and $\mathcal{E}_m = 4 \times 10^{-6}$ are dipolar as are those on the weak-field branch for $\mathcal{E}_m = 8.3 \times 10^{-6}$. Perhaps surprisingly, at $\mathcal{E}_m = 4 \times 10^{-6}$, nondipolar solutions are found on the weak-field branch. These solutions will be the focus of subsequent sections. The behavior at $\mathcal{E}_m = 10^{-4}$ is typical of solutions from much of the previous
<table><caption>TABLE I. Values of input and key output parameters for the simulations discussed in detail in Sec. III. Values of output parameters are time-averaged over the post-transient phase of each simulation and are given to three significant figures. A complete dataset including all simulations carried out for this work is also available [33].</caption><thead><tr><th>Solution</th><th>ε</th><th>ε<sub>m</sub></th><th>q</th><th>R̃a'</th><th>R<sub>m</sub></th><th>Λ</th><th>Λ'</th><th>f<sub>dip</sub></th><th>ℓ<sub>B</sub></th></tr></thead><tbody><tr><td>Weak-field dipolar</td><td>10<sup>-4</sup></td><td>8.3×10<sup>-6</sup></td><td>12</td><td>2.1</td><td>216</td><td>0.133</td><td>0.00615</td><td>0.834</td><td>0.101</td></tr><tr><td>Strong-field dipolar</td><td>10<sup>-4</sup></td><td>8.3×10<sup>-6</sup></td><td>12</td><td>5</td><td>456</td><td>24.4</td><td>1.10</td><td>0.857</td><td>0.0485</td></tr><tr><td>Weak-field hemispherical</td><td>10<sup>-4</sup></td><td>4×10<sup>-6</sup></td><td>25</td><td>1.4</td><td>193</td><td>0.243</td><td>0.0163</td><td>0.532</td><td>0.0772</td></tr><tr><td>Weak-field nondipolar</td><td>10<sup>-4</sup></td><td>4×10<sup>-6</sup></td><td>25</td><td>1.6</td><td>262</td><td>0.139</td><td>0.00738</td><td>0.282</td><td>0.0717</td></tr><tr><td>Strong-field dipolar</td><td>10<sup>-4</sup></td><td>4×10<sup>-6</sup></td><td>25</td><td>1.8</td><td>334</td><td>9.70</td><td>0.518</td><td>0.725</td><td>0.0561</td></tr></tbody></table>

literature, e.g., Ref. [7]: strongly dipole-dominated solutions at weak driving are replaced by fluctuating multipolar solutions for large enough $R̃a$.

In the following sections we begin with a brief review of the essential properties of the known weak- and strong-field dipolar solutions previously found in direct numerical simulations such as those we perform here (i.e., those at $ε_m = 8.3×10^{-6}$). We then discuss our key new findings at $ε_m = 4×10^{-6}$ where nondipolar solutions show clear lengthscale separation between the flow and magnetic field. Time-averaged output parameters for the simulations discussed in detail are shown in Table I.

### A. Weak- and strong-field solutions at $ε_m = 8.3×10^{-6}$

The theory of weak- and strong-field dipolar branches [25,39] and their identification in DNS, including bistability between branches [10,26], have been discussed previously. The bistability of solutions was first identified at $ε = 3×10^{-4}$ [10] and also, recently, at $ε = 10^{-4}$ [19,40]. Here we briefly discuss bifurcations of the dynamo onset at $ε_m = 8.3×10^{-6}$ (i.e., $P_m = 12$).

For these input parameters, onset of convection occurs at $Ra_c ≈ 69.7$ (where, by definition, $R̃a' = 1$), after which dynamo onset on the weak-field branch appears at $R̃a' ≈ 2$, before becoming unstable to the strong-field branch at larger $R̃a'$. Typical plots for $B_r$ of the two types of solutions are shown in Fig. 2. These solutions are similar to those conducted by Ref. [10] (at $ε = 3×10^{-4}$) and by Ref. [19] (at $ε = 10^{-4}$, albeit at different values of $R̃a$). The strongly symmetric nature of the weak-field solution can be seen in $B_r$, where a dipolar solution with 4-fold symmetry in azimuth is preferred. The strong-field solution has a noticeably more complicated field structure. Nevertheless

<figure><img src="image_5.png" /><figcaption>FIG. 2. Snapshots of $B_r$ on the outer boundary for (a) a typical weak-field branch solution ($R̃a' = 2.1$, $ε_m = 8.3×10^{-6}$ at $t ≈ 1.13$) and (b) a typical strong-field branch solution ($R̃a' = 5$, $ε_m = 8.3×10^{-6}$ at $t ≈ 0.23$).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 3. Plots for a hemispherical solution ($\tilde{Ra}' = 1.4$, $E_m = 4 \times 10^{-6}$). Snapshots at $t \simeq 3.61$ of (a) $B_r$ on a spherical surface at the outer boundary and (b) isocontours of $\omega_z$ colored with $u_z \omega_z$, the helicity associated with the columnar flow, blue negative and green positive.</figcaption></figure>

both solutions are predominantly dipolar and do not show clear separation between the scale of the flow and the scale of the magnetic field (see further discussion in the Appendix).

## B. Scale separation

We now reduce the magnetic Ekman number to $E_m = 4 \times 10^{-6}$ (or, equivalently, increase the magnetic Prandtl number to $P_m = 25$), retaining $E = 10^{-4}$. The dynamo regimes and branches at onset under these conditions is rather different to the simple transition from weak-field dipolar to strong-field dipolar seen in Sec. III A. At this reduced $E_m$, scale separation between the velocity and magnetic fields now becomes a key feature of (at least) two different regimes found in the viscously dominated (i.e., weak-field) solutions. At $E_m = 4 \times 10^{-6}$, we find dynamo onset at $\tilde{Ra}' \lesssim 1.4$; specifically, our first dynamo solution is found at $\tilde{Ra}' = 1.4$ whereas a run at $\tilde{Ra}' = 1.3$, yields a failed (decaying) dynamo. Several different regimes (discussed below) are found as $\tilde{Ra}$ is increased until we once more obtain a transition to a strong-field dipolar branch at $\tilde{Ra}' = 1.8$.

### 1. Hemispherical solutions

The dynamo solution found at lowest supercriticality ($\tilde{Ra}' = 1.4$) is hemispherical in nature [Fig. 3(a)]. The flow takes the form of convection columns attached to the lower boundary with $m = 7$ symmetry [Fig. 3(b)]. The solution is similar to those that eventually decay (i.e., those found at lower $\tilde{Ra}$). However, this solution has been time-integrated for over 3.3 diffusion times and shows no sign of decaying. There is clear scale separation between $E_K$ and $E_M$ at low $l$ with the magnetic and velocity fields peaking at $l \sim 2-3$ and $l \sim 7-10$, respectively [Fig. 4(a)]. The corresponding spectra in $m$ [Fig. 4(b)] each exhibit ~7-fold symmetry with an additional peak in the axisymmetric component. We note that these solutions are of a different nature to previously found hemispherical solutions, e.g., Ref. [41] which have been found at larger driving and are typically induced through favourable choices of boundary and/or symmetry conditions.

Interestingly, for this solution, the magnetic energy is an order of magnitude greater than the kinetic energy [Fig. 5(a)] despite the modest (“weak-field”) value of the Elsasser number ($\Lambda' \approx 0.017$). The viscous force is important [Fig. 4(c)] confirming this remains a viscously dominated solution on a weak-field branch, despite being nondipolar. Indeed, the hierarchy of forces is similar to viscous HD runs at low driving (cf. Fig. 1 from Ref. [19]) albeit with the Lorentz force also entering the primary balance. This leads to a “VMAC-balance” across a range of intermediate lengthscales and a complicated balance influenced by symmetry at the largest scales.
<figure><img src="image_7.png" /><figcaption>FIG. 4. Plots for a hemispherical solution ($\tilde{R}a' = 1.4$, $\epsilon_m = 4 \times 10^{-6}$). Snapshots at $t \approx 3.61$ of (a) the $l$ and (b) $m$ spectra of the kinetic and magnetic energy. For panel (b) the $m = 0$ contributions are plotted as crosses on the vertical axis at $m = 1$. (c) The $l$ spectra of solenoidal forces, averaged over the time interval $1.40 \le t \le 1.61$.</figcaption></figure>

The hemispherical nature of the solution is produced by a particular selection of preferred modes, which can be observed when the Lowes-Mauersberger spectrum is separated into its equatorially symmetric (where $l+m$ is even) and asymmetric (where $l+m$ is odd) parts (Fig. 6). The hemispherical solution [Fig. 6(a)] displays an alternating pattern of symmetric and asymmetric modes as $l$ increases; this results in a cancellation in one hemisphere. Conversely, a dipolar solution [Fig. 6(b)], is dominated by equatorially asymmetric modes only. The effect is perhaps most obvious in the lowest order modes of each type, $g_1^0$ and $g_2^0$ (i.e., the axial dipole and axial quadrupole components). These modes are plotted alongside energies in Fig. 5 as a function of time. For the hemispherical solution both modes settle to nonzero values whereas the quadrupolar component vanishes in the dipolar solution. These plots also demonstrate the steady nature of the solutions.

## 2. Nondipolar solutions

With a further increase in the Rayleigh number we remain on the viscously dominated branch but find alternative nondipolar solutions. Starting at $\tilde{R}a' \approx 1.5$, the solutions found are somewhat unusual in that they often are able to spend extended periods of time (often a magnetic diffusion time or longer) in a transitional solution before a reconfiguration to a final state.

Figure 7(a) shows the time series of the kinetic and magnetic energy of a typical simulation (for $\tilde{R}a' = 1.6$). This simulation was initialized from a strong-field dipolar solution in nearby parameter space. After an initial short transient period, by $t \approx 0.15$ the simulation first develops a large scale

<figure><img src="image_8.png" /><figcaption>FIG. 5. Time series plots of the kinetic and magnetic energies and Gauss coefficients for (a) a weak-field hemispherical dynamo ($\tilde{R}a' = 1.4$, $\epsilon_m = 4 \times 10^{-6}$) and (b) a weak-field dipolar dynamo ($\tilde{R}a' = 2.1$, $\epsilon_m = 8.3 \times 10^{-6}$). Note the separate vertical scales on each plot: the left hand scale corresponds to energies and the right hand scale to the axial dipole, $g_1^0$, and axial quadrupole, $g_2^0$.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

FIG. 6. Time averaged Lowes-Mauersberger spectrum, $g_l$, for a (a) weak-field hemispherical solution ($\tilde{R}a' = 1.4$, $\varepsilon_m = 4 \times 10^{-6}$) and (b) weak-field dipolar solution ($\tilde{R}a' = 2.1$, $\varepsilon_m = 8.3 \times 10^{-6}$). The red curve shows the complete spectra, $g_l$, whereas blue/purple points show only the equatorially asymmetric/symmetric parts of $g_l$, respectively.

magnetic field with peaks at $l = 4$ and $m = 7$ in the associated spectra. The clear nondipolar nature of the solution can be observed in the radial magnetic field at the outer boundary [Fig. 8(a)]. The preference of modes is confirmed in Figs. 9(a) and 9(b), which show the energy spectra at a typical time during this early phase of the simulation. The large scale magnetic field is generated despite the peak of the kinetic energy being in the $l = 10$ component. The solution is reasonably long-lived (for ~0.6 magnetic diffusion times).

At $t \approx 0.8$ the solution changes such that the magnetic energy strength reduces, the kinetic energy increases, and the end state of the solution develops equipartition between the energies. At this time, there is an associated reconfiguration of the field morphology. The radial magnetic field at the outer boundary [Fig. 8(c)] retains its earlier nondipolar nature albeit with reduced field strength and a more patchy distribution with regions of dominant flux. Figures 9(d) and 9(e) show the energy spectra in the final state of the simulation. There remains a clear scale separation in the solution between the flow (peaking at $l = 9$) and magnetic field (now peaking at $l = 2$).

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

FIG. 7. Time series plots of the kinetic and magnetic energy for a typical weak-field, nondipolar solution ($\tilde{R}a' = 1.6$, $\varepsilon_m = 4 \times 10^{-6}$). (a) simulation initialized from a strong-field solution in nearby parameter space; (b) simulation initialized from the end point of the simulation displayed in panel (a) with the field amplitudes reduced to 1% of their original values.
<figure><img src="image_11.png" /><figcaption>FIG. 8. Plots for a weak-field nondipolar solution ($\tilde{R}a' = 1.6$, $\epsilon_m = 4 \times 10^{-6}$). Snapshots of (a) $B_r$ on a spherical surface at the outer boundary and (b) isocontours of $\omega_z$ colored with $u_z\omega_z$, the helicity associated with the columnar flow (blue negative and green positive), during the early phase of the solution (at $t \simeq 0.45$); (c), (d) as panels (a), (b) but during the late phase of the solution (at $t \simeq 2.05$).</figcaption></figure>

The solenoidal force hierarchy shown in Figs. 9(c) (for the early part of the simulation) and 9(f) (for the end state of the simulation) provide additional information on the state the solution. In both figures, the Coriolis force peaks at the same scale as the flow. The role of the Lorentz force is the main difference between the plots at the two different times. At early times the Lorentz force enters the primary balance. While similar to strong-field dipolar solutions [19] that exhibit MAC-balance at large scales, this state additionally contains the viscous force in the primary balance, thus forming a "VMAC-balance" at large scale. In particular, this demonstrates that a balance involving more than three terms in the vorticity equation is viable for a working dynamo. This highlights the fact that our solutions are found in an under-explored region of the parameter space which uses a modest value of $\mathcal{E}$, but unusually small $\mathcal{E}_m$. This amounts to approaching Earth-like values through a different limit to most previous works. The force hierarchy (and field strength and spectra) are similar to that of the solutions found in Sec. III B 1. However, in contrast to that hemispherical solution, at later times the solution here has evolved to a thermal-wind (i.e., AC-) balance at large scales with the Lorentz force and viscous force weakening. This later hierarchy is reminiscent of weak-field dipolar solutions (e.g., Figs. 3(a)–3(c) of Ref. [19]). However, strong symmetry found in the buoyancy force at larger $\mathcal{E}_m$ (see the Appendix) is lost perhaps indicating its pivotal role for generating dipolar solutions on the weak-field branch. The small value of $\mathcal{E}_m$ means the inertial effects are too small to break the symmetry as might be the case in other regimes. Here the large Lorentz force is able to couple the symmetries, breaking the dominant symmetry in the flow.

The question remains as to whether the long-lived transitional phase seen in these solutions is a result of the choice of initial conditions or intrinsic to the solution. To confirm the preference for a
<figure><img src="image_12.png" /><figcaption>FIG. 9. Spectra plots for a weak-field nondipolar solution ($\tilde{R}a' = 1.6$, $\epsilon_m = 4 \times 10^{-6}$). Plots of (a) the $l$ and (b) $m$ spectra of the kinetic and magnetic energy (snapshots at $t \simeq 0.45$) and (c) the $l$ spectra of solenoidal forces (averaged over $0.45 \le t \le 0.46$); all during the early phase of the solution. (d)–(f) As (a)–(c) but during the late phase of the solution with $t \simeq 2.05$ [for panels (d), (e)] and averaged over $1.93 \le t \le 2.10$ [for panel (f)]. For panels (b), (e) the $m = 0$ contributions are plotted as crosses on the vertical axis at $m = 1$.</figcaption></figure>

chaotic transitional phase, a weak-field nondipolar case (specifically $\tilde{R}a' = 1.6$) was rerun from an initial state of weak seeds for both the velocity and magnetic fields. The fields of the final state of the solution were reduced to 1% of their amplitude and the simulation was rerun. Figure 7(b) shows the time series of the kinetic and magnetic energy over the full duration of this rerun simulation.

From the weak seed fields the magnetic energy grows and the field develops a morphology with peaks at $l = 2$ and $m = 7$. This state is unstable to a second solution (beyond $t \approx 2$) with larger magnetic field strength and a dominant $l = 4$ and $m = 7$ morphology, the same state found early in the original simulation. As found before, this solution itself is unstable (beyond $t \approx 2.5$) and transitions to a solution with spectra peaks at $l = 2$ and $m = 5$. Ultimately, this means the two runs (the original run and the run initialized with weak seed fields) achieve the same final state and the intermediate regime is a requirement.

Behaviour such as this, characterized by a dynamo mode which can be sustained for a significant amount of time, but which eventually destabilises to another mode with no fluctuations is not new. However, such behavior is usually met in the case of “self-killing” dynamos: dynamo solutions which can survive for very long transients, but eventually decay to the purely hydrodynamic solution [9,42].

### C. Drift rates

The flow patterns in weak-field solutions are dominated by columnar structures reminiscent of those found close to the onset of (nonmagnetic) convection. For $\mathcal{E} = 10^{-4}$, the preferred mode is $m = 7$ [cf. Fig. 3(b)]. These solutions could therefore be argued to be (weakly) magnetically adjusted hydrodynamic solutions. This raises the natural question of the effect of the magnetic field
<figure><img src="image_13.png" /><figcaption>FIG. 10. (a) Drift rate of convection columns, $\omega$, and (b) the wave-flow parameter, $W$, as a function of supercriticality, $\tilde{Ra}'$. Each point in each plot represents a single simulation and data has been time-averaged. The following symbols are used for the given regime: circles for weak-field dipolar; triangles for weak-field hemispherical; squares for weak-field nondipolar; crosses for failed dynamos; plus symbols for hydrodynamic (HD) convection. The solitary red circle represents the predicted value of the critical frequency, $\omega_c$, from linear theory [29]. Strong-field solutions have been omitted from this figure.</figcaption></figure>

on the drift rate and preferred mode of the convection columns and whether dynamo action is driven primarily by waves or flows. Figure 10(a) shows that the drift rate of convection columns, $\omega$, is adjusted in the presence of dynamo action near to the onset of convection. Branches for the dynamo solutions (colored points) peel off the main curve for hydrodynamic convection (black plus symbols) and maintain faster drift rates than would otherwise be the case in nonmagnetic runs. Indeed, the presence of the magnetic field actually increases the drift rate as $\tilde{Ra}'$ is increased, whereas nonlinear (nonmagnetic) convection suppresses the drift rate. This is in contrast to the linear theory for the onset of magnetoconvection [22,30], which actually predicts slower drift rates.

For $\mathcal{E}_m = 4 \times 10^{-6}$, failed dynamos (colored crosses) recover the values of $\omega$ for hydrodynamic convection, as might be expected. Nevertheless, the situation at $\mathcal{E}_m = 8.3 \times 10^{-6}$ is more intriguing since failed dynamo solutions maintain the $m = 4$ symmetry of the weak-field dipolar solution into the hydrodynamic state. This remains the case over many magnetic diffusion times and is thus unlikely to be a transient feature. Consequently, the drift rates of columns in the failed dynamo runs do not match with their counterparts of the purely hydrodynamic runs. Naturally, this outcome is a result of the choice of an initial condition with 4-fold symmetry but demonstrates bistability of hydrodynamic solutions. Interestingly, this effect is not observed at $\mathcal{E}_m = 4 \times 10^{-6}$ because, in that case, the preferred mode of the weak-field solutions matches that of the nonmagnetic convection (i.e., $m = 7$).

Figure 10(b) demonstrates a clear transition in the regime of hydrodynamic convection (black plus symbols) between solutions dominated by waves ($W < 1$) and flows ($W > 1$) as driving increases. Dynamo solutions (colored points), at both $\mathcal{E}_m = 8.3 \times 10^{-6}$ and $\mathcal{E}_m = 4 \times 10^{-6}$, delay this transition. In particular, solutions at $\mathcal{E}_m = 4 \times 10^{-6}$ have $W < 1$ suggesting the scale-separated hemispherical and nondipolar solutions found in this region of parameter space (see Sec. III B) are driven by waves. Conversely, at $\mathcal{E}_m = 8.3 \times 10^{-6}$, the dynamo onset is at sufficiently high driving to allow flows (since $W > 1$) to drive dynamos and produce dipolar solutions without scale separation.
### D. Connection with mean-field theory

The possibility of driving a kinematic mean-field dynamo utilizing the convective flow near onset in plane layer geometry has been considered [43]. A nonlinear solution for this problem was later obtained and highlighted the possibility of complex dynamics on such dynamos [20]. In a further extension to the approach, [44] considered the flow near onset in a model that partly replicated the geometry of Earth's core. The work relied on a perturbative approach, the so-called annulus model, as previously introduced in Ref. [37]. The model in effect applies only to the weak-field branch of dynamo action (since the Lorentz force is treated as a perturbation). The flow of Ref. [44] has properties consistent with our solutions. In particular, the helicity plots presented in Figs. 3(b) and 8(b), 8(d) highlight that most of the vertical flow here accounts for the sloping boundaries at the top and bottom of the convection columns. The average helicity thus vanishes at leading order in his flow. Despite the presence of some Ekman pumping, this is also the case in our realized solutions.

As highlighted by Busse [44], the stark contrast between the horizontal lengthscale and the vertical length-scale of the flow implies that a magnetic Reynolds number based on the small scale only is not relevant. In the presence of scale separation, the product of the small-scale and large-scale magnetic Reynolds numbers provide the relevant criterion for a mean-field dynamo. This is relevant to the simulations portrayed in Figs. 4(a) and 9(a). In both cases, the dominant lengthscale of the flow, $\ell$, corresponds to a spherical harmonic degree $l = 7$ and the magnetic Reynolds number corresponding to the flow at this scale is, respectively, 17.2 and 26.7. While these values are very low, they could, in theory, maintain a dynamo at this scale. However, a mean-field dynamo in the above sense will involve a lengthscale $\sqrt{r_o \ell}$ [1,45], which will be about three times larger, highlighting the importance of scale separation.

One of the puzzling observations in the results we report is the self-killing nature of the dynamo solution obtained at early times for $\tilde{Ra}' = 1.6$, $\mathcal{E}_m = 4 \times 10^{-6}$. Busse [44] finds a scale-separated dynamo with a toroidal field of the same order of magnitude as the poloidal field. Interestingly, the initial solution obtained with $\tilde{Ra}' = 1.6$ [see Figs. 8(a) and 8(b)], as well as the solution obtained at weaker forcing ($\tilde{Ra}' = 1.4$, see Fig. 3), are within a factor of two of equipartition between the toroidal, $E_M^{\text{Tor}}$, and poloidal, $E_M^{\text{Pol}}$, parts of the magnetic energy. Specifically, $E_M^{\text{Tor}}/E_M^{\text{Pol}} = 1.89$ and $E_M^{\text{Tor}}/E_M^{\text{Pol}} = 1.67$ for $\tilde{Ra}' = 1.6$ and $\tilde{Ra}' = 1.4$, respectively. The final, and stable, solution obtained at later time (corresponding to the final state in Fig. 7), is however characterized by a much stronger toroidal field, with $E_M^{\text{Tor}}/E_M^{\text{Pol}} = 4.95$. This mode therefore corresponds to a different sort of weak-field dynamo, which could involve a mean-field $\Omega$-effect.

### E. Strong-field solutions at $\mathcal{E}_m = 4 \times 10^{-6}$

A further increase in the driving to $\tilde{Ra}' = 1.8$ leads to the emergence of a strong-field dipolar branch for $\mathcal{E}_m = 4 \times 10^{-6}$, as seen originally in Fig. 1. The characteristics of these solutions are broadly similar to strong-field branches identified previously (e.g., Sec. III A and Refs. [10,19]). In particular, these solutions do not arise via a continuous path from the solutions at lower driving; rather, there is a sharp change in the solution for a infinitesimally small change in the value of the Rayleigh number. The sharp transition is identified by a jump in the value of $\Lambda'$ and key differences in the energy and force hierarchy spectra of the strong-field solutions (compared with weak-field solutions).

The field on the outer boundary and energy spectra for a typical solution on the strong-field branch are shown in Fig. 11. Both the field strength and morphology and the energy spectra [Figs. 12(a) and 12(b)] are typical for a strong-field solution. A plot of the solenoidal forces [Fig. 12(c)] also demonstrates a MAC balance at large lengthscales transitioning to an MC-balance for smaller scales with viscous forces only coming into the leading balance at the smallest scales. These are clear signatures of a strong-field dipolar branch solution [19] confirming its existence at $\mathcal{E}_m = 4 \times 10^{-6}$ despite the different nature of the solutions on the associated weak-field branch presented earlier. Interestingly, once the solution has transitioned to a strong-field case, the ratio
<figure><img src="image_14.png" /><figcaption>FIG. 11. Plots for a strong-field dipolar solution ($\tilde{R}a' = 1.8$, $\mathcal{E}_m = 4 \times 10^{-6}$). Snapshots at $t \simeq 1.54$ of (a) $B_r$ on a spherical surface at the outer boundary and (b) isocontours of $\omega_z$ colored with $u_z \omega_z$, the helicity associated with the columnar flow.</figcaption></figure>

of the toroidal to poloidal magnetic energy is strongly reduced compared to the final state of the weak-field case at $\tilde{R}a = 1.6$; specifically, $E_M^{\text{Tor}}/E_M^{\text{Pol}} = 1.73$ for $\tilde{R}a' = 1.8$.

## IV. DISCUSSION

We have presented results of dynamo simulations in a setup relevant to Earth's core although, as with all DNS models, the input parameters used are restricted (due to numerical constraints) to values far from the true physical values. Our aim was to explore the weak-field branch within dynamo simulations in greater detail than previous studies [10,26]. To achieve this a sufficiently large value of the magnetic Prandtl number is required because of computational limitations on the value of $\mathcal{E}$. While this may seem counterintuitive (since the value of $P_m$ in Earth's core is very small), the value of $P_m$ at which the existence of distinct weak- and strong-field solutions can be found is expected to reduce through a distinguished limit as $\mathcal{E}$ also becomes small [10]. It is therefore appropriate to consider dynamo solutions at “large” $P_m$ to understand their behavior at the moderate values of $\mathcal{E}$ available for numerical study [46]. It can be insightful to instead present results in terms of the magnetic Ekman number (rather than $P_m$), as we choose to do here. This is because $\mathcal{E}_m$ needs to be small both for giving rise to weak- and strong-field branches in DNS and for realistic core values. To adequately broaden the survey of weak-field solutions in parameter space,

<figure><img src="image_15.png" /><figcaption>FIG. 12. Plots for a strong-field dipolar solution ($\tilde{R}a' = 1.8$, $\mathcal{E}_m = 4 \times 10^{-6}$). Snapshots at $t \simeq 1.54$ of (a) the $l$- and (b) the $m$ spectra of the kinetic and magnetic energy. For (b) the $m = 0$ contributions are plotted as crosses on the vertical axis at $m = 1$. (c) The $l$ spectra of solenoidal forces, averaged over $1.16 \le t \le 1.23$.</figcaption></figure>
we have focused on single (moderate) values of the Ekman number and Prandtl number. Future work could consider the scaling with $\mathcal{E}$, of solutions presented here. Our work has three main discussion points: the variety of solutions on the weak-field branch; the introduction of solutions with scale separation on the weak-field branch; the persistence of the strong-field branch and identification of weak- versus strong-field solutions. We address each in turn.

Fundamentally, the results presented demonstrate that the viscously dominated weak-field branch is richer in regimes than previously identified. At small enough magnetic Ekman number ($\mathcal{E}_m$) not only does branching materialise with weak-field dipolar solutions (WD, see Ref. [26]), but the weak-field branch can admit nondipolar solutions. Hemispherical dynamos and solutions that are dominated by a combination of magnetic quadrupolar and octupolar components are found in the latter category of weak-field solutions. A decreased $\mathcal{E}_m$ leads to dynamo onset moving ever closer (in $\mathcal{R}a$ space) to the onset of nonmagnetic convection (i.e., $\mathcal{R}a' = 1$). The result is dynamo regimes that share characteristics with solutions near the onset of thermal convection. In particular, we have shown that the dynamos at $\mathcal{E}_m = 4 \times 10^{-6}$ (equivalently $P_m = 25$ for $\mathcal{E} = 10^{-4}$), where the weak-field branch becomes nondipolar, are driven primarily by wave motion rather than flows. Such solutions share the same mode of convection as the hydrodynamic state, unlike the dipolar solutions at $\mathcal{E}_m = 8.3 \times 10^{-6}$ (equivalently $P_m = 12$ for $\mathcal{E} = 10^{-4}$). The dominance of larger scale flows, symmetry in the buoyancy force, and larger modes of convection are likely the cause of the development of dipolar solutions. The hemispherical dynamos discussed exist firmly on the weak-field branch of dynamo action close to the onset of (nondecaying) dynamo solutions, in a model using fixed temperature conditions and differential heating. This is in contrast to previous studies which typically employ fixed temperature flux conditions (combined with internal heating) to promote hemispherical solutions. Such conditions are considered potentially favourable since hemispherical dynamos require some form of north-south symmetry breaking which may be aided by stronger variation in temperature gradients between hemispheres. Indeed, the studies most focused on generating hemispherical dynamos often apply north-south asymmetric (i.e., heterogeneous) boundary conditions or forcing directly to enhance the effect further, e.g., Ref. [47]. Furthermore, in all such scenarios, the solutions are found at larger driving than the weak-field solutions presented in our work; hence our solutions are a new form of hemispherical regime. For example, a recent study [48], uses fixed flux boundary conditions and internal heating to find hemispherical wave solutions located in $\mathcal{R}a$ space between the dipolar and multipolar regimes.

Solutions on the weak-field branch at $\mathcal{E}_m = 4 \times 10^{-6}$ are also distinguished from their dipolar counterparts at $\mathcal{E}_m = 8.3 \times 10^{-6}$ by a significant separation of scales between the velocity and magnetic fields. This scale separation offers a connection with mean-field theory. Interestingly, the scale separation is possible despite the small value of $\mathcal{E}_m$, implying here $P_m > 1$, which one might naively expect to create small scale magnetic field. It also confirms that dynamos with small values of $\mathcal{E}_m$ are of relevance to natural dynamos. The separation of scales exhibited here may not be as distinct as mean-field theory envisages, but it nevertheless demonstrates that a large scale magnetic field can be driven in weak-field solutions by velocity perturbations on a clearly different scale. As the Ekman number is reduced below the $10^{-4}$ value we have studied here, the resultant fundamental mode of convection will grow and it can therefore be speculated that the separation of scales will widen. Simulations at $\mathcal{E} < 10^{-4}$, combined with small $\mathcal{E}_m$, are computationally challenging and beyond the scope of this initial study. Previous studies have shown that scale separation remains possible at lower Ekman number [49] although the effect was not prominent at larger values of the magnetic Prandtl number ($P_m > 0.4$). However, it is unclear which regime/dynamo branch is being studied in that work especially given an alternative setup of internal heating and a fixed flux thermal boundary condition on the outer boundary.

Given the existence of weak-field nondipolar solutions at $\mathcal{E}_m = 4 \times 10^{-6}$, there was no guarantee of the emergence (at this $\mathcal{E}_m$) of the dipolar strong-field regime and transition between the weak- and strong-field branches known from $\mathcal{E}_m = 8.3 \times 10^{-6}$. Our work has shown, however, that the strong-field regime persists at $\mathcal{E}_m = 4 \times 10^{-6}$ and shares the same fundamental characteristics of
its counterpart at larger $\mathcal{E}_m$. The similarity of the kinetic and magnetic energy spectra in our scale-separated solutions and those found previously (see, e.g., Fig. 3 of Ref. [49]) points to a common regime of all these solutions (i.e., weak-field). However, the lower $\mathcal{E}$ simulations are more dynamically active (i.e., differing from our steady solutions) and have $E_K < E_M$, often cited as an indicator of “strong-field behavior.” Nevertheless, the danger in using this criterion to label a solution as “strong-field” is twofold. First, as we have seen in this work, solutions that are indisputably on the weak-field branch can possess $E_K < E_M$ combined with a Lorentz force with significant magnitude; crucially, however, viscous force remains in the primary balance and the Elsasser number remains small. Evidently, the magnitude of the magnetic energy is not (alone) a reliable determining factor of weak- versus strong-field solutions. Second, the variety of weak-field branches now identified as $\mathcal{E}_m$ is varied indicates a richness of regimes that may extend to dynamically active solutions at lower $\mathcal{E}$ and/or $\mathcal{E}_m$. Ultimately, if strong-field solutions are desired (e.g., due to their relevance to natural dynamos), then it would be prudent to investigate the balance of (solenoidal) forces and/or check for the existence of distinct weak- and strong-field branches of dynamo action.

## ACKNOWLEDGMENTS

The authors thank the anonymous reviewers for their useful comments. This work was supported by the Science and Technology Facilities Council (Grant No. ST/Y00146X/1); the DiRAC Project No. ACTP245. Part of the numerical work for this publication was performed using the DiRAC Memory Intensive service (Cosma8/Cosma7/Cosma6 [*]) at Durham University, managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility [50]. The DiRAC service at Durham was funded by BEIS, UKRI, and STFC capital funding, Durham University, and STFC operations grants. DiRAC is part of the UKRI Digital Research Infrastructure.

## DATA AVAILABILITY

The data that support the findings of this article are openly available [33], embargo periods may apply.

## APPENDIX: SPECTRA OF WEAK- AND STRONG-FIELD BRANCHES AT $\mathcal{E}_m = 8.3 \times 10^{-6}$

We present here further discussion of the nature of the solutions on the weak- and strong-field branches at $\mathcal{E}_m = 8.3 \times 10^{-6}$ (mentioned in Sec. III A), to aid comparison with our key results found at $\mathcal{E}_m = 4 \times 10^{-6}$. As such the spectra plots displayed in Fig. 13 are results from the simulations presented in Fig. 2.

The saw-tooth pattern in the energy spectrum of the weak-field solution [Fig. 13(a)] demonstrates that particular modes are preferred, producing highly symmetric flow and field. The 4-fold symmetry in azimuth of the weak-field dipolar solution (originally seen in its radial field at the outer boundary) is further revealed in Fig. 13(b). The strong-field solution [Figs. 13(d) and 13(e)] has smoother and broader energy spectra in both $l$ and $m$ with no preferred symmetry. Neither the weak-field nor strong-field solutions at $\mathcal{E}_m = 8.3 \times 10^{-6}$ display clear scale separation in their spectra. This is in contrast with the results of Sec. III B.

The sharp distinction between the weak- and strong-field solutions is also underscored by the stark difference in their lengthscale dependent force hierarchy [Figs. 13(c) and 13(f)]. The viscous force remains important at all scales in the weak-field solution and, in particular, enters a primary balance with the Coriolis force for odd values of $l$. Due to the symmetric nature of the solution, the buoyancy force is only in primary balance (with the Coriolis force) for even $l$. The Lorentz force is mostly a secondary contributor throughout except for $l = 3$ where it (and inertia) come into balance with a VC balance that exists at all other odd $l$. The inertial force is also important at $l = 1$ and both the inertial and Lorentz terms enter the balance again at very small scales. Elsewhere they are secondary contributors throughout. Conversely, the Lorentz force enters the primary balance in the
<figure><img src="image_16.png" /></figure>

FIG. 13. *l* and *m* spectra of energy for: (a), (b) a typical weak-field branch solution ($\tilde{R}a' = 2.1$, $\epsilon_m = 8.3 \times 10^{-6}$ at $t \simeq 1.13$); (c), (e) a typical strong-field branch solution ($\tilde{R}a' = 5$, $\epsilon_m = 8.3 \times 10^{-6}$ at $t \simeq 0.23$). These correspond to the same parameters as in Fig. 2. For panels (b), (e) the *m* = 0 contributions are plotted as crosses on the vertical axis at *m* = 1. Time-averaged *l* spectra of solenoidal forces are presented in panels (c), (f).

strong-field solution with the viscous force weak at all but the smallest scales. The inertial force is a
secondary contributor, which is expected at low supercriticalities. Hence the weak-field solution can
be predominantly described as VAC and the strong-field solution as MAC (at least at large scales).

[1] E. Dormy, Rapidly rotating magnetohydrodynamics and the geodynamo, Annu. Rev. Fluid Mech. **57**, 335 (2025).

[2] C. Guervilly, P. Cardin, and N. Schaeffer, Turbulent convective length scale in planetary cores, *Nature* (London) **570**, 368 (2019).

[3] M. Steenbeck, F. Krause, and K.-H. Rädler, Berechnung der mittleren Lorentz-Feldstärke für ein elektrisch leitendes Medium in turbulenter, durch Coriolis-Kräfte beeinflußter Bewegung, Z. Naturforsch. **21**, 369 (1966) [English Translation: Roberts & Stix (1971), pp. 29–47].

[4] F. Krause and K.-H. Rädler, *Mean-Field Magnetohydrodynamics and Dynamo Theory* (Pergamon Press, Oxford, UK, 1980).

[5] H. K. Moffatt and E. Dormy, *Self-Exciting Fluid Dynamos*, Cambridge Texts in Applied Mathematics (Cambridge University Press, Cambridge, UK, 2019).

[6] K. Mizerski, *Non-Equilibrium Hydromagnetic Dynamos* (IOP Publishing, Bristol, UK, 2023).

[7] U. Christensen, P. Olson, and G. Glatzmaier, Numerical modelling of the geodynamo: A systematic parameter study, *Geophys. J. Int.* **138**, 393 (1999).

[8] C. Kutzner and U. Christensen, From stable dipolar towards reversing numerical dynamos, *Phys. Earth Planet. Int.* **131**, 29 (2002).
[9] V. Morin and E. Dormy, The dynamo bifurcation in rotating spherical shells, *Int. J. Mod. Phys. B* **23**, 5467 (2009).

[10] E. Dormy, Strong-field spherical dynamos, *J. Fluid Mech.* **789**, 500 (2016).

[11] N. Schaeffer, D. Jault, H.-C. Nataf, and A. Fournier, Turbulent geodynamo simulations: A leap towards Earth's core, *Geophys. J. Int.* **211**, 1 (2017).

[12] J. Aubert, State and evolution of the geodynamo from numerical models reaching the physical conditions of Earth's core, *Geophys. J. Int.* **235**, 468 (2023).

[13] C. Jones, A dynamo model of jupiter's magnetic field, *Icarus* **241**, 148 (2014).

[14] M. Schrinner, K.-H. Rädler, D. Schmitt, M. Rheinhardt, and U. R. Christensen, Mean-field concept and direct numerical simulations of rotating magnetoconvection and the geodynamo, *Geophys. Astrophys. Fluid Dyn.* **101**, 81 (2007).

[15] M. Schrinner, L. Petitdemange, and E. Dormy, Oscillatory dynamos and their induction mechanisms, *Astron. Astrophys.* **530**, A140 (2011).

[16] K. Soderlund, A. Sheyko, E. King, and J. Aurnou, The competition between Lorentz and Coriolis forces in planetary dynamos, *Prog. Earth Planet. Sci.* **2**, 24 (2015).

[17] T. Schwaiger, T. Gastine, and J. Aubert, Force balance in numerical geodynamo simulations: A systematic study, *Geophys. J. Int.* **219**, S101 (2019).

[18] T. Schwaiger, T. Gastine, and J. Aubert, Relating force balances and flow length scales in geodynamo simulations, *Geophys. J. Int.* **224**, 1890 (2020).

[19] R. J. Teed and E. Dormy, Solenoidal force balances in numerical dynamos, *J. Fluid Mech.* **964**, A26 (2023).

[20] A. Soward, A convection-driven dynamo I. The weak field case, *Philos. Trans. R. Soc. London A* **275**, 611 (1974).

[21] A. Soward, Convection driven dynamos, *Phys. Earth Planet. Inter.* **20**, 134 (1979).

[22] D. R. Fearn, Thermal and magnetic instabilities in a rapidly rotating fluid sphere, *Geophys. Astrophys. Fluid Dyn.* **14**, 103 (1979).

[23] S. Horn and J. M. Aurnou, The Elbert range of magnetostrohic convection. I. Linear theory, *Proc. R. Soc. A* **478**, 20220313 (2022).

[24] S. Horn and J. M. Aurnou, The Elbert range of magnetostrophic convection. II. Comparing linear theory to nonlinear low-$Rm$ simulations, *Proc. R. Soc. A* **481**, 20240016 (2025).

[25] P. Roberts, Pure dynamo theory and geomagnetic dynamo theory, in *Proceedings of the 1st International Workshop on Dynamo Theory and the Generation of the Earth's Magnetic Field*, edited by I. Cupal (Czech Academy of Sciences, Prague, 1979), pp. 7–12.

[26] E. Dormy, L. Oruba, and L. Petitdemange, Three branches of dynamo action, *Fluid Dyn. Res.* **50**, 011415 (2018).

[27] R. D. Simitev and F. H. Busse, Bistability and hysteresis of dipolar dynamos generated by turbulent convection in rotating spherical shells, *Europhys. Lett.* **85**, 19001 (2009).

[28] L. Petitdemange, Systematic parameter study of dynamo bifurcations in geodynamo simulations, *Phys. Earth Planet. Inter.* **277**, 113 (2018).

[29] E. Dormy, A. Soward, C. Jones, D. Jault, and P. Cardin, The onset of thermal convection in rotating spherical shells, *J. Fluid Mech.* **501**, 43 (2004).

[30] C. A. Jones, A. I. Mussa, and S. J. Worland, Magnetoconvection in a rapidly rotating sphere: The weak-field case, *Proc. R. Soc. London A* **459**, 773 (2003).

[31] A. Willis, B. Sreenivasan, and D. Gubbins, Thermal core–mantle interaction: Exploring regimes for “locked” dynamo action, *Phys. Earth Planet. Inter.* **165**, 83 (2007).

[32] “Leeds spherical dynamo” code, https://github.com/Leeds-Spherical-Dynamo.

[33] R. J. Teed, Data for: Toward scale separated weak-field spherical dynamos (2025), https://doi.org/10.5525/gla.researchdata.1935.

[34] P. Mauersberger, Das mittel der energiedichte des geomagnetischen hauptfeldes an der Erdoberfläche und seine säkulare anderung, *Gerlands Beitr. Geophys.* **65**, 207 (1956).

[35] F. Lowes, Spatial power spectrum of the main geomagnetic field, and extrapolation to the core, *Geophys. J. Int.* **36**, 717 (1974).
[36] P. Roberts, On the thermal instability of a rotating-fluid sphere containing heat sources, *Philos. Trans. R. Soc. London A* **263**, 93 (1968).

[37] F. H. Busse, Thermal instabilities in rapidly rotating systems, *J. Fluid Mech.* **44**, 441 (1970).

[38] E. Dormy, Modélisation numérique de la dynamo terrestre, Ph.D. thesis, Institut de Physique du Globe de Paris (1997).

[39] P. Roberts and A. Soward, Dynamo theory, *Annu. Rev. Fluid Mech.* **24**, 459 (1992).

[40] A. Guseva, L. Petitdemange, and S. M. Tobias, Run-away transition to turbulent strong-field dynamo, *JGR Planets* **130**, e2024JE008496 (2025).

[41] W. Dietrich and J. Wicht, A hemispherical dynamo model: Implications for the martian crustal magnetization, *Phys. Earth Planet. Inter.* **217**, 10 (2013).

[42] H. Fuchs, K.-H. Rädler, and M. Rheinhard, On self-killing and self-creating dynamos, *Astron. Nachr.* **320**, 129 (1999).

[43] F. H. Busse, Generation of magnetic fields by convection, *J. Fluid Mech.* **57**, 529 (1973).

[44] F. H. Busse, A model of the geodynamo, *Geophys. J. R. Astron. Soc.* **42**, 437 (1975).

[45] C. Guervilly and E. Dormy, The cross-over from viscous to inertial lengthscales in rapidly-rotating convection, *Geophys. Res. Lett.* **52**, e2024GL111593 (2025).

[46] Indeed, if $P_m$ is chosen to be too small (for a given $\mathcal{E}$), then dynamo action ceases altogether!

[47] S. Stanley, L. Elkins-Tanton, M. T. Zuber, and E. M. Parmentier, Mars' paleomagnetic field as the result of a single-hemisphere dynamo, *Science* **321**, 1822 (2008).

[48] Y. Wu, W. Dietrich, and X. Tao, Parameter regimes of hemispherical dynamo waves in a spherical shell from 3D MHD simulations, *JGR Planets* **129**, e2023JE007976 (2024).

[49] A. Sheyko, C. Finlay, J. Favre, and A. Jackson, Scale separated low viscosity dynamos and dissipation within the Earth's core, *Sci. Rep.* **8**, 12566 (2018).

[50] www.dirac.ac.uk