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

HAL
open science

# Scaling of Strong-Field Spherical Dynamos

Robert J Teed, Emmanuel Dormy

► To cite this version:

Robert J Teed, Emmanuel Dormy. Scaling of Strong-Field Spherical Dynamos. Geophysical Research Letters, 2025, 52 (20), pp.e2025GL118078. 10.1029/2025GL118078. hal-05407953

HAL Id: hal-05407953

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

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
# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2025GL118078

**Key Points:**

* The Geodynamo is strongly influenced by rotation and magnetic forces, conditions which numerical models struggle to achieve

* A new output parameter, the Curwsasser number, is introduced to help identify appropriate numerical solutions to study the Geodynamo

* Strong-field solutions in rapid rotation provide a region of numerical parameter space to target future Geodynamo simulations

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

R. J. Teed and E. Dormy,
Robert.Teed@glasgow.ac.uk;
Emmanuel.Dormy@ens.fr

**Citation:**

Teed, R. J., & Dormy, E. (2025). Scaling of strong-field spherical dynamos. *Geophysical Research Letters*, *52*, e2025GL118078. https://doi.org/10.1029/2025GL118078

Received 17 JUL 2025
Accepted 29 SEP 2025

## Scaling of Strong-Field Spherical Dynamos

Robert J. Teed¹ and Emmanuel Dormy²

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

¹School of Mathematics & Statistics, University of Glasgow, Glasgow, UK, ²Département de Mathématiques et Applications, UMR-8553, École Normale Supérieure, CNRS, PSL University, Paris, France

**Abstract** Numerical experiments of dynamo action designed to understand the generation of Earth's magnetic field produce different regime branches identified within bifurcation diagrams. Notable are distinct branches where the resultant magnetic field is either weak or strong. Weak-field solutions are identified by the prominent role of viscosity (and/or inertia) on the motion, whereas the magnetic field has a leading-order effect on the flow in strong-field solutions. We demonstrate the persistence of the strong-field branch, preserving the expected force balance of Earth's core, and provide scaling laws governing its onset as parameters move toward values appropriate for the Geodynamo. We introduce a new output parameter, based on dynamically important parts of rotational and magnetic forces, that captures expected $O(1)$ values of strong-field solutions throughout input parameter space. This new measure of the field strength and our bounds on scaling laws can guide future studies in locating strong-field dynamos in parameter space.

**Plain Language Summary** Earth's magnetic field is produced by turbulent motions of an electrically conducting fluid within the core of the planet. Since Earth's deep interior is not directly observable, numerical simulations are routinely used to model the fluid flow and magnetic field. Due to present-day computational limitations, such models cannot replicate the conditions of Earth's core but can guide us on the underlying processes. Key to obtaining informative simulations in the future is to understand how solutions with relevant properties change within the currently accessible numerical parameter ranges. In this work we study numerical solutions that display key attributes of Earth's magnetic field. In particular, we focus on models that exhibit both the magnetic field's north-south dipole dominance and its strong influence on the flow. We find that these most suitable solutions persist as the conditions of simulations are moved toward those of Earth's core. However, this is only the case when simulation inputs are carefully selected and so we provide new “tools” to identify and maintain the solutions. Our work will guide future numerical experiments seeking to study the generation of Earth's magnetic field as computing power advances.

## 1. Introduction

The inability to sample conditions within Earth's core directly means direct numerical simulations (DNS) can provide an essential role in progressing our understanding of the Geodynamo mechanism (as well as dynamos of other planets). Varying the key control parameters of the system ($\mathcal{E}$, $\mathcal{E}_m$, $q$, and $\tilde{R}a$ - each is introduced formally later) maps out part of the solution space. Yet, due to a lack of computational resolution, the regimes identified in DNS do not operate under conditions of the core. In particular, viscosity must be set many orders of magnitude too large. Key to making progress is identifying regions of solution space most appropriate to Earth's core at the moderate input parameters available computationally.

Early pioneering work using simulations (Christensen et al., 1999; Kutzner & Christensen, 2002) discovered two distinct regimes of Geodynamo action in the accessible parameter space. Dipolar solutions were found at the onset of dynamo action (above the onset of convection), transitioning to a multipolar regime at higher convective driving. Later work (Morin & Dormy, 2009) focused on bifurcations in $\tilde{R}a$-space and showed the existence of subcritical bifurcations and isola. A systematic route to studying subcritical dynamos has recently been introduced by Mannix et al. (2022). More recent work identified different branches of dynamo action solely within the dipolar regime. It was shown that a bifurcation diagram with two dipolar branches is possible in which the dynamo onsets on a weak-field dipolar branch before transitioning to a strong-field dipolar branch rather than a multipolar branch (Dormy, 2016). For a range of $\tilde{R}a$, the weak- and strong-field branches are bistable resulting in a double fold and a cusp singularity in the parameter space (see Dormy, 2025, for a recent review of these ideas). Dormy et al. (2018) then considered lower $\mathcal{E}$; however, with no tracking of bifurcations, the reported strong-field

<footnote>© 2025. The Author(s).
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>
dynamos were based on unconfirmed scaling laws, not on observation of separate branches and bistability. Menu et al. (2020) further showed that the strong-field solution exists at a level of forcing for which the weak-field destabilizes to a multipolar state. Recently, Guseva et al. (2025) were able to provide a mathematical framework for studying the destabilization of the weak-field solution to the strong-field regime.

The importance of the strong-field branch lies in its solutions providing an anticipated balance of forces for Earth's core, in contrast to the other branches discussed. MAC-balance, between magnetic (i.e., Lorentz), Archimedean (i.e., buoyancy), and Coriolis forces, is expected to prevail in Earth's core (see also discussions in Buffett, 2000; Horn & Aurnou, 2022, 2025; Mason et al., 2022; Soderlund et al., 2012, 2015). The strong-field regime is thus a part of solution space requiring further investigation to check its existence and scaling as parameters are gradually moved toward those of Earth's core.

In this work we demonstrate the persistence of the strong-field branch as $\mathcal{E}$ and $\mathcal{E}_m$ are lowered toward values appropriate to the Geodynamo. We investigate scaling laws, as a function of $\mathcal{E}$, for the critical value of $\mathcal{E}_m$ below which distinct weak- and strong-field branches can be found. Furthermore, we introduce an alternative definition of the Elsasser number based on the ratio of curls of the Lorentz and Coriolis forces. In contrast to previous definitions, this new output parameter is able to capture the expected $O(1)$ values for strong-field solutions and maintain this across a wide range of values of $\mathcal{E}$.

## 2. Methods

### 2.1. Governing Equations

We consider a spherical shell (bounded between $r = r_i$ and $r = r_o$ with aspect ratio, $\chi = r_i/r_o$) filled with a conducting Boussinesq fluid of density, $\rho_0$, viscosity, $\nu$, thermal diffusivity, $\kappa$, magnetic diffusivity, $\eta$, and rotating with rate $\Omega$ about the vertical. A temperature difference of $\Delta T$ is applied between the inner and outer boundary and gravity takes the form $g = -g_0\mathbf{r}$ where $\mathbf{r} = (1-\chi)r \mathbf{e}_r$. The variables are nondimensionalized using lengthscale, $d = r_o-r_i$, timescale, $d^2/\eta$, temperature scale, $\Delta T$, and magnetic scale, $\sqrt{\mu_0\rho_0\eta\Omega}$. The evolution of the velocity of the fluid, $\mathbf{u}$, the temperature $T$, and the magnetic induction, $\mathbf{B}$, are governed by the following set of non-dimensionalized equations

$$ \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 $p$ is the pressure. The nondimensional parameters in Equations 1a–1e are the Ekman number, $\mathcal{E}$, magnetic Ekman number, $\mathcal{E}_m$, modified Rayleigh number, $\tilde{R}a$, and Roberts number, $q$, defined as

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

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

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

$$ q = \frac{\kappa}{\eta}, \quad (2d) $$
where $\alpha$ is the coefficient of thermal expansion. The Rayleigh number used in this work is related to the “rotational” Rayleigh number, $Ra = \alpha g_0 \Delta T d / \Omega \kappa$, and the classical (non-rotational) Rayleigh number, $R = \alpha g_0 \Delta T d^3 / \nu \kappa$, through $\tilde{Ra} = qRa = qER$. For convenience, the Rayleigh number is often measured by its supercriticality, $\tilde{Ra}' = \tilde{Ra}/qRa_c = Ra/Ra_c$, where $Ra_c$ is the value of $Ra$ for the onset of non-magnetic convection (for a given $E$). The governing Equations 1a–1e, subject to impenetrable, rigid, electrically-insulating, isothermal boundary conditions with $\chi = 0.35$, are solved using the Leeds spherical dynamo code (Willis et al., 2007). Further details of our methods are given in Text S1 in Supporting Information S1.

## 2.2. Elsasser Number and Curlsasser Number

The Elsasser number, $\Lambda$, is a measure of the ratio of the Lorentz force to the Coriolis force. It is often introduced to determine whether a particular dynamo solution has a “strong” magnetic field with $\Lambda \sim O(1)$, suggestive of a dominant balance between Lorentz and Coriolis forces as part of a wider MAC-balance.

In its broadest definition the Elsasser number can be written as

$$ \Lambda = \frac{\{\mathbf{F}_L\}}{\{\mathbf{F}_C\}} = \frac{\{(\mu_0 \rho_0)^{-1} (\nabla \times \mathbf{B}) \times \mathbf{B}\}}{\{2 \Omega \mathbf{e}_z \times \mathbf{u}\}} = \frac{B^2}{2 \mu_0 \rho_0 \Omega U \ell_B}, \quad (3) $$

where $F_L$ and $F_C$ are the Lorentz and Coriolis forces, respectively, and $\{\cdot\}$ is an “order-of-magnitude” operator indicating some process of approximating typical values and averaging. In Equation 3, $U$ and $B$ are typical (e.g., rms) values of the velocity and magnetic field, respectively and $\ell_B$ is a lengthscale appropriate for the magnetic field.

One can then make choices for measures of $U$, $B$, and $\ell_B$, leading to two common versions of the Elsasser number:

$$ \tilde{\Lambda} = \frac{B^2}{2\mu_0\rho_0\Omega\eta} \qquad \text{for } U \sim \frac{\eta}{d}, \text{ with } \ell_B \sim d, \quad (4) $$

$$ \Lambda' = \frac{B^2 d}{2\mu_0\rho_0\Omega\eta R_m \ell_\eta} = \tilde{\Lambda} \frac{d}{R_m \ell_\eta} \qquad \text{for } U \sim \frac{R_m}{d/\eta}, \text{ with } \ell_B \sim \ell_\eta. \quad (5) $$

Here, $R_m$ is the magnetic Reynolds number based on the rms value of the velocity and

$\ell_\eta^2 = \int_V \mathbf{B}^2 dV / \int_V (\nabla \times \mathbf{B})^2 dV$ is a measure of the magnetic dissipation lengthscale. The classical definition of the Elsasser number is given by Equation 4 whereas $\Lambda'$ is a “modified” (Dormy, 2016) or “dynamic” (Soderlund et al., 2012) Elsasser number.

The leading order balance of forces is however partly disguised by the large pressure gradient force. Recent work (Teed & Dormy, 2023) has shown that the leading order balance of forces can be readily recovered by filtering out dynamically unimportant gradient parts of all forces, forming solenoidal forces. If the purpose of the Elsasser number is to indicate balance between Lorentz and Coriolis forces, another logical option for quantifying the influence of the Lorentz force is to take the ratio of the curls of the original forces. We therefore introduce a new, previously unused measure, the “Curlsasser number,” given by

$$ \sigma = \frac{\{\nabla \times \mathbf{F}_L\}}{\{\nabla \times \mathbf{F}_C\}} = \frac{\{(\mu_0 \rho_0)^{-1} \nabla \times ((\nabla \times \mathbf{B}) \times \mathbf{B})\}}{\{2\Omega \nabla \times (\mathbf{e}_z \times \mathbf{u})\}} = \frac{B^2 \ell_z}{2\mu_0 \rho_0 \Omega U \ell_B \ell_C} = \Lambda \frac{\ell_z}{\ell_C}, \quad (6) $$

where $\ell_z$ is the vertical lengthscale of the flow and $\ell_C$ is some correlation lengthscale coupling the magnetic field and its curl. In practice the quantity on the RHS of Equation 6 will be difficult to calculate because of the uncertainty in measuring $\ell_C$, although, in a similar vein to $\ell_\eta$, one could approximate $\ell_z^2 = \int_V (\mathbf{e}_z \times \mathbf{u})^2 dV / \int_V (\nabla \times (\mathbf{e}_z \times \mathbf{u}))^2 dV$.
The uncertainties in the definitions given by Equations 3 and 6 can be avoided by calculating the relevant forces and/or curls of forces directly. In this work we therefore present $\Lambda$ based on the ratio of the directly calculated forces and $\sigma$ based on the ratio of the directly calculated curls of forces; that is we set

$$ \Lambda = \sqrt{\frac{\int_V ((\nabla \times \mathbf{B}) \times \mathbf{B})^2 dV}{\int_V (2\mathbf{e}_z \times \mathbf{u})^2 dV}}, \qquad (7a) $$

$$ \sigma = \sqrt{\frac{\int_V (\nabla \times ((\nabla \times \mathbf{B}) \times \mathbf{B}))^2 dV}{\int_V (2\nabla \times (\mathbf{e}_z \times \mathbf{u}))^2 dV}}. \qquad (7b) $$

The measures given by Equations 7a–7b provide a direct estimate of the relevant force ratio unlike those given by Equations 3–6, which rely on choices being made for the relevant lengthscales and flow/field magnitudes. Both $\Lambda$ and $\sigma$ are derived from measurable quantities and can be readily calculated from simulation data at no extra computational cost compared to $\tilde{\Lambda}$ and $\tilde{\Lambda}'$. Indeed, Soderlund et al. (2015) calculated the ratio of the forces directly (i.e., $\Lambda$) and compared various other definitions of the Elsasser number against this. The approach provided in Equation 7b could also be applied to other ratios of forces and may therefore aid the study of further regime transitions (e.g., Oruba & Dormy, 2014) of the Geodynamo or other natural dynamos, combined with other recent advances (Soderlund et al., 2025).

To obtain a balance representing the key dynamics in the bulk, it can be important to account for the contribution in boundary layers (Naskar et al., 2025; Schwaiger et al., 2019). Boundary layers have been removed in our calculations and we note that the curl operation removes the geostrophic flow, which is a large contribution to the axisymmmetric flow. Further details, including on the determination of lengthscale dependent solenoidal forces, $\hat{C}_l$, are presented in the Text S1 in Supporting Information S1.

### 3. Strong-Field Dynamos at Low $\mathcal{E}$

We have identified weak- and strong-field branches of dynamo action at all Ekman numbers tested (including at the very numerically demanding $\mathcal{E} = 10^{-6}$). On both the weak- and strong-field branch the surface magnetic field is dominantly dipolar (Figures 1a and 1b)), but the flow differs drastically on the two branches (Figures 1c and 1d)).

#### 3.1. Bifurcation Diagrams

We track bifurcations as a function of supercriticality, $\tilde{Ra}'$; the values of $Ra_c$ used in our study are listed in Table 1. Bifurcations in $(\tilde{Ra}', \Lambda')$-space are shown in Figure 2, where each panel displays averaged simulation data at a different value of $\mathcal{E}$. In each plot we find weak-field solutions (with small $\Lambda'$) above some critical value of $\tilde{Ra}'$ below which dynamo action is lost. However, such solutions only exist if $\mathcal{E}_m < \mathcal{E}_{mc}$ (or, equivalently, if $P_m > P_{mc}$ where $P_m = \nu/\kappa = \mathcal{E}/\mathcal{E}_m$ is the magnetic Prandtl number). Here, $P_{mc} = \mathcal{E}/\mathcal{E}_{mc}$ is the critical value for the onset of dynamo action; see Table 1, where values in parentheses are estimated from Figure 1 of Christensen and Aubert (2006). If $\mathcal{E}_m$ is small enough the weak-field solution becomes unstable (as $\tilde{Ra}'$ is increased) leading to a transition to a separate strong-field branch characterized by large $\Lambda'$, with bistability between the branches for an interval of $\tilde{Ra}'$. This behavior is not seen in the branches at larger $\mathcal{E}_m$ and its existence is used to identify the weak- and strong-field branches with the bistable interval shaded in Figure 2. It is therefore immediately apparent from the plots of Figure 2 that the general picture identified by Dormy (2016) remains valid as the Ekman number is reduced. That is, for a given $\mathcal{E}$, separate weak- and strong-field branches can only be found for small enough $\mathcal{E}_m$.

For each $\mathcal{E}$, it is possible to approximate a critical value of the magnetic Ekman number, $\mathcal{E}_{mc}^S$, below which the strong-field branch exists (or, equivalently, for $P_m > P_{mc}^S$). For the purposes of our numerical work we define $\mathcal{E}_{mc}^S$ to be the largest value of $\mathcal{E}_m$ for which bistability is found (see Table 1). Tracking of branches at further values of $\mathcal{E}_m$ would be required to improve the accuracy of these critical values.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Snapshots of the radial magnetic field, $B_r$, on a spherical surface at the outer boundary for: (a) a weak-field solution ($\tilde{Ra}' = 3, \mathcal{E} = 10^{-5}, \mathcal{E}_m = 2 \times 10^{-6}$); (b) a strong-field solution ($\tilde{Ra}' = 5, \mathcal{E} = 10^{-5}, \mathcal{E}_m = 2 \times 10^{-6}$). (c–d) As (a–b) but displaying isocontours of the vertical vorticity, $\omega_z = [\nabla \times \mathbf{u}]_z$, colored with the temperature field, $T$.</figcaption></figure>

For a given $\mathcal{E}$, both the separation between the branches in $\Lambda'$-space and the range of $\tilde{Ra}'$ for which bistability is observed reduces as $\mathcal{E}_m$ increases. This continues until the double turning point (see Dormy, 2025) connecting the individual branches unwinds and the branches merge at $\mathcal{E}_m \sim \mathcal{E}_{m_c}^S$. This marks the point in $\mathcal{E}_m$-space above which the strong-field branch no longer exists. These properties are supportive of the existence of a cusp singularity (Dormy, 2016) where the regions of bistability in Figure 2 become vanishingly small. Separately, bistability between the branches also moves to larger $\tilde{Ra}$ as $\mathcal{E}_m$ increases reflecting the increased forcing required to drive dynamos at higher magnetic diffusivity.

For the values of $\mathcal{E}$ tested, the strong-field branch was not observed for $\mathcal{E} < \mathcal{E}_m$ (i.e., for $P_m < 1$) although it does exist for $\mathcal{E} = \mathcal{E}_m = 10^{-6}$ (i.e., $P_m = 1$). Hence below this value of $\mathcal{E}$ the strong-field branch will presumably exist for $P_m < 1$, allowing for its preservation as both $\mathcal{E}_m$ and $P_m$ are reduced to small values appropriate for Earth's core. This further supports the idea of a distinguished limit (see Section 4) proposing $\mathcal{E}_m$ (or $P_m$) should scale with $\mathcal{E}$ to maintain strong-field solutions.

An unexpected feature of the results presented in Figure 2 is the reduction in $\Lambda'$ as $\mathcal{E}$ is lowered. Its value both at the lower end of the strong-field branch and for strong-field solutions more generally is found to drop below an $O(1)$ value (particularly so at $\mathcal{E} = 10^{-6}$). This unanticipated attribute prompts a reassessment of the definitions used for the Elsasser number and introduction of the Curlsasser number.

<table><caption>**Table 1**

Critical Values of Key Parameters for Values of the Ekman Number Considered in This Study</caption><thead><tr><th>ε</th><th>Ra<sub>c</sub></th><th>ε<sub>m c</sub></th><th>P<sub>m c</sub></th><th>ε<sub>m c</sub><sup>S</sup></th><th>P<sub>m c</sub><sup>S</sup></th></tr></thead><tbody><tr><td>3 × 10<sup>-4</sup></td><td>60.8</td><td>2 × 10<sup>-4</sup></td><td>1.5 (~ 1.5)</td><td>2.1 × 10<sup>-5</sup></td><td>14</td></tr><tr><td>10<sup>-4</sup></td><td>69.7</td><td>2 × 10<sup>-4</sup></td><td>0.5 (~ 0.45)</td><td>2 × 10<sup>-5</sup></td><td>5</td></tr><tr><td>3 × 10<sup>-5</sup></td><td>84.9</td><td>1.2 × 10<sup>-4</sup></td><td>0.25 (~ 0.18)</td><td>6 × 10<sup>-6</sup></td><td>5</td></tr><tr><td>10<sup>-5</sup></td><td>106</td><td>5 × 10<sup>-5</sup></td><td>0.2 (~ 0.09)</td><td>5 × 10<sup>-6</sup></td><td>2</td></tr><tr><td>10<sup>-6</sup></td><td>179</td><td>3.3 × 10<sup>-5</sup></td><td>0.03 (~ 0.045)</td><td>10<sup>-6</sup></td><td>1</td></tr></tbody></table>

*Note.* The (rotational) critical Rayleigh number, $Ra_c$, is for the onset of (non-magnetic) convection. No superscript indicates (approximate) maximal $\mathcal{E}_m$ values (minimal $P_m$) required for weak-field dipolar dynamo action. The superscript “S” indicates equivalent quantities for the strong-field dipolar branch.
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Modified Elsasser number, $\Lambda'$, as a function of supercriticality, $\tilde{Ra}'$, for (a) $\mathcal{E} = 10^{-4}$, (b) $\mathcal{E} = 3 \times 10^{-5}$, (c) $\mathcal{E} = 10^{-5}$, and (d) $\mathcal{E} = 10^{-6}$. Each point represents a single simulation where data has been volume averaged and time-averaged. Filled/unfilled circles respectively indicate strong-/weak-field dipolar; diamonds for fluctuating multipolar; cross symbols for failed dynamos (with data points moved to the lower horizontal axis). Shaded regions indicate observed bistability between the weak- and strong-field branches.</figcaption></figure>

Bifurcations in $(\tilde{Ra}', \sigma)$-space are shown in Figure 3. A key difference is found in the magnitude of $\sigma$ when compared with the values of $\Lambda'$ from Figure 2. Strong-field solutions cluster around $\sigma = 1$. Crucially this remains the case as $\mathcal{E}$ or $\mathcal{E}_m$ are lowered, in contrast to the behavior of $\Lambda'$. This suggests that $\sigma$ is a more robust measure of the ratio of the Lorentz to Coriolis effects that remains $O(1)$ as input parameters are moved to more realistic values. Naturally, the (non-)existence of separate branches and positions of bistability in $\tilde{Ra}'$-space are unchanged from Figure 2. However, separation between branches in $\sigma$-space differs from the equivalent in $\Lambda'$-space.

## 3.2. Strong-Field Master Curve

Data points for all values of Ekman number are presented in Figures 4a–4c, as function of $\tilde{Ra}$. Solutions on the strong-field branch are shown as filled symbols whereas all other regimes are displayed with unfilled symbols. Figure 4b is an updated version of Figure 9b from Dormy et al. (2018) now incorporating several values of the Ekman number, across several orders of magnitude. In contrast with that work, a curve of strong-field solutions
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Curl sasser number, $\sigma$, as a function of supercriticality, $\overline{Ra}'$, for (a) $\mathcal{E} = 10^{-4}$, (b) $\mathcal{E} = 3 \times 10^{-5}$, (c) $\mathcal{E} = 10^{-5}$, and (d) $\mathcal{E} = 10^{-6}$. See Figure 2 for explanations of symbols and shaded regions.</figcaption></figure>

with $\Lambda' \sim O(1)$ is no longer observed. As discussed earlier, this results from a reduction in $\Lambda'$-space of the location of the strong-field branch as $\mathcal{E}$ is reduced. Conversely, the Curl sasser number shown in Figure 4c captures the strong-field solutions on a striking master curve with $\sigma \sim 1$. The plots of Figures 4a–4c incorporate simulations across a range of all input parameters, further highlighting the robustness of $\sigma$ in identifying strong-field solutions.

Some solutions are “strong-ish” in the sense that they have $\Lambda'$ and/or $\sigma$ close to 1 but are not found above and beyond a turning point in the bifurcation diagrams of Figures 2 and 3. Such solutions are typically very close to the cusp singularity in parameter space, for which the bistability is replaced by a smooth transition (see discussion in Dormy, 2025). For these dynamos, the Lorentz force does not enter the primary balance and inertia and/or viscous effects remain important in the solution (see Section 3.3). Such solutions also destabilize to a multipolar regime at larger convective driving, unlike the genuine strong-branch solutions (Menu et al., 2020).

### 3.3. Force Balance

Solutions on the strong-field branch have previously been shown to be in MAC-balance or MC-balance across a range of lengthscales (Teed & Dormy, 2023) and broadly by position in the domain (Dormy, 2016). Some solutions from Schwaiger et al. (2019) are also suggestive of a large-scale MAC balance although distinct branches
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Dependence on the Rayleigh number, $\tilde{Ra}$, of (a) $\Lambda$, the Elsasser number (directly from forces), (b) $\Lambda'$, the modified Elsasser number, and (c) $\sigma$, the Curlsasser number (directly from solenoidal forces). (d) Location of strong-field solutions in Ekman number space. Blue lines show extrema for possible scalings of $\mathcal{E}_m$ with $\mathcal{E}$ for preservation of strong-field solutions; the red line is a proposed scaling for the onset of the weak-field branch (Christensen & Aubert, 2006). See Figure 2 for explanations of symbols.</figcaption></figure>

of dynamo action are not identified; they instead use $M$ (the ratio of the magnetic to kinetic energy) as a proxy to determine the influence of the magnetic field. Values of $M$ from our solutions on the strong-field branch are of a similar order (Teed & Dormy, 2025). Furthermore, since we consider solenoidal forces here, the purely geostrophic contribution, along with all gradient parts, is filtered out.

Figure 5 shows how MAC-balance (at the largest lengthscales) and MC-balance (at intermediate to small lengthscales) is preserved on the strong-field branch as $\mathcal{E}$ is reduced. Figures 5a–5b are for simulations chosen at $\mathcal{E} = 10^{-5}$ and $\mathcal{E} = 10^{-6}$ to demonstrate the balance at the lowest Ekman numbers available. Similar results are found at larger $\mathcal{E}$; indeed, an equivalent plot for $\mathcal{E} = 10^{-4}$ can be found in Figure 3e of Teed and Dormy (2023). Conversely, solutions that are “strong-ish” (Figure 5c) have a diminished role for the Lorentz force, no fully formed MAC-balance at large scales, and have a greatly increased role for inertial and/or viscous forces.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Solenoidal forces, $\hat{C}_l$, as a function of lengthscale, $l$ for simulations with: (a) $\mathcal{E} = 10^{-5}$, $\mathcal{E}_m = 5 \times 10^{-6}$, and $\overline{\mathcal{R}a'} = 10$ (strong-field); (b) $\mathcal{E} = 10^{-6}$, $\mathcal{E}_m = 10^{-6}$, and $\overline{\mathcal{R}a'} = 10$ (strong-field); (c) $\mathcal{E} = 10^{-4}$, $\mathcal{E}_m = 10^{-4}$, and $\overline{\mathcal{R}a'} = 10$ ("strong-ish" field).</figcaption></figure>

## 4. Scaling of the Strong-Field Branch With $\mathcal{E}$ for a Distinguished Limit

Dormy (2016) suggested that, since $\mathcal{E}$ and $\mathcal{E}_m$ are both small parameters in Earth's core, they could be related in numerical studies via a distinguished limit, involving only one small parameter, $\epsilon$. This ensures that as $\mathcal{E}$ and $\mathcal{E}_m$ tend to zero, they do so in such a way that the nature of the solution (e.g., its strong-field property), is preserved in the limiting process. Following Dormy (2016), we therefore set $\mathcal{E} = \mathcal{E}_0\epsilon$ and $\mathcal{E}_m = \mathcal{E}_{m0}\epsilon^\alpha$ for some constants $\mathcal{E}_0$, $\mathcal{E}_{m0}$, and $\alpha$, to be determined. The values of $\mathcal{E}_0$ and $\mathcal{E}_{m0}$ are selected to set the desired regime in parameter space and $\alpha$ sets the scaling relationship to preserve the behavior as $\epsilon \to 0$. Alternatively, one can pose the same scaling using $P_m$, in which case $P_m = \mathcal{E}/\mathcal{E}_m = P_{m0}\epsilon^{(1-\alpha)}$ where $P_{m0} = \mathcal{E}_0/\mathcal{E}_{m0}$.

Unlike Dormy (2016), the data presented in our study uses values of $\mathcal{E}$ ranging several orders of magnitude. It can therefore provide a numerical estimation of $\alpha$ for the preservation of strong-field solutions with MAC-balance in the limit of small $\mathcal{E}$. Assuming that our values of $\mathcal{E}$ and $\mathcal{E}_m$ are approaching some asymptotic behavior, the values of $\mathcal{E}_{m0}^S$ (Table 1) offer a range of admissible exponents, $\alpha \in [0.35, 0.8]$, as shown in Figure 4d. This range includes the value $\alpha = 1/2$ used by Aubert et al. (2017). Yet, strikingly, the large values of $\alpha$ in this range are associated with a fast decrease of $\mathcal{E}_m$ with $\mathcal{E}$. At an approximate core value of $\mathcal{E} \sim 10^{-15}$, such $\alpha$ would yield a minimum $\mathcal{E}_m$ value that is too small for Earth's core where $\mathcal{E}_m \sim 10^{-9}$. This points to the importance of further pursuing this (very demanding) parameter survey.

## 5. Discussion

We have tracked the weak-field and strong-field branches as the effect of rotation is increased (via decreasing Ekman number). The emergence (or not) of the strong-field branch is dependent on the value of the magnetic Ekman number. A sufficiently low $\mathcal{E}_m$ (or high $P_m$) is required for the strong-field branch to be realized. The critical value of the magnetic Ekman number, $\mathcal{E}_{m0}$, below which strong-field solutions are found has been found to scale with the Ekman number as $\mathcal{E}_m \propto \mathcal{E}^\alpha$ where $\alpha \in [0.35, 0.8]$. The value of $\alpha$ could be refined further with more simulation data and, while its value may depend weakly on the precise mathematical set-up (e.g., magnetic or thermal boundary conditions), the qualitative picture of Figure 4 should remain intact. On the strong-field branch, one would expect an approximately $O(1)$ value for the Elsasser number given the balance of forces is MAC at large lengthscales and MC at a range of intermediate and small scales (Schwaiger et al., 2019; Teed & Dormy, 2023). Forming the ratio of the solenoidal forces (i.e., curls of forces) accurately captures the effects of the Lorentz force across a range of rotation rates. We therefore advocate for a “Curlsasser number” as a new measure for identifying solutions as strong-field candidates.

The motivation for investigating the strong-field solutions is their replication of the expected force hierarchy of Earth's core. However, along and between individual strong-field branches, other important geophysical properties naturally evolve. For instance, $\mathcal{R}_m$ typically varies from $O(10^2)$ at onset (at low driving) to the more Earth-like $O(10^3)$ (at larger values of driving). Similarly, although strong-field solutions are dipolar, the level of dipolarity on the core-mantle boundary can vary considerably. Strong-field solutions typically maintain dipolarity
values closer to those of present-day Earth than their weak-field counterparts (Teed & Dormy, 2025). As a result, strong-field solutions are more likely to exhibit features such as reversed magnetic flux patches, yet full dipole reversals are absent in our solutions in the time-frame we have studied. Further work studying strong-field solutions across parameter space is obviously needed to identify the most geophysically suitable branches; for example, Jones and Tsang (2025) recently suggested varying *q* may be a fruitful endeavor for obtaining solutions with both strong dipolar fields and Earth-like dipole reversals. In this work, we have cemented the importance of the strong-field branch (and bistability between dipolar branches) and demonstrated its persistence and scaling as input parameters are gradually moved toward realistic values for the Geodynamo.

## Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

## Data Availability Statement

The data set for this research is available on the University of Glasgow's "Enlighten: Research Data" repository (Teed & Dormy, 2025).

## Acknowledgments

R.J.T. acknowledges support from the UK Science & Technology Facilities Council under Grant ST/Y00146X/1. R.J.T. also acknowledges support from DiRAC Project ACTP245. This work used the DiRAC Memory Intensive service (Cosma8) at Durham University, managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). 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. This work was also granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. The authors would like to thank two anonymous referees for their useful sets of comments, which improved the manuscript.

## References

Aubert, J., Gastine, T., & Fournier, A. (2017). Spherical convective dynamos in the rapidly rotating asymptotic regime. *Journal of Fluid Mechanics*, 813, 558-593. https://doi.org/10.1017/jfm.2016.789

Buffett, B. (2000). Earth's core and the geodynamo. *Science*, 288(5473), 2007–2012. https://doi.org/10.1126/science.288.5473.2007

Christensen, U., & Aubert, J. (2006). Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields. *Geophysical Journal International*, 166(1), 97–114. https://doi.org/10.1111/j.1365-246x.2006.03009.x

Christensen, U., Olson, P., & Glatzmaier, G. (1999). Numerical modelling of the geodynamo: A systematic parameter study. *Geophysical Journal International*, 138(2), 393–409. https://doi.org/10.1046/j.1365-246x.1999.00886.x

Dormy, E. (2016). Strong-field spherical dynamics. *Journal of Fluid Mechanics*, 789, 500–513. https://doi.org/10.1017/jfm.2015.747

Dormy, E. (2025). Rapidly rotating magnetohydrodynamics and the geodynamo. *Annual Review of Fluid Mechanics*, 57(1), 335–362. https://doi.org/10.1146/annurev-fluid-031224-121649

Dormy, E., Oruba, L., & Petitdemange, L. (2018). Three branches of dynamo action. *Fluid Dynamics Research*, 50(1), 011415. https://doi.org/10.1088/1873-7005/aa769c

Guseva, A., Petitdemange, L., & Tobias, S. (2025). Run-away transition to turbulent strong-field dynamo. *Journal of Geophysical Research: Planets*, 130(1), e2024JE008496. https://doi.org/10.1029/2024je008496

Horn, S., & Aurnou, J. M. (2022). The Elbert range of magnetostrophic convection. I. Linear theory. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 478(2264), 20220313. https://doi.org/10.1098/rspa.2022.0313

Horn, S., & Aurnou, J. M. (2025). The Elbert range of magnetostrophic convection. II. Comparing linear theory to nonlinear low-Rm simulations. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 481(2310), 20240016. https://doi.org/10.1098/rspa.2024.0016

Jones, C. A., & Tsang, Y.-K. (2025). Low inertia reversing geodynamos. *Physics of the Earth and Planetary Interiors*, 360, 107303. https://doi.org/10.1016/j.pepi.2024.107303

Kutzner, C., & Christensen, U. (2002). From stable dipolar towards reversing numerical dynamos. *Physics of the Earth and Planetary Interiors*, 131(1), 29–45. https://doi.org/10.1016/s0031-9201(02)00016-x

Mannix, P. M., Ponty, Y., & Marcotte, F. (2022). Systematic route to subcritical dynamo branches. *Physical Review Letters*, 129(2), 024502. https://doi.org/10.1103/physrevlett.129.024502

Mason, S. J., Guervilly, C., & Sarson, G. R. (2022). Magnetoconvection in a rotating spherical shell in the presence of a uniform axial magnetic field. *Geophysical & Astrophysical Fluid Dynamics*, 116(5-6), 458–498. https://doi.org/10.1080/03091929.2022.2107202

Menu, M., Petitdemange, L., & Galtier, S. (2020). Magnetic effects on fields morphologies and reversals in geodynamo simulations. *Physics of the Earth and Planetary Interiors*, 307, 106542. https://doi.org/10.1016/j.pepi.2020.106542

Morin, V., & Dormy, E. (2009). The dynamo bifurcation in rotating spherical shells. *International Journal of Modern Physics B*, 23(28n29), 5467–5482. https://doi.org/10.1142/s021797920906378x

Naskar, S., Davies, C., Mound, J., & Clarke, A. (2025). Force balances in spherical shell rotating convection. *Journal of Fluid Mechanics*, 1009, A70. https://doi.org/10.1017/jfm.2025.259

Oruba, L., & Dormy, E. (2014). Transition between viscous dipolar and inertial multipolar dynamos. *Geophysical Research Letters*, 41(20), 7115–7120. https://doi.org/10.1002/2014gl062069

Schwaiger, T., Gastine, T., & Aubert, J. (2019). Force balance in numerical geodynamo simulations: A systematic study. *Geophysical Journal International*, 219(Supplement_1), S101–S114. https://doi.org/10.1093/gji/ggz192

Soderlund, K., King, E., & Aurnou, J. (2012). The influence of magnetic fields in planetary dynamo models. *Earth and Planetary Science Letters*, 333, 9–20. https://doi.org/10.1016/j.epsl.2012.03.038

Soderlund, K., Sheyko, A., King, E., & Aurnou, J. (2015). The competition between Lorentz and Coriolis forces in planetary dynamos. *Progress in Earth and Planetary Science*, 2(1), 1–10. https://doi.org/10.1186/s40645-015-0054-5

Soderlund, K., Wulff, P., Käpylä, P., & Aurnou, J. (2025). Magnetohydrodynamic control of differential rotation and dynamo transitions: Rise of the local magnetic Rossby number. *Monthly Notices of the Royal Astronomical Society*, 541(2), 1816–1830. https://doi.org/10.1093/mnras/staf1081

Teed, R., & Dormy, E. (2023). Solenoidal force balances in numerical dynamos. *Journal of Fluid Mechanics*, 964, A26. https://doi.org/10.1017/jfm.2023.332
Teed, R., & Dormy, E. (2025). Scaling of strong-field spherical dynamos [Dataset]. University of Glasgow "Enlighten: Research Data" repository. https://doi.org/10.5525/gla.researchdata.2006

Willis, A., Sreenivasan, B., & Gubbins, D. (2007). Thermal core–mantle interaction: Exploring regimes for “locked” dynamo action. *Physics of the Earth and Planetary Interiors*, 165(1–2), 83–92. https://doi.org/10.1016/j.pepi.2007.08.002