# Geomagnetic dipole stability and zonal flows controlled by mantle heat flux heterogeneities

Thomas Frasson, Nathanaël Schaeffer, Henri-Claude Nataf, Stéphane Labrosse

► To cite this version:

Thomas Frasson, Nathanaël Schaeffer, Henri-Claude Nataf, Stéphane Labrosse. Geomagnetic dipole stability and zonal flows controlled by mantle heat flux heterogeneities. Geophysical Journal International, 2025, 240, pp.1481-1504. 10.1093/gji/ggae457. hal-04780728

HAL Id: hal-04780728

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

Submitted on 21 Mar 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Geomagnetic dipole stability and zonal flows controlled by mantle heat flux heterogeneities

T. Frasson<sup>1</sup>, N. Schaeffer<sup>1</sup>, H-C. Nataf<sup>1</sup> and S. Labrosse<sup>1</sup>

Accepted 2024 December 18. Received 2024 November 18; in original form 2024 June 20

## SUMMARY

Palaeomagnetic evidence shows that the behaviour of the geodynamo has changed during ge- ological times. These changes are visible through variations in the strength and stability of the magnetic dipole. Variations in the heat flux at the core–mantle boundary (CMB) due to mantle convection have been suggested as one possible mechanism capable of driving such a change of behaviour. This work aims at acquiring a more complete understanding of how lateral het- erogeneities of the CMB heat flux affect the geodynamo while other relevant parameters are pushed towards realistic values. For this purpose, we ran geodynamo simulations with degree 1 and 2 spherical harmonic patterns of heat flux at the CMB. Several geodynamo models are used, ranging from standard numerical dynamos to more extreme parameters, including strong field cases and turbulent cases. We show that heat flux heterogeneities with amplitudes compatible with our knowledge of mantle convection history can favour multipolar dynamos. The multipolar transition is associated with a disruption of westward flows either through eastward thermal winds or through a loss of equatorial symmetry. Strong field dynamo models are found to have larger westward flows and are less sensitive to heat flux heterogeneities. Fur- thermore, we find that the dipolar fraction of the magnetic field correlates with $M_{Za}^* = \frac{\Lambda_{Za}}{Rm_{Za}^2}$ where $\Lambda_{Za}$ is the zonal antisymmetric Elsasser number and $Rm_{Za}$ is the zonal antisymmetric magnetic Reynolds number. Importantly, $M_{Za}^*$ estimated for the Earth’s core is consistent with a reversing dipolar magnetic field. Within the range of $M_{Za}^*$ susceptible to reversals, breaking the equatorial symmetry or forcing eastward zonal flows through an equatorial cooling of the core consistently triggers reversals or a transition towards multipolar dynamos in our simu- lations. Our results support that time variations of heat-flux heterogeneities driven by mantle convection through Earth’s history are capable of inducing the significant variations in the reversal frequency observed in the palaeomagnetic record.

**Key words:** Dynamo: theories and simulations; Reversals: process, time scale, magne- tostratigraphy; Numerical modelling; Core–mantle boundary; Heat flow.

# 1 INTRODUCTION

## 1.1 Magnetic reversals and dipolar–multipolar transition in dynamo models

The magnetic field of the Earth is dominated by its dipole com- ponent. This magnetic dipole undergoes chaotic polarity rever- sals, with a reversal frequency that varies through geological times (Lowrie & Kent 2013). Reversals of the magnetic dipole have been reproduced in numerical dynamo models since the first attempts to simulate the geodynamo (Glatzmaier & Roberts 1995). These reversals are characterized by a transient decrease in the strength and stability of the magnetic dipole, which regrows and stabilizes

with the opposite polarity (Olson *et al.* 2011). Numerical mod- els of reversing dynamos can be found at the transition between strongly dipolar dynamos and weakly dipolar (or multipolar) dy- namos (Kutzner & Christensen 2002; Olson & Christensen 2006; Wicht & Tilgner 2010). In the available set of dynamo models, this dipolar–multipolar transition seems to be controlled by the relative importance of inertia in the force balance (Christensen & Aubert 2006; Sreenivasan & Jones 2006). The strength of inertial forces relative to Coriolis forces has notably been found to successfully grasp this transition in early dynamo models, larger inertia being associated with multipolar dynamos (Christensen & Aubert 2006; Olson & Christensen 2006; Christensen 2010; Wicht & Tilgner 2010). In dynamo simulations operating in a strong magnetic field

<footnote n="1"><sup>1</sup>Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERre, F-38000 Grenoble, France. E-mail: t.a.frasson@leeds.ac.uk</footnote>

<footnote n="2"><sup>2</sup>ENS de Lyon, Université Lyon-1, LGL-TPE, 46 allée d'Italie, F-69007 Lyon, France</footnote>
regime, in which the Lorentz force plays an important role as ex-
pected for the Earth, the relative strength of Coriolis forces and
inertia does not accurately grasp the transition (Menu et al. 2020;
Tassin et al. 2021). Tassin et al. (2021) showed that the dipolar–
multipolar transition can be more accurately described by the ratio
of inertia to Lorentz forces, approximated by the kinetic over mag-
netic energy ratio. According to this criterion, reversing dynamos
should be obtained when the magnetic and kinetic energies within
the core are of similar magnitudes. However, the Earth’s dynamo
seems to operate differently, with a magnetic energy several orders
of magnitude higher than the kinetic energy. Though reversing dy-
namos producing a higher magnetic energy than kinetic energy on
average have been obtained (Driscoll & Olson 2009a; Nakagawa &
Davies 2022), inertial forces are less than one order of magnitude
weaker than Lorentz forces and tend to strengthen during reversals,
suggesting that inertia is probably still important to drive reversals
in these simulations.

## 1.2 Mantle control on the behaviour of the magnetic dipole

Based on numerical dynamo models, different mechanisms have
been proposed to explain the observed variations in the reversal
frequency of the magnetic dipole. Driscoll & Olson (2009a), for
example, suggest that variations in the control parameters of the
geodynamo, such as the rotation rate or the power input, could ex-
plain regime changes. This mechanism requires the Earth to have
been operating close to the inertia-driven dipolar-multipolar tran-
sition since the onset of dynamo action in the Earth's core. Using
these assumptions, Driscoll & Olson (2009b) showed that variable
reversal frequencies can be obtained in dynamo models with a sim-
ple core evolution model.

In this regard, the heat flux at the core–mantle boundary (CMB)
can have an important role by modulating convective activity at the
top of the core. Moreover, the CMB heat flux is expected to be
strongly heterogeneous (Nakagawa & Tackley 2008). Large-scale
heat flux heterogeneities have been shown to have a significant im-
pact on the behaviour of the dynamo (Glatzmaier et al. 1999; Olson
& Christensen 2002; Kutzner & Christensen 2004; Olson et al.
2010; Olson & Amit 2014; Amit et al. 2015; Sahoo et al. 2016).
However, no clear trend has emerged. For instance, large heat flux
heterogeneities have been shown to suppress dynamo action (Olson
& Christensen 2002), while weak heterogeneities could on the con-
trary help sustain a dynamo (Sahoo et al. 2016). Interestingly, the
reversal frequency has been shown to be very sensitive to heat flux
heterogeneities. Various effects have been reported so far. Notably,
large heat flux heterogeneities have been found to increase the rever-
sal frequency, particularly for high equatorial heat fluxes (Kutzner &
Christensen 2004; Olson et al. 2010). Stabilizing effects of the mag-
netic dipole due to heat flux heterogeneities are also observed when
polar cooling exceeds equatorial cooling (Glatzmaier et al. 1999;
Kutzner & Christensen 2004).

Recent models of the geodynamo showed that the present-day
heat flux pattern, as inferred from seismic tomography, can explain
the longitudinal structure of the observed magnetic field (Mound
& Davies 2023), and could trigger reversals from locally high heat
flux regions (Terra-Nova & Amit 2024). These findings suggest that
the present-day heat flux pattern at the CMB plays an important role
in the dynamics of the core. Though the present-day lateral hetero-
geneities of temperatures in the lower mantle are known to a certain
extent (Trampert 2004; Mosca et al. 2012), the past CMB heat flux

is highly uncertain. On large scales, the present-day CMB heat flux
is controlled by the positions of the observed large low velocity
provinces (LLVPs) below Africa and the Pacific (Garnero & McNa-
mara 2008). These structures are interpreted as regions hotter than
the surrounding mantle, thus locally extracting a low heat flux from
the core (Nakagawa & Tackley 2008). Long-term stability of the ob-
served LLVPs have been suggested based on records of hotspot erup-
tion sites (Burke et al. 2008; Dziewonski et al. 2010; Torsvik et al.
2010), suggesting that the present-day heat flux pattern could be rel-
evant for at least the past 300 Myr. However, mantle convection mod-
els show that basal mantle structures are pushed by subducting slabs
and thus could have significantly been displaced through geological
times (Zhang & Zhong 2011; Flament et al. 2022). These variations
in the CMB heat flux are expected to occur on timescales similar
to the longest timescales of variations in the magnetic field, corre-
sponding to changes in the reversal frequency (Biggin et al. 2012;
Hounslow et al. 2018). Understanding how the dynamo can react to a
change in the CMB heat flux geometry is thus of importance to con-
strain the present-day geodynamo and its evolution through Earth’s
history.

Mantle convection models can be used to have an estimate of
the CMB heat flux heterogeneities in the past (Zhang & Zhong
2011). Using self-consistent and plate-driven mantle convection
models that reproduce ~ 1 Gyr of mantle convection, Frasson et al.
(2024) showed how episodes of true polar wander could signifi-
cantly alter the CMB heat flux distribution on timescales that can
be shorter than the convective timescale in the mantle. The strength
of equatorial and polar cooling of the core is notably shown to
depend on the way true polar wander rotates chemical structures
at the base of the mantle. The mantle convection scenarios cor-
rected from true polar wander presented in Frasson et al. (2024)
can thus be used to obtain estimates of heat flux heterogeneity
amplitudes.

## 1.3 Aim of the study

Our current understanding of how heat flux heterogeneities affect
the geodynamo is limited by the parameter regime reached in numer-
ical models. In the last decade, new geodynamo simulations have
reached a rapidly rotating turbulent regime, in which the Lorentz
force dominates inertial and viscous forces as expected for the Earth
(Yadav et al. 2016; Aubert et al. 2017; Schaeffer et al. 2017). Dy-
namos operating in this regime are computationally demanding,
making parametric studies challenging. This challenge becomes
even more important when rare events such as reversals are tar-
geted. For this reason, previous studies used dynamo models that
are only moderately turbulent. Inertia and viscous forces have been
shown to be non-negligible in these types of dynamos (Soderlund
et al. 2012; Dormy 2016), which are therefore not representative of
the Earth's core.

In this study, we aim at exploring the effect of heat flux het-
erogeneities on dynamo models that extend from standard dynamo
regimes towards a more Earth-like rapidly rotating turbulent regime.
We first study the effect of large-scale heat flux patterns (spher-
ical harmonic degrees 1 and 2) on a standard dynamo reference
case. We then focus on the effect of polar cooling and equato-
rial cooling on more turbulent dynamos, exploring more Earth-
like parameter regimes. We notably focus on the reversing be-
haviour and the dipolar-multipolar transition triggered by heat flux
heterogeneities.
## 2 METHODS

### 2.1 Numerical model of the geodynamo

The Earth's liquid outer core is modelled as a spherical shell of inner radius $r_i$ and outer radius $r_o$. The thickness of the shell is defined as $D = r_o - r_i$. The aspect ratio is chosen to match the present-day Earth by fixing $r_i/r_o = 0.35$. The shell spins along the z axis at a rate $\Omega$. The acceleration of gravity is defined as $g = -g_0 \frac{r}{r_o} \mathbf{e}_r$ with $g_0$ the acceleration of gravity at the core surface and $\mathbf{e}_r$ the unit vector in the radial direction. We use the Boussinesq approximation in which the electrically conducting liquid contained in the shell is treated as an incompressible fluid of density $\rho_0$ and viscosity $\nu$. The electrical properties of the fluid are defined by its magnetic permeability $\mu_0$ and its magnetic diffusivity $\eta$. Convection in the Earth's outer core is the result of both thermal and compositional heterogeneities, which translate into density heterogeneities. Here, we use the codensity approximation defined in Braginsky & Roberts (1995), which is exact when both thermal and compositional diffusivities are equal, and when the same boundary conditions apply to temperature and composition fields. The thermal and compositional sources of buoyancy are then combined into one single variable called codensity, with diffusivity $\kappa$. If we call $T$ the temperature field, $\chi$ the concentration in light element and $\rho_0$ the reference density, the codensity is defined as $C = \rho_0(\alpha T + \beta \chi)$ with $\alpha = -\frac{1}{\rho_0} \frac{\partial \rho}{\partial T}$ and $\beta = -\frac{1}{\rho_0} \frac{\partial \rho}{\partial \chi}$ the thermal and compositional expansion coefficients, respectively. Convection in the shell is driven by imposed buoyancy fluxes at the inner and outer boundaries. The buoyancy fluxes per unit surface (in kg m$^{-2}$ s$^{-1}$) at the inner and outer boundaries are called $q_i$ and $q_o$, respectively, and are defined by

$$q_i = -\kappa \partial_r C(r_i), \quad (1)$$

$$q_o = -\kappa \partial_r C(r_o). \quad (2)$$

The buoyancy flux entering the core at the inner boundary $F_i$ and the buoyancy flux leaving the core at the outer boundary $F_o$ can then be obtained by integration over the boundaries:

$$F_i = \int_{S_i} q_i \, dS_i, \quad (3)$$

$$F_o = \int_{S_o} q_o \, dS_o, \quad (4)$$

where $S_i$ and $S_o$ are the inner and outer boundaries, respectively.

The geodynamo equations consists in the coupling of Navier-Stokes and induction equations. We solve this set of equations in their dimensionless form using the Boussinesq approximation. The thickness of the shell $D$ is used as the length scale. The time is scaled by the viscous timescale $\tau_v = \frac{D^2}{\nu}$ and the magnetic field by $\sqrt{\rho_0 \Omega \eta \mu_0}$. The codensity is scaled by $\frac{F_i}{4\pi D \nu}$. The codensity field $C$ is decomposed into a radially dependent base profile $C_0$ and a perturbation $C'$, such that $C = C_0 + C'$. Calling $\mathbf{u}$ the velocity field, $P$ the pressure and $\mathbf{B}$ the magnetic field, the set of dimensionless equations is as follow:

$$\partial_t \mathbf{u} + (\nabla \times \mathbf{u}) \times \mathbf{u} + \frac{2}{E} \mathbf{e}_z \times \mathbf{u} = -\nabla P + \nabla^2 \mathbf{u} \\ + \frac{1}{Pm} \frac{1}{E} (\nabla \times \mathbf{B}) \times \mathbf{B} + \frac{Ra}{Pr} C' \frac{r}{r_o} \mathbf{e}_r, \quad (5)$$

$$\partial_t C' + \mathbf{u} \cdot \nabla (C' + C_0) = \frac{1}{Pr} \nabla^2 (C' + C_0) + S_c, \quad (6)$$

$$\partial_t \mathbf{B} = \nabla \times (\mathbf{u} \times \mathbf{B}) + \frac{1}{Pm} \nabla^2 \mathbf{B}, \quad (7)$$

$$\nabla \cdot \mathbf{u} = 0 \quad (8)$$

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

where $S_c$ is a sink term to compensate for the imbalance between $F_i$ and $F_o$. The dimensionless parameters $E, Pm, Pr$ and $Ra$ are, respectively, the Ekman number, the magnetic Prandtl number, the Prandtl number and the Rayleigh number. They are defined using the physical parameters of the system by

$$E = \frac{\nu}{\Omega D^2}, \quad (10)$$

$$Pm = \frac{\nu}{\eta}, \quad (11)$$

$$Pr = \frac{\nu}{\kappa}, \quad (12)$$

$$Ra = \frac{g_0 F_i D^2}{4\pi \kappa \nu^2 \rho_0}. \quad (13)$$

No slip conditions are applied at both boundaries for the velocity. The outer boundary is electrically insulating in all our simulations. The inner core is treated as a conducting sphere in all but one of our dynamo models, with a conductivity of the inner core equal to that of the outer core. An insulating inner core is used for the most numerically demanding simulations. A heterogeneous flux condition at the outer boundary can be imposed. This flux pattern is described by an amplitude $\delta q_o$ and a pattern $p(\theta, \phi)$, which depends on the colatitude $\theta$ and the longitude $\phi$. We impose that the heterogeneous flux has a zero space average and set $F_o = 0$.

In this codensity formalism, the flux imposed at either of the boundaries is a non-distinguishable combination of a heat flux and a light-element flux. This buoyancy flux is expected to be of thermal origin at the top of the core, with a null compositional flux, and mostly of compositional origin at the inner boundary. In the following, we will thus refer to the flux at the top of the core as a heat flux. As we use the Boussinesq framework, the imposed fluxes at the boundaries correspond to superisentropic fluxes, and the isentropic flux would have to be added to obtain the total fluxes entering and leaving the core. This means notably that a negative flux at the top of the core tends to act against convection by stratifying the fluid.

The equations are solved using the XSHELLS code benchmarked in Matsui et al. (2016) and already used in previous studies (Schaeffer et al. 2017; Guervilly et al. 2019). Finite differences are used in the radial direction, while the fields are expanded in spherical harmonics in the horizontal directions. The spherical harmonics transforms are performed using SHTns (Schaeffer 2013). We use hyperdiffusivity to damp small scales and stabilize the simulations in some cases as done in Guervilly et al. (2019) (see Appendix A for more details).

### 2.2 Output parameters of interest

The simulation outputs are given in the following as a function of time scaled by the magnetic dipole diffusion timescale rather than the viscous timescale used to solve the equations. The dipole diffusion timescale is defined by

$$\tau_{\eta} = \frac{1}{\pi^2} \frac{r_o^2}{\eta} \quad (14)$$

and can be estimated to be approximately 40 kyr for the Earth using $\eta = 1 \text{ m}^2 \text{ s}^{-1}$ (Olson 2015), though the exact value of $\eta$ is prone to some debate (Pozzo et al. 2014; Holdenried-Chernoff & Buffett
2022). The conversion between the simulation time t and the output time t<sub>η</sub> is then simply t<sub>η</sub> = t (τ<sub>v</sub>/τ<sub>η</sub>).

Several output parameters are used to compare our simulations between each other and to the Earth. The strength of the magnetic dipole is quantified by the dipolar fraction f<sub>dip</sub> defined at the CMB as

$$
f_{\text{dip}}(t) = \sqrt{\frac{\sum_{m=0}^{1} (2 - \delta_{m0}) |b_r(l=1, m, t)|^2}{\sum_{l,m \ge 0} (2 - \delta_{m0}) |b_r(l, m, t)|^2}}, \quad (15)
$$

where *l* and *m* are the spherical harmonic degree and order, *b*<sub>*r*</sub>(*l*, *m*) is the spherical harmonic coefficient of degree *l* and order *m* of the radial component of the magnetic field at the CMB, and *δ*<sub>*lm*</sub> is the Kronecker delta. In order to compare our results to the Earth's magnetic field, we also compute an alternative dipolar fraction based only on the degrees *l* ≤ 12 of the magnetic field:

$$
f_{\text{dip}}^{\dagger}(t) = \sqrt{\frac{\sum_{m=0}^{1} (2 - \delta_{m0}) b(l=1, m, t)^2}{\sum_{l,m \ge 0} (2 - \delta_{m0}) b(l \le 12, m, t)^2}}. \quad (16)
$$

This definition of the dipolar fraction enables to compare with mea-
surements of the magnetic field that are limited to a maximum
spherical harmonic degrees *l* < 13 and is commonly used as out-
puts in dynamo simulations (Christensen & Aubert 2006). The di-
mensionless kinetic and magnetic energy integrated over the whole
spherical shell are, respectively, defined as

$$
E_u(t) = \frac{1}{2} \int_{V_f} \mathbf{u}(t)^2 \mathrm{d}V, \qquad (17)
$$

$$
E_b(t) = \frac{1}{2} \int_{V_c} \left( \frac{\mathbf{B}(t)}{\sqrt{E P m}} \right)^2 \mathrm{d}V, \quad (18)
$$

where $V_f$ is the volume of the fluid outer core and $V_c$ the volume of the conducting shell (corresponding to the full core if the inner core is conducting, or only the fluid outer core if the inner core is insulating). The ratio of these two quantities is called $M$ in the following, and is defined as

$$
M = E_b / E_u. \tag{19}
$$

This ratio is used as a proxy of the ratio between the Lorentz force and the inertial forces by Tassin *et al.* (2021), and is found to discriminate between dipole-dominated dynamos and multipolar dynamos, with a transition around *M* ≃ 1. The amplitudes of the velocity field and of the magnetic field are quantified by the magnetic Reynolds number *Rm* = *U*<sub>0</sub>*D*/*η* and the Elsasser number *Λ* = *B*<sup>2</sup>*ρ*<sub>0</sub>*Ω*<sub>μ0</sub>*η* where *U* and *B* are the flow and magnetic field RMS amplitudes, respectively. Using the dimensionless outputs of the simulation, those numbers are computed as

$$
Rm(t) = Pm \sqrt{\frac{2E_u(t)}{V_f}}, \qquad (20)
$$

$$
\Lambda(t) = \frac{2PmE}{V_f} E_b(t). \qquad (21)
$$

The strength of inertial force relative to Coriolis force can be ap-
proximated using the local Rossby number (Christensen & Aubert
2006)

$$
Ro_l(t) = \frac{E l_c(t)}{\pi} \sqrt{\frac{2E_u(t)}{V_f}}, \quad (22)
$$

where $l_c$ is the characteristic spherical harmonic degree of the flow.
This characteristic degree is defined by

$$
l_c(t) = \frac{\sum_l l E_{u,l}(t)}{E_u(t)} \qquad (23)
$$

with E<sub>u,l</sub> the kinetic energy at degree *l*. Finally, we quantify the stability of the magnetic dipole by computing the average period spent by the dipole in a given hemisphere. This period, called T<sub>rev</sub>, is computed for the reversing dynamos as

$$
T_{\text{rev}} = \frac{\Delta t_{\text{sim}}}{N_{\text{cross}}}, \qquad (24)
$$

where Δt<sub>sim</sub> is the duration of the equilibrated simulation and N<sub>cross</sub> is the number of times the magnetic dipole crosses the equator. We note that the value of T<sub>rev</sub> is not very robust when N<sub>cross</sub> is small, which can happen for short simulations or when the reversal rate is low.

Because the simulations start from an initial condition that is
different from their statistically steady state, the transient evolu-
tion from the initial condition has to be removed before computing
statistics on the diagnostics. We ensure that this transitional state is
not affecting the results by removing the evolution corresponding
to the first two dipole diffusion times.

2.3 Explored parameter space

2.3.1 Reference dynamo models

We select five reference dynamo cases with homogeneous heat flux
conditions at the top of the core, to which heterogeneous heat flux
patterns are applied. In Fig. 1 these dynamos are compared in the
E − Pm parameter space with previous studies that focused on the
effect of heterogeneous heat flux on the dynamo. Our reference
case E1e-4.hRm is comparable to dynamo models already used
in previous studies, with E = 10⁻⁴ and Pm = 3. Our cases E1e-
5.hRm, E1e-5_SF1 and E1e-5_SF2 have similar magnetic Prandtl
number, but have a lower viscosity than previous studies. Finally,
case E1e-6.hRm is getting closer to the parameter regime of the
Earth’s core with both a lower Ekman number and a lower magnetic
Prandtl number (E = 10⁻⁶, Pm = 0.2).

The values of the parameters introduced in Section 2.2 can be
estimated to a certain extent for the Earth. Values for the mag-
netic Reynolds number and the Elsasser number can be estimated
to Rm = 400 – 4000 and Λ = 5 – 50 (Aubert et al. 2017). Sim-
ilarly the ratio of magnetic to kinetic energy is likely in the range
M = 10³ – 10⁴ (Tassin et al. 2021). For the dipolar fraction of the
magnetic field, we will here take the range f<sup>†</sup><sub>dip</sub> = 0.35 – 0.75 sug-
gested by Davies et al. (2022). Table 1 compares the time-averaged
values of the dipolar fraction, the energy ratio, the Elsasser num-
ber and the magnetic Reynolds number of the reference cases with
the Earth estimates. All the reference cases have energy ratios that
are far too small compared to the expectation for the Earth. Cases
E1e-5_SF1 and E1e-5_SF2 nevertheless have M ≫ 10, the overline
standing for time average. For this reason, these two cases are called
‘SF’ for ‘Strong Field’ (Schwaiger et al. 2021). The magnetic over
kinetic energy ratio is of the order of 1 or larger for the three other
cases. However, these three cases have magnetic Reynolds numbers
that match the values expected for the Earth, while the two strong
field dynamos have values of Rm that are too low or just at the limit
of the Earth’s range. We thus call these three cases ‘hRm’ for ‘high
Rm’. Cases E1e-4_hRm, E1e-5_hRm and E1e-5_SF1 have dipolar
fractions within the Earth’s range, while E1e-5_SF2 and E1e-6_hRm
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Ekman number *E* and magnetic Prandtl number *Pm*, for the five reference dynamos E1e-4_hRm, E1e-5_hRm, E1e-6_hRm, E1e-5_SF1 and E1e-5_SF2 compared to previous studies focusing on the effect of large scale heterogeneous heat flux (Olson & Christensen 2002; Olson & Amit 2014; Sahoo *et al.* 2016). All the dynamo models are obtained for *Pr* = 1. The Rayleigh number of the dynamo models are 1.5 × 10⁸, 5 × 10⁹, 9.6 × 10¹¹, 5 × 10⁸ and 1.5 × 10⁸ for E1e-4_hRm, E1e-5_hRm, E1e-6_hRm, E1e-5_SF1 and E1e-5_SF2, respectively. The expected values for the Earth are *E* ≃ 10⁻¹⁵ and *Pm* ≃ 10⁻⁶.</figcaption></figure>

<table><caption>**Table 1.** Characteristics of the reference dynamos compared to estimated values for the Earth.</caption><thead><tr><th>Name</th><th>f&#773;<sub>dip</sub></th><th>f&#773;<sup>&dagger;</sup><sub>dip</sub></th><th>M&#773;</th><th>&Lambda;&#773;</th><th>Rm&#773;</th></tr></thead><tbody><tr><td>E1e-4_hRm</td><td>0.45</td><td>0.58</td><td>0.95</td><td>13.1</td><td>649</td></tr><tr><td>E1e-5_hRm</td><td>0.38</td><td>0.53</td><td>3.49</td><td>16.8</td><td>987</td></tr><tr><td>E1e-6_hRm</td><td>0.72</td><td>0.80</td><td>3.53</td><td>5.8</td><td>575</td></tr><tr><td>E1e-5_SF1</td><td>0.61</td><td>0.70</td><td>32.12</td><td>16.9</td><td>399</td></tr><tr><td>E1e-5_SF2</td><td>0.83</td><td>0.79</td><td>53.50</td><td>8.13</td><td>214</td></tr><tr><td>Earth</td><td>-</td><td>0.35-0.75<sup>a</sup></td><td>10<sup>3</sup>-10<sup>4</sup><sup>b</sup></td><td>5-50<sup>c</sup></td><td>400-4000<sup>c</sup></td></tr></tbody></table>

have larger dipolar fractions. All the cases have Earth-like Elsasser numbers. It can be noted that the strong field cases have low magnetic Reynolds numbers, while the high *Rm* dynamo cases have more comparable magnetic and kinetic energies with *M* ≳ 1. The inability to meet both criteria for a single dynamo model is a consequence of the distance between the parameter regime of dynamo models and the parameter regime of the Earth's core (Aubert *et al.* 2017; Nataf & Schaeffer 2024). In this study, we consider dynamo models that satisfy one criterion or the other in order to examine whether these two types of dynamo models behave differently in the presence of heat flux heterogeneities. This choice enables us to study the effect of heat flux heterogeneities on models close to the dipolar multipolar transition (high *Rm* models) and on models away from the transition (strong field models).

### 2.3.2 Heat flux patterns

We study the effect of heterogeneous heat flux patterns at the top of the core applied to the reference dynamo cases described above. We

focus on large scale heterogeneity patterns with spherical harmonics of degree 1 or 2. The patterns are shown in Fig. 2. Note that, in this figure, areas shown in red and blue correspond to a positive and negative heat flux, respectively, that is, an increased and decreased heat flux from the core to the mantle, respectively. Pattern $Y_{1,0}$ corresponds to a positive heat flux in the northern hemisphere and a negative heat flux in the southern hemisphere. Pattern $Y_{1,1}$ is the same as $Y_{1,0}$, but tilted by 90° in latitude. $+Y_{2,0}$ imposes a positive heat flux at the poles and a negative heat flux at the equator. $-Y_{2,0}$ is the opposite of $+Y_{2,0}$, imposing a negative heat flux at the poles and a positive heat flux at the equator. Finally, $Y_{2,2}$ imposes two antipodal heat flux highs at the equator, surrounded by two antipodal heat flux lows. Negative amplitudes for the $Y_{1,1}$ and $Y_{2,2}$ patterns only implies a rotation in longitude of the patterns that is irrelevant for the dynamo. Neither does a sign change of the $Y_{1,0}$ pattern affect the behaviour of the dynamo due to the symmetry of the dynamo equations. These heat flux patterns directly affect the flow through the generation of temperature gradients at the top of the core. The $Y_{1,1}$ and $Y_{2,2}$ patterns have been shown to generate an equatorial circulation with localized downwellings and upwellings offset from

<footnote><sup>a</sup>Davies *et al.* (2022); <sup>b</sup>Tassin *et al.* (2021); <sup>c</sup>Aubert *et al.* (2017).</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2. Heterogeneous heat flux patterns $p(\theta, \phi)$ used in the dynamo simulations. The positive and negative heat flux areas are shown in red and blue, respectively. The dashed lines show the positions where the patterns are equal to zero. The patterns are shown in a Mollweide projection.</figcaption></figure>

the high and low heat flux patches (Zhang & Gubbins 1992; Dietrich et al. 2016), which is verified in our simulations. The $Y_{1,0}$ and $Y_{2,0}$ patterns generate thermal winds that are equatorially antisymmetric for $Y_{1,0}$ and dominantly westward/eastward for $+Y_{2,0}$ and $-Y_{2,0}$, respectively (Amit et al. 2011; Dietrich et al. 2017). The imposed heat flux at the CMB is given by the product of an amplitude $\delta q_o$ and a pattern $p$,

$$q_o(\theta, \phi) = \delta q_o p(\theta, \phi). \qquad (25)$$

The normalization of the pattern $p$ ensures that $2\delta q_o = \max(q_o) - \min(q_o)$. We choose to scale the heterogeneity amplitude using the buoyancy flux at the inner boundary by defining a dimensionless amplitude $\delta q_o^* = \delta q_o/q_i$. Our scaling differs from most previous studies, where heat flux heterogeneities have been more commonly scaled by the mean outer boundary flux (Olson et al. 2010; Olson & Amit 2014; Mound & Davies 2023; Dannberg et al. 2024; Terra-Nova & Amit 2024), but it has the advantage of being well-defined even when $F_o = 0$. It is also relevant for the Earth, for which the main forcing is thought to occur at the inner core boundary in the form of a compositional flux (Loper 1978), while the radial temperature profile could be close to isentropic ($F_o \simeq 0$). The same scaling was used by Aubert et al. (2013) and Nakagawa & Davies (2022). The heterogeneity amplitude normalized this way can be expressed for the Earth as

$$\delta q_o^* = \frac{\alpha \delta q_o^{th}}{c_p \dot{r}_i \left( \Delta \rho + \frac{\alpha \rho_c L}{c_p} \right)}, \qquad (26)$$

where $2\delta q_o^{th}$ is the peak-to-peak amplitude of the CMB heat flux, $c_p$ is the heat capacity of the liquid iron, $\dot{r}_i$ is the growth rate of the inner core, $\Delta \rho$ is the compositional density jump at the inner core boundary, $\rho_c$ is the mean density of the core and $L$ is the latent heat of liquid iron crystallization. Table 2 shows estimates of the maximum values of $\delta q_o^*$ using mantle convection simulations corrected for true polar wander by Frasson et al. (2024). The amplitude of the total heat flux pattern heterogeneities reach up to $\delta q_o^* = 10.3$ per cent. The heterogeneities for a given spherical harmonic coefficient are lower, but can be as high as $\delta q_o^* = 3.4$ per cent for $Y_{2,2}$. Note that we impose $F_o = 0$ and use the Boussinesq approximation in our dynamo models, so that the average heat flux extracted by the mantle is identified with the isentropic heat flux and the flux of light elements between the core and the mantle is null. The heterogeneous flux will thus be stabilizing in areas where $q_o(\theta, \phi) < 0$ (i.e. lower than the isentropic flux).

All five patterns have been applied to the reference case E1e-4_hRm. The four other cases require much more computing resources due to their lower Ekman numbers. We thus focus on the $\pm Y_{2,0}$ patterns for cases E1e-5_hRm, E1e-6_hRm, E1e-5_SF1 and E1e-5_SF2.

# 3 RESULTS

## 3.1 Reference homogeneous dynamos

We shortly describe here the reference dynamo simulations with homogeneous heat flux. The five reference dynamos cover a wide range of parameters, with Ekman numbers between $10^{-4}$ and $10^{-6}$ and magnetic Prandtl numbers between 3 and 0.2 (Fig. 1). The Rayleigh number is also significantly varied among the reference cases.

In Fig. 3 are shown the time evolution of the dipolar fraction $f_{dip}$, the ratio of magnetic to kinetic energy $M$ and the dipole latitude (angle between the equator and the magnetic dipole) $\theta_{dip}$. All the reference cases have strong dipoles with $f_{dip} \ge 0.39$. The E1e-4_hRm reference case displays magnetic reversals. The ratio of magnetic to kinetic energy is around unity. This reference dynamo runs relatively quickly, and is thus useful for extensive explorations. Case E1e-5_hRm has a similar reversing behaviour to case E1e-4_hRm. Lowering the Ekman number allows for a higher magnetic to kinetic energy ratio ($\bar{M} = 3.5$). Dynamo E1e-5_SF1 has the same Ekman number as E1e-5_hRm but a ten times smaller Rayleigh number. The magnetic Prandtl number is also slightly higher ($Pm = 2$ for E1e-5_hRm, $Pm = 3$ for E1e-5_SF1). Lowering the Rayleigh number and increasing the magnetic Prandtl number result in a magnetic to kinetic energy ratio larger than 20. The dipolar fraction also significantly increases and the dipole becomes more stable. No reversals are observed in this case for the duration of the simulation. The E1e-5_SF2 is obtained by decreasing again the Rayleigh number by 30 per cent. This case has the largest magnetic to kinetic energy ratio ($\bar{M} = 53$). It is also the reference case with the highest dipolar fraction. No reversals are observed in this case. Finally E1e-6_hRm is obtained using $E = 10^{-6}$. The magnetic energy is larger than the kinetic energy, with a similar value on average to case E1e-5_hRm. No reversals are observed.

Though all the reference cases have a dipole-dominated magnetic field, they lie at a different distance from the dipolar-multipolar transition. Tassin et al. (2021) found a transition from strong dipoles towards multipolar magnetic field when the ratio of magnetic to kinetic energy becomes lower than about 1. Following this argument, the three high Rm cases (E1e-4_hRm, E1e-5_hRm and E1e-6_hRm) should be close to the transition with a magnetic energy similar or slightly larger than the kinetic energy. Cases E1e-4_hRm and E1e-5_hRm are reversing, suggesting indeed that they lie close to the transition. Case E1e-6_hRm does not reverse polarity for the simulated time, though the eventual occurrence of reversals cannot be excluded due to the limited duration of the simulation. Cases E1e-5_SF1 and E1e-5_SF2 have magnetic energies significantly higher than the kinetic energies, and can thus be expected to be far from the dipolar-multipolar transition according to the energy ratio criterion.
<table><caption>**Table 2.** Maximal heterogeneity amplitude of CMB heat flux in the mantle convection simulations corrected for true polar wander by Frasson *et al.* (2024). The values below label $\delta q_o^*$ give the maximum peak-to-peak amplitude in the total heat flux patterns for the different cases. The values below $\delta q_o^*(Y_{1,m})_{\text{max}}$ give the maximum peak-to-peak amplitude of the degree *l* and order *m* of the CMB heat flux time-series in the different cases. For the $Y_{2,0}$ pattern, the first value gives the minimum negative amplitude (corresponding to $-Y_{2,0}$) and the second gives the maximum positive amplitude (corresponding to $+Y_{2,0}$). In each column, we emphasize with bold font the highest value reached among the mantle convection cases (all defined in Frasson *et al.* 2024).</caption><thead><tr><th>Case</th><th>δq<sub>o</sub><sup>*</sup><sub>max</sub><br>(per cent)</th><th>δq<sub>o</sub><sup>*</sup>(Y<sub>1,0</sub>)<sub>max</sub><br>(per cent)</th><th>δq<sub>o</sub><sup>*</sup>(Y<sub>1,1</sub>)<sub>max</sub><br>(per cent)</th><th>δq<sub>o</sub><sup>*</sup>(-Y<sub>2,0</sub>)<sub>max</sub><br>(per cent)</th><th>δq<sub>o</sub><sup>*</sup>(+Y<sub>2,0</sub>)<sub>max</sub><br>(per cent)</th><th>δq<sub>o</sub><sup>*</sup>(Y<sub>2,2</sub>)<sub>max</sub><br>(per cent)</th></tr></thead><tbody><tr><td>MF1</td><td><b>10.3</b></td><td>1.6</td><td>2.4</td><td><b>2.9</b></td><td>2.6</td><td>3.2</td></tr><tr><td>MF2</td><td><b>10.3</b></td><td><b>1.7</b></td><td>2.3</td><td>0</td><td><b>2.7</b></td><td>3.1</td></tr><tr><td>MF*</td><td><b>10.3</b></td><td>1.2</td><td><b>2.6</b></td><td>2.4</td><td>1.8</td><td><b>3.4</b></td></tr><tr><td>MC1</td><td>7.1</td><td>0.9</td><td>1.3</td><td>1.3</td><td>0</td><td>1.0</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Time-series of the dipolar fraction ($f_{dip}$), the ratio of magnetic to kinetic energy ($M$) and of the dipole latitude for the five reference geodynamo simulations considered in this study. Note that the horizontal scale for time is different for each case. The two first dipole diffusion time of the simulations have been removed for E1e-6_hRm, E1e-5_SF1 and E1e-5_SF2. The simulations for E1e-4_hRm and E1e-5_hRm have been started from an already converged state.</figcaption></figure>

## 3.2 Categories of dynamo behaviours

Before looking in details into the effect of heat flux heterogeneities, we will first define three distinct dynamical behaviours: (1) ‘Dipolar stable’, if $\overline{f}_{dip} > 0.25$ (or $\overline{f}_{dip}^\dagger > 0.35$) and no reversals or excursions are observed during the simulation; (2) ‘Reversing’, if $\overline{f}_{dip} > 0.25$ (or $\overline{f}_{dip}^\dagger > 0.35$) and reversals or excursions are observed during the simulation; (3) ‘Multipolar’, if $\overline{f}_{dip} < 0.25$ (or $\overline{f}_{dip}^\dagger < 0.35$). The dipolar stable behaviour corresponds to dynamos that are far from a dipolar–multipolar transition. The reversing behaviour corresponds to dynamos that are still dipolar but closer to a transition. In some cases, the dynamo spends a significant amount of time in both the dipolar and the multipolar regime. Using the aforementioned criteria, these bistable dynamos would be categorized as reversing due to their relatively high dipolar fraction on average. A bistable behaviour is however different from the reversing behaviour expected for the Earth, for which we expect the dipolar regime to dominate. We thus define a fourth behaviour: (4) ‘Bistable’, if $\overline{f}_{dip} > 0.25$ and if $\overline{f}_{dip} < 0.25$ for at least 20 per cent of the simulation time.

The separation between the dipolar and multipolar behaviour based on the time-averaged value of $f_{dip}$ is similar to what is done in previous studies (Christensen & Aubert 2006; Soderlund *et al.* 2012; Tassin *et al.* 2021). A separation between non-reversing and

reversing dynamos based on the relative time spent by the dipole below a certain latitude has been previously used (Wicht & Meduri 2016; Meduri *et al.* 2021). Though this kind of separation could work for our simulations, we chose to use our simpler classification that relies only on the statistics of $f_{dip}$. We moreover note that, for bistable dynamos, the bistability appears most clearly in the distribution of $f_{dip}$ while it is much less clear in the distribution of the dipole tilt.

Depending on the geological period, the behaviour of the Earth’s dynamo would be dipolar stable or reversing. With these definitions, the separation between a dipolar and a multipolar behaviour depends directly on the dipolar fraction at the CMB. Although dipolar stable dynamos and reversing dynamos are not discriminated through their dipolar fraction, we can expect the variations of $f_{dip}$ to be higher in the reversing dynamo models, as larger variations of the dipole moment are required for reversals to occur (Buffett *et al.* 2013). We quantify this variation by defining the relative standard variations of the dipolar fraction as $\tilde{\sigma}(f_{dip}) = \frac{\sigma(f_{dip})}{f_{dip}}$ where $\sigma(f_{dip})$ is the standard variation of the dipolar fraction. In Fig. 4 are shown the values of $\overline{f}_{dip}$ and $\tilde{\sigma}(f_{dip})$ for all the simulations presented in this study, including simulations with heterogeneous heat fluxes. We indeed find a separation between dipolar stable dynamos and reversing dynamos based on the relative fluctuations of the dipolar fraction. The
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Fluctuation $\tilde{\sigma}(f_{\text{dip}})$ and mean value $\overline{f_{\text{dip}}}$ of the dipolar fraction at the CMB in all the performed simulations (including simulations with heterogeneous heat fluxes). The marker style shows the reference case from which each dynamo is derived. The colour shows the dynamo behaviour, as defined in Section 3.2. All multipolar dynamos plot left of the vertical dashed line at $\overline{f_{\text{dip}}} = 0.25$, while the horizontal dashed line at $\tilde{\sigma}(f_{\text{dip}}) = 0.16$ separates the dipolar stable and reversing dynamos.</figcaption></figure>

dipolar stable dynamos all have $\tilde{\sigma}(f_{\text{dip}}) < 0.16$, while all reversing dynamos (and multipolar dynamos) have $\tilde{\sigma}(f_{\text{dip}}) > 0.16$.

Using these definitions we classify our reference cases E1e-4_hRm and E1e-5_hRm as reversing dynamos. The three other cases are classified as dipolar stable dynamos. However, we stress that observing no occurrence of a rare event such as reversals in a finite duration simulation cannot guarantee that it would never happen.

## 3.3 Effect of heat flux heterogeneities for the E1e-4_hRm case

In this section, we will explore how the large-scale heat flux patterns presented in Fig. 2 affect the behaviour of the dynamo case E1e-4_hRm. As described in Section 3.1, this reference case lies close to the dipolar–multipolar transition, showing a reversing behaviour. This reference case can thus be expected to be very sensitive to changes in the flux boundary conditions.

Violin plots in Fig. 5 display the probability distribution of the dipolar fraction, of the ratio of magnetic to kinetic energy, and of the dipole latitude as a function of the pattern amplitude for the $Y_{1,0}$, $Y_{1,1}$ and $Y_{2,2}$ patterns. The probability distributions for the $\pm Y_{2,0}$ patterns are shown in Fig. 6. The red shaded area shows the expected amplitudes obtained from mantle convection models given in Table 2. The grey shaded area shows the expected amplitudes for the total heat flux pattern also given in Table 2. The behaviours described in Section 3.2 are well illustrated by the probability distributions. The dipolar stable dynamos have the highest values of $\overline{M}$ and have low fluctuations of both $f_{\text{dip}}$ and $M$. They are obtained for moderate amplitudes of the $Y_{1,1}$ pattern, or using the $+Y_{2,0}$ pattern. The dynamo that has the lowest amplitude of the $Y_{2,2}$ pattern also shows a dipolar stable behaviour. The fluctuations of the dipolar fraction

in this case are however large, making it very likely for reversals or excursions to eventually occur. The multipolar dynamos show a dipolar fraction that fluctuates around low values, though the dipolar fraction can occasionally reach high values similar to the dipolar dynamos in some cases. The distribution of the dipole latitude covers all the range between -90° and +90°. They are obtained using the $Y_{1,0}$ pattern or the $-Y_{2,0}$ pattern. A multipolar solution is also obtained using a high amplitude of the $Y_{2,2}$ pattern that lies within the expected range for a total heat flux pattern, but that is higher than the expectations for a $Y_{2,2}$ pattern only. The reversing dynamos have similar probability distributions of $f_{\text{dip}}$ and $M$ to the dipolar stable dynamos, but with larger tails towards low values. The distribution of the dipole latitude shows two maxima near $\pm 90°$. In some cases, the dynamo spends a significant amount of time in both the dipolar ($f_{\text{dip}} > 0.25$) and the multipolar ($f_{\text{dip}} < 0.25$) state. These bistable dynamos are found with the $Y_{2,2}$ pattern for intermediate amplitudes within the estimated range and with the $Y_{1,1}$ pattern for an amplitude larger than the expected range. The bistable nature of these dynamos is well visible in the cases obtained using the $Y_{2,2}$ pattern. In both bistable cases, the probability distribution of the dipolar fraction and the energy ratio is clearly bimodal. The switching between a dipolar state and a multipolar state is visible in the time-series of the dipolar fraction as shown in Fig. A1. The bistable dynamo obtained using the $Y_{1,1}$ pattern shows probability distributions that are less clearly bimodal. Nevertheless, this case spent a significant amount of time in a multipolar state because of a relatively low dipolar fraction on average and very large fluctuations. A reversing dynamo is obtained for a weak amplitude of the $Y_{1,0}$ pattern, with a behaviour similar to the reference case, before the dynamo transition towards a multipolar behaviour for higher amplitudes.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Probability distributions of the dipolar fraction at the CMB ($f_{\text{dip}}$), of the ratio of magnetic to kinetic energy ($M$), and of the dipole latitude ($\theta_{\text{dip}}$) for the geodynamo simulations derived from the reference case $\text{E1e-4\_hRm}$. The horizontal width of the violins is a measure of the likelihood. The probability distributions are computed using 100 estimation points, each with a relative bandwidth of 0.1. The distributions are given as a function of the amplitude of the heat flux patterns normalized by the inner boundary flux. The colours encode the dynamos category with the same colour code as in Fig. 4. The heat flux patterns are (a) $Y_{1,0}$, (b) $Y_{1,1}$ and (c) $Y_{2,2}$. The expected range of amplitude for the Earth are shown in red for each pattern and in grey for a full pattern (see Table 2).</figcaption></figure>

The opposite effects of the $+Y_{2,0}$ pattern and the $-Y_{2,0}$ pattern is particularly striking. As shown in Fig. 6, the $+Y_{2,0}$ pattern preserves the dipole-dominated structure of the magnetic field, even for large amplitudes, and slightly increases the magnetic energy relatively to the kinetic energy. In contrast, the $-Y_{2,0}$ significantly decreases the energy ratio and triggers a transition towards a multipolar behaviour even for the lowest heterogeneity amplitudes. The increase of the magnetic field intensity or the flow intensity when the $+Y_{2,0}$ or $-Y_{2,0}$ patterns are applied can be seen in more details through the variations of the Elsasser and magnetic Reynolds numbers shown in Fig. A2.

## 3.4 Systematic study on the effect of equatorial cooling and polar cooling

The $+Y_{2,0}$ pattern imposes a negative heat flux at low latitudes, and a positive heat flux at higher latitudes. This pattern thus favours a polar cooling of the core, while the equatorial region tends to be stratified. Conversely, the $-Y_{2,0}$ pattern favours an equatorial cooling of the core and acts against convection in higher latitudes. In this section, we will study in more details the effect of both situations on a wide range of dynamo models.

### 3.4.1 Effect on high Rm dynamos

The probability distributions of $f_{\text{dip}}$ and $M$ are shown as a function of the amplitude of the $Y_{2,0}$ pattern in Fig. 7. The geodynamo models $\text{E1e-5\_hRm}$ (Fig. 7b) and $\text{E1e-6\_hRm}$ (Fig. 7c) are affected by the $Y_{2,0}$ pattern in a very similar way to model $\text{E1e-4\_hRm}$ (Fig. 7a). In those three high $Rm$ models, the dynamo becomes multipolar when

$a-Y_{2,0}$ heat flux pattern is applied, with a decrease of the magnetic to kinetic energy ratio. For these three cases, the transition occurs for amplitudes of the $-Y_{2,0}$ pattern within the expected range for the Earth. The $+Y_{2,0}$ pattern preserves a strong stable dipole and increases the magnetic energy, even for very large amplitudes. It is noteworthy that very similar effects on the dynamo behaviour are observed for three geodynamo models with Ekman numbers that vary by two orders of magnitude.

### 3.4.2 Effect on strong field dynamos

The effect of the $Y_{2,0}$ pattern is very different for the two strong field cases $\text{E1e-5\_SF1}$ (Fig. 7d) and $\text{E1e-5\_SF2}$ (Fig. 7e). The amplitude of the magnetic energy relative to the kinetic energy decreases when the amplitude of the pattern is increased, regardless of its sign. The decrease of the magnetic energy is however faster with a $-Y_{2,0}$ pattern. A more prominent difference compared to the high $Rm$ cases exists in the response of the dipolar fraction to the heat flux pattern. For moderate amplitudes ($|\delta q_o^*| \le 5$ per cent) the dipolar fraction is higher for a $-Y_{2,0}$ pattern than for a $+Y_{2,0}$ pattern, opposite to the three high $Rm$ cases. For larger amplitudes, the fluctuations of $f_{\text{dip}}$ and $M$ significantly increase with the $-Y_{2,0}$ pattern, while the fluctuations do not significantly vary with the $+Y_{2,0}$ pattern. In both the $\text{E1e-5\_SF1}$ case and the $\text{E1e-5\_SF2}$ case, the dynamo behaviour is not significantly affected by the $Y_{2,0}$ pattern for amplitudes within the expected range for the Earth. The dynamo becomes bistable in the $\text{E1e-5\_SF1}$ case and reversing in the $\text{E1e-5\_SF2}$ case for amplitudes of the $-Y_{2,0}$ pattern larger than the expected values for this pattern alone but within the expected range for the total heat flux. The dynamo becomes multipolar in the $\text{E1e-5\_SF1}$ case for larger
<figure><img src="image_8.png" /><figcaption>**Figure 6.** As in Fig. 5 for $Y_{2,0}$. Positive amplitudes $\delta q_o^*$ correspond to a polar cooling, while negative amplitudes correspond to an equatorial cooling.</figcaption></figure>

amplitudes of the $-Y_{2,0}$ pattern that are probably unrealistic for the Earth. No multipolar dynamos are obtained in the E1e-5_SF2 case, even for unrealistically large amplitudes. In this latter case, the dynamo stays reversing when higher amplitudes of the $-Y_{2,0}$ pattern are used. Interestingly, the transition induced by the $-Y_{2,0}$ pattern from a strongly dipolar dynamo to a multipolar dynamo in E1e-5_SF1 occurs while the magnetic to kinetic energy ratio stays significantly larger than unity. For example, for the multipolar dynamo obtained with $\delta q_o^* = -14.5$ per cent, we have $\bar{M} = 6$ and $M > 2$ for the whole simulation. The value of $\bar{M}$ stays also significantly larger than 1 for the reversing dynamos obtained with the E1e-5_SF2 model, though the minimum value falls closer to unity in these cases. Such a behaviour is very different from the one observed for the high $Rm$ models, for which the transition to a multipolar state is associated with a decrease of $\bar{M}$ to values that are of order 1 or below.

The eight reversing dynamos in the E1e-5_SF2 case have fluctuations of the dipolar fraction significantly larger than the reversing cases obtained with the high $Rm$ models. In this respect, they are more similar to the bistable dynamos obtained in the E1e-5_SF1 case and in the E1e-4_hRm case using the $Y_{1,1}$ pattern (see Fig. 5b). The dipolar fraction in these cases can range from almost 0 to 0.9. Similarly, the ratio of magnetic to kinetic energy fluctuates over more than one order of magnitude. They are nevertheless classified as reversing using our criteria as they have a large dipolar fraction on average and they spend a limited amount of time in a multipolar state. In details, the shape of the probability distributions vary between these reversing dynamos. Some of these reversing dynamos spend most of their time in a strongly dipolar state, in which the magnetic energy largely dominates the kinetic energy

[see $\delta q_o^* = -12.2$ per cent in Fig. 7(e)]. In other cases, the dipolar fraction decreases more frequently and the magnetic energy, while always dominating the kinetic energy, exhibits larger fluctuations [see $\delta q_o^* = -11.0$ per cent in Fig. 7(e)].

## 3.5 Relation between mean zonal flows and the dipole stability

### 3.5.1 Equatorially antisymmetric heat flux ($Y_{1,0}$)

In Fig. 8 we plot the meridional cuts of the azimuthal component of the magnetic field and of the velocity, averaged both in time and in the azimuthal direction, for the simulations using the E1e-4_hRm model and the $Y_{1,0}$ and $Y_{2,2}$ patterns. The field lines of the poloidal magnetic field and the streamlines of the meridional velocity field averaged in time and in the azimuthal direction are also displayed. In the homogeneous simulation (Fig. 8a), the zonal flow is strong within the cylinder tangent to the inner core (tangent cylinder) and weaker outside. The zonal flow is highly symmetric with respect to the equator and the magnetic field lines show a clear dipolar geometry. The $Y_{1,0}$ pattern breaks the equatorial symmetry in the zonal flow even for the lowest amplitude through the imposed thermal winds (Fig. 8b, see also Amit et al. 2011). The antisymmetry becomes even clearer for larger amplitudes, and the magnetic field becomes multipolar. The $Y_{2,2}$ pattern does not affect significantly the flow for moderate amplitudes, and the magnetic field stays dipolar (Figs 8 f–h). However, for the largest amplitude, the magnetic field becomes multipolar and the zonal flows become strongly antisymmetric (Fig. 8i). Interestingly, zonal flows adopt a geometry
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Probability distributions of the dipolar fraction at the CMB ($f_{\text{dip}}$) and of the ratio of magnetic to kinetic energy ($M$) for the geodynamo models using the $Y_{2,0}$ pattern. The reference dynamo models are: (a) E1e-4_hRm, (b) E1e-5_hRm, (c) E1e-6_hRm, (d) E1e-5_SF1, (e) E1e-5_SF2. The expected range of amplitude for the Earth are shown in red for each pattern and in grey for the full pattern (see Table 2).</figcaption></figure>

very similar to the ones forced by the $Y_{1,0}$ pattern, with an eastward flow in the north of the tangent cylinder and a westward flow in the south.

### 3.5.2 Equatorial and polar cooling of the core

We similarly show in Fig. 9 meridional cuts of the azimuthal velocity and magnetic field for cases E1e-6_hRm and E1e-5_SF1 using the $Y_{2,0}$ heat flux pattern. Case E1e-6_hRm is representative of the three high $Rm$ models, while case E1e-5_SF1 is representative of the two strong field models. The same meridional cuts for the other dynamo models are shown in Fig. A3. Without heat flux heterogeneities, the zonal flow in E1e-6_hRm (Fig. 9b) is similar to the one in E1e-4_hRm (Fig. 8a). The $+Y_{2,0}$ pattern increases the strength of the westward azimuthal flow outside the tangent cylinder, which becomes of a similar amplitude to the azimuthal flow inside the tangent cylinder (Fig. 9c). This results in a larger azimuthal component of the magnetic field outside the tangent cylinder, which can be explained by an increased $\omega$-effect. For an equatorial cooling, the dynamo becomes multipolar and the magnetic field is thus dominated by smaller scales (Fig. 9a). The equatorial symmetry of zonal flows is broken outside the tangent cylinder, though the antisymmetry remains lower than what is obtained with the $Y_{1,0}$ pattern and

varies depending on the model (see Fig. A3). The flow is dominated by a strong eastward flow in the northern hemisphere and a more localized westward flow in the southern hemisphere. The reinforcement of the westward flow for a polar cooling and the occurrence of a mostly eastward flow for an equatorial cooling is consistent with thermal winds driven by the heat flux heterogeneities (Dietrich et al. 2017).

The E1e-5_SF1 case without heat flux heterogeneities (Fig. 9e) shows a dipolar magnetic field, similar to the E1e-6_hRm case. The westward azimuthal flow outside the tangent cylinder is however much stronger in the E1e-5_SF1 case. As in the E1e-6_hRm case, the $+Y_{2,0}$ pattern strengthens the westward azimuthal flow and the azimuthal component of the magnetic field, while the poloidal magnetic field keeps a dipolar geometry (Fig. 9f). As in the E1e-6_hRm case again, the $-Y_{2,0}$ pattern imposes an eastward azimuthal flow outside the tangent cylinder and the magnetic field becomes multipolar (Fig. 9d). Contrary to the high $Rm$ cases, the equatorial symmetry of the zonal flows is maintained. In both high $Rm$ models and strong field models a polar cooling of the core is found to decrease the dipolar fraction of the magnetic field, despite maintaining a non-reversing dipole.

All the reference dynamos considered in this study have a dominantly westward zonal flow outside the tangent cylinder when averaged in time. The high $Rm$ and the strong field dynamos differ in
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Meridional cuts showing the azimuthal component of the magnetic field and of the velocity averaged in time and in the azimuthal direction relative to the root mean square of the fields for the dynamo cases using the $Y_{1,0}$ and the $Y_{2,2}$ patterns. (a) E1e-4_hRm reference case without heat flux heterogeneities. (b–e) Dynamo simulations using the $Y_{1,0}$ pattern with an increasing amplitude. (f–i) Dynamo simulations using the $Y_{2,2}$ pattern with an increasing amplitude. For each figure, the azimuthal magnetic field is shown on the left and the azimuthal velocity field is on the right. The grey contours show the poloidal magnetic field lines and the streamlines of the meridional circulation. The flow is clockwise where the field lines are plain and anticlockwise where the field lines are dashed.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Meridional cuts showing the azimuthal component of the magnetic field and of the velocity averaged in time and in the azimuthal direction relative to the root mean square of the fields for several dynamo cases using a $Y_{2,0}$ heat flux pattern. The time-averaged dipolar fraction of the magnetic field is given in each case. Figs a, b and c correspond to the E1e-6_hRm case. Figs d, e and f correspond to the E1e-5_SF1 case. Figs a and d are obtained for a negative $Y_{2,0}$ heat flux pattern (equatorial cooling). Figs c and f are obtained for a positive $Y_{2,0}$ pattern (polar cooling). Figs b and e are obtained without heat flux heterogeneities. For each figure, the magnetic field is shown on the left and the velocity field is on the right. The grey contours show the field lines of the poloidal magnetic field and the contours of the meridional circulation. The flow is clockwise where the field lines are plain and anticlockwise where the field lines are dashed. Note that the inner core is insulating for E1e-6_hRm and conducting for E1e-5_SF1. The colour scale for the magnetic field in the E1e-6_hRm case with the $-Y_{2,0}$ pattern is between ±10 per cent. All the other colour scales are between ±100 per cent. The heterogeneity amplitudes are: (a) $\delta q_o^* = -9.1$ per cent, (b) $\delta q_o^* = 0$ per cent, (c) $\delta q_o^* = +9.1$ per cent, (d) $\delta q_o^* = -14.5$ per cent and (e) $\delta q_o^* = +14.5$ per cent.</figcaption></figure>

the strength of this westward zonal flow, the strong field dynamos showing a stronger zonal flow. In order to study the variations of this zonal flow we introduce the mean angular velocity of the fluid relative to the RMS velocity defined as

$$ \Omega_f(t) = \frac{L_z(t) r_o}{I_0 U_{\text{RMS}}}, \qquad (27) $$

where $L_z$ is the angular momentum of the liquid core around the spin axis, and $I_0$ is the axial moment of inertia of a homogeneous spherical shell. We chose $\Omega_f$ as a simple proxy for the zonal flows outside the tangent cylinder, as the angular momentum $L_z$ is dominated by flows far from the rotation axis. The dependence of $\Omega_f$ on the $Y_{2,0}$ pattern amplitude is shown in Fig. 10. The average time $T_{\text{rev}}$ spent by the dipole in a given hemisphere (defined in eq. 24) is
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Mean angular velocity of the flow $\bar{\Omega}_f$ as a function of the $Y_{2,0}$ heat flux pattern amplitude. The colour scale shows the average time spent by the magnetic dipole in a single hemisphere (i.e. the average duration of magnetic chrons). The dipolar stable dynamos that do not reverse polarity are shown in blue. The multipolar dynamos showing a weak dipole are shown in orange. Only multipolar dynamos plot above the orange horizontal dashed line at $\bar{\Omega}_f = 0.13$, while only dipolar stable dynamos plot below the blue horizontal dashed line at $\bar{\Omega}_f = -0.13$. The grey shaded area highlights the transition between these two regimes, in which dipolar and multipolar dynamos are observed.</figcaption></figure>

shown with the colour scale. The dipolar stable dynamos that do not reverse are shown in blue, and the multipolar dynamos are shown in orange. Note that not all the dynamo simulations are shown as the angular momentum has not been recorded for all the simulations (see Appendix A). The difference between the amplitudes of the zonal flows between the high $Rm$ and the strong field reference cases is visible by looking at the value of $\bar{\Omega}_f$ for $\delta q_o^* = 0$. The mean angular velocity is only a few percents of the RMS velocity in the high $Rm$ cases. $\bar{\Omega}_f$ is negative for both E1e-4_hRm and E1e-5_hRm, corresponding to a dominantly westward azimuthal flow. The mean angular velocity is slightly positive for the E1e-6_hRm case despite the westward zonal flow at the equator shown in Fig. 9(b). The strengthening of the westward flows caused by a polar cooling of the core is clearly visible in the decrease of $\bar{\Omega}_f$ for all the models. Only non-reversing dipolar dynamos are obtained for strong westward flows ($\bar{\Omega}_f < -0.13$). Conversely, only multipolar dynamos are associated with strong eastward flows ($\bar{\Omega}_f > 0.13$). Weaker zonal flows ($-0.13 \le \bar{\Omega}_f \le 0.13$) are associated with both dipolar and multipolar dynamos.

## 3.6 Capturing the dipolar–multipolar transition

Our results show that dipolar dynamos can become multipolar when heterogeneous heat flux patterns are applied at the top of the core. The value of $\delta q_o^*$ for which the dynamo becomes multipolar depends on the heat flux pattern (see Figs 5 and 6) and the dynamo models (see Fig. 7). Several authors suggested that the transition was controlled by the local Rossby number, $Ro_l$ (defined in eq. 22), as they observed a transition from dipole-dominated dynamos towards multipolar dynamos when $Ro_l \gtrsim 0.12$ (Christensen & Aubert 2006; Olson & Christensen 2006; Christensen 2010; Wicht & Tilgner 2010). This criteria however fails to capture the transition for strong field dynamos (Menu et al. 2020; Tassin et al. 2021), a regime which

holds at least for our reference cases E1e-5_SF1 and E1e-5_SF2. Tassin et al. (2021) found that the transition can be better constrained by the ratio between the magnetic and kinetic energy and occurs for $\bar{M} \lesssim 1.1$. In Fig. 11 we show, for all our dynamo models, the dipolar fraction $f_{dip}$ as a function of the local Rossby number, the ratio of magnetic to kinetic energy and the magnetic Reynolds number. As expected, there is a tendency of a decreased dipolar fraction when $Ro_l$ or $Rm$ is increased and when $\bar{M}$ is decreased. We quantify this correlation by computing a linear fit between the dipolar fraction and the logarithm of the various criteria. The fit is shown by the grey dashed lines in Fig. 11. This fit shows a poor correlation between the dipolar fraction and the local Rossby number with a correlation coefficient $r^2 = 0.45$. The correlation is slightly improved with the ratio of magnetic to kinetic energy ($r^2 = 0.57$) and with the magnetic Reynolds number ($r^2 = 0.61$). None of these three criteria accurately capture the dipolar–multipolar transition. Both dipolar and multipolar dynamos are obtained for $0.04 < Ro_l < 0.1$ and for $0.6 < \bar{M} < 5.9$. This shared range is more limited but also exists for the magnetic Reynolds number, with both dipolar and multipolar dynamos found for $558 < Rm < 988$. Note that the transition is well captured by all the three criteria when a single reference dynamo case is considered. The transition nevertheless occurs for different critical values depending on the reference case. All the multipolar dynamos have values of $Ro_l$ that are lower than the threshold of $Ro_l \simeq 0.1$ initially suggested. The transition towards multipolar dynamos occurs at lower values of $Ro_l$ when the magnetic to kinetic energy ratio increases. Similarly, the transition towards multipolar dynamos occurs at higher values of $\bar{M}$ when $Ro_l$ decreases, yielding multipolar dynamos with $\bar{M}$ as high as 6. In Fig. 12, we introduce the parameter $M_{Za}^*$, which is defined as

$$M_{Za}^* = \frac{\Lambda_{Za}}{Rm_{Za}^2}, \qquad (28)$$

where $\Lambda_{Za}$ is the Elsasser number based on the zonal antisymmetric component of the magnetic energy and $Rm_{Za}$ is the magnetic Reynolds number based on the zonal antisymmetric component of the kinetic energy. Equivalently, it is also the ratio between the zonal antisymmetric part of the magnetic energy ($E_{bZa}$) and the zonal antisymmetric part of the kinetic energy ($E_{uZa}$) multiplied by the magnetic Ekman number $E_\eta = \frac{E}{Pm}$, such that $M_{Za}^* = \frac{E_{bZa}}{E_{uZa}} E_\eta$. In Fig. 12 we show $f_{dip}$ as a function of $\bar{M}_{Za}^*$ for our set of simulations. We find a better description of the evolution of the dipolar fraction using this parameter, with a correlation coefficient $r^2 = 0.79$. Only multipolar dynamos are found for $M_{Za}^* < 10^{-5}$, while only dipolar stable dynamos are found for $M_{Za}^* > 10^{-3}$. In between these two values, dipolar and multipolar dynamos are found. This is also where reversing dynamos are found.

# 4 DISCUSSION

## 4.1 Effect of heat flux patterns on the stability of the magnetic dipole

### 4.1.1 Destabilizing effect of the $Y_{1,0}$ and $-Y_{2,0}$ patterns

The magnetic dipole is affected in different ways by the five heat flux patterns considered in this study. The $Y_{1,0}$ and $-Y_{2,0}$ heat flux geometries very efficiently destabilize and weaken the magnetic dipole in the E1e-4_hRm case, while the $+Y_{2,0}$ pattern tends to stabilize the magnetic dipole. The effect obtained for the $Y_{2,0}$ pattern
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Dipolar fraction as a function of (a) the local Rossby number, (b) the magnetic to kinetic energy ratio and (c) the magnetic Reynolds number. The horizontal dashed line at a dipolar fraction of 0.25 separates the multipolar dynamos (below) from the dipolar dynamos (above). The transitional value of $Ro_l = 0.12$ initially suggested by Christensen & Aubert (2006) and of $M = 1.1$ suggested by Tassin *et al.* (2021) are drawn as the vertical red lines in figs a and b. The markers correspond to the dynamo model used and have the same meaning as in Fig. 4. The colours mark the amplitude of either $\bar{M}$ or $Ro_l$, as indicated at the colour bar.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Dipolar fraction of the magnetic field at the CMB as a function of $\bar{M}_{Za}^*$. The marker corresponds to the different reference dynamo models, and the colour gives the behaviour of the dynamo models (same as in Fig. 4). The orange and blue vertical dashed lines for $\bar{M}_{Za}^* = 10^{-5}$ and $\bar{M}_{Za}^* = 10^{-3}$ delimitate the range of $M_{Za}^*$ in which we obtain reversing dynamos.</figcaption></figure>

is consistent with previous results that showed a higher reversal frequency when the heat flux is higher at the equator than at the poles (Glatzmaier *et al.* 1999; Kutzner & Christensen 2004; Olson *et al.* 2010). Note however that Olson *et al.* (2010) obtained a higher reversal frequency with the $+Y_{2,0}$ pattern than with uniform heat flux conditions, while in our simulations the $+Y_{2,0}$ pattern seems to stabilize the magnetic dipole compared to uniform conditions. We notably confirm that the $-Y_{2,0}$ pattern, which imposes a high heat flux at the equator, is the pattern that affects the most the dynamo behaviour by weakening the magnetic field and destabilizing the magnetic dipole (i.e. it pushes towards a multipolar behaviour). The $Y_{1,0}$ pattern has also been consistently found to moderately destabilize and weaken the magnetic dipole (Glatzmaier *et al.* 1999;

Kutzner & Christensen 2004; Olson *et al.* 2010), as observed in our study.

### 4.1.2 Effect of non-zonal heat flux patterns

Olson *et al.* (2010) found an increased reversal frequency with the $Y_{1,1}$ pattern. Sahoo *et al.* (2016) on the contrary did not obtain reversals with weak amplitudes of the $Y_{1,1}$ pattern. Our results can reconcile both observations: dipolar stable dynamos are obtained for weak amplitudes of the $Y_{1,1}$ pattern, while a bistable dynamo is obtained for larger amplitudes. The $Y_{2,2}$ pattern has been shown to have a more moderate effect on the dynamo behaviour (Kutzner
& Christensen 2004; Olson et al. 2010; Sahoo et al. 2016). Nev-
ertheless, Glatzmaier et al. (1999) obtain a weaker and less stable
dipole with this pattern. We show that the effect of the Y₂,₂ pat-
tern largely depends on its amplitude. We find bistable dynamos
for moderate amplitudes, and multipolar dynamos for larger am-
plitudes. This suggests that the Y₂,₂ pattern slightly destabilizes the
dynamo, which becomes multipolar only for a large amplitude of
the pattern.

4.1.3 Comparison between high Rm and strong field models

The effects discussed above hold for the E1e-4_hRm case, which is
a standard case compared to previous dynamo studies as shown in
Fig. 1. We show in Fig. 7 that the effects of the ±Y2.0 patterns are
the same in the two other high Rm cases. The effect however differs
in the strong field cases. The –Y2.0 pattern is destabilizing in the
two strong field dynamos only for very large amplitudes, which are
probably unrealistic for the Earth. For more moderate amplitudes,
an equatorial cooling of the core strengthens the magnetic dipole.
The +Y2.0 pattern decreases the dipolar fraction without triggering
reversals. Olson & Amit (2014) showed that an increase in the re-
versal frequency can be obtained either for an equatorial cooling
or for a polar cooling depending on the geodynamo model. They
suggest that an equatorial cooling tends to bring the magnetic flux
towards the equator, thus destabilizing the axial dipole, while a po-
lar cooling tends to reinforce the underlying convection and thus
increase the role of inertia. The first effect, termed as ‘geographic
control’, would then hold for dynamos close to the onset of dy-
namo action (low Rm), while the second effect, termed as ‘inertial
control’, would hold for more turbulent dynamos (high Rm). We
find that an equatorial cooling destabilizes the magnetic dipole in
both high Rm models and in strong field models, suggesting that
a geographic control could be at play. In strong field dynamos, a
polar cooling of the core favours convection as shown by the in-
crease in the magnetic Reynolds number with the +Y2.0 pattern in
Table A4 for E1e-5_SF1 and Table A5 for E1e-5_SF2. However, this
strengthening of convection is not associated with a destabilization
of the dipole. This lack of inertial control on the dipole stability can
be explained by the fact that the strong field models are far from
the dipolar-multipolar transition in the dipolar regime, while the
models used by Olson & Amit (2014) lie closer to the transition.

## 4.2 Relation between thermal winds and the magnetic dipole stability

4.2.1 Interaction between zonal flows and the magnetic dipole

We find that the $Y_{2,0}$ pattern affects the zonal flow outside the tan-
gent cylinder through thermal winds, in agreement with previously
published models (Dietrich *et al.* 2017): a polar cooling ($+Y_{2,0}$)
favours a westward flow outside the tangent cylinder while an equa-
torial cooling ($-Y_{2,0}$) favours an eastward flow outside the tangent
cylinder. Previous studies showed that multipolar dynamo solutions
are associated with an increase in the zonal flows (Busse & Simitev
2011; Gastine *et al.* 2012; Schrinner *et al.* 2012). This observation
can be explained by the role of Maxwell stress: a strong dipolar mag-
netic field tends to counteract zonal flows and thus decrease the shear
that bends the field lines. Conversely, strong zonal flows increase
the shear of the poloidal magnetic field, generating an azimuthal
component of the magnetic field through an ω-effect (Roberts &
Glatzmaier 2000). We find that a strong polar cooling of the core

strengthen the ω-effect, which is associated with a decrease of the
dipolar fraction (Fig. 7). For the two strong field dynamo models, the
dipolar fraction is increased for a moderate equatorial cooling. This
could be associated in this case with a weakened ω-effect induced
by the weakening of the westward flows.

4.2.2 Disruption of zonal flows and dipolar–multipolar transition

The decrease of the dipolar fraction when increasing the amplitude of the +Y<sub>2,0</sub> pattern is smooth and is not associated with a destabilization of the magnetic dipole as observed with the −Y<sub>2,0</sub> pattern. Stanley (2010) showed that a degree 2 order 0 thermal heterogeneity can stabilize or destabilize the dipole if thermal winds generated by the boundary heterogeneity, respectively, strengthens or counteracts the zonal flows occurring without thermal heterogeneities. This is also what we observe: the dynamo becomes abruptly multipolar under the action of the −Y<sub>2,0</sub> pattern when the eastward thermal winds outside the tangent cylinder are strong enough to disrupt the otherwise dominantly westward zonal flows. As shown in Fig. 8, the transition imposed by the Y<sub>1,0</sub> and Y<sub>2,2</sub> patterns is also associated with a disruption of the westward flows. For these two patterns, the disruption mechanism does not rely on the weakening of westward flows but on the loss of equatorial symmetry in the zonal flows. We find that the −Y<sub>2,0</sub> pattern also favours antisymmetric flows for the high Rm dynamos, though the antisymmetry is less pronounced in the mean flows than with the Y<sub>1,0</sub> and Y<sub>2,2</sub> patterns. Such a loss of equatorial symmetry is expected with the Y<sub>1,0</sub> pattern that imposes equatorially antisymmetric thermal winds (Amit et al. 2011). However, the symmetry breaking is spontaneous in the case of the Y<sub>2,2</sub> and −Y<sub>2,0</sub> patterns as those patterns do not impose any antisymmetry. Solutions not dominated by equatorial symmetry have been previously obtained using the −Y<sub>2,0</sub> pattern by Cao et al. (2014) for dynamo models driven by a volumetric buoyancy source. Without solid inner-core, a spontaneous loss of equatorial symmetry of the mean flow has been observed with homogeneous boundary conditions in dynamo simulations (Landeau & Aubert 2011), and in double-diffusive convection (Monville et al. 2019).

To quantify the effect of the heat flux pattern on the equatorial
symmetry of the flow we define the ratio between the RMS anti-
symmetric zonal flow and the RMS of the full velocity field, called
⟨U<sup>a</sup><sub>Φ</sub⟩. In Fig. 13 we illustrate the two mechanisms that can lead
to a decreased stability of the magnetic dipole associated with the
disruption of the westward flows in our set of simulations. A first
mechanism is obtained with the Y<sub>1,0</sub> pattern (rightward triangles in
Fig. 13) that forcibly breaks the equatorial symmetry of the zonal
flows without affecting the space-averaged amplitude of the zonal
flow. A similar but spontaneous loss of equatorial symmetry is ob-
tained for a large amplitude of the Y<sub>2,2</sub> pattern (leftward triangles in
Fig. 13), associated with a transition to a multipolar behaviour (see
also Fig. 8). A second mechanism is obtained with the −Y<sub>2,0</sub> pat-
tern (dashed upward triangles in Fig. 13), which favours dynamos
that have weaker westward flows and are reversing or multipolar. In
high Rm dynamos, zonal flows become dominantly eastward under
the action of the −Y<sub>2,0</sub> pattern, with a significant loss of equatorial
symmetry in most cases. Note that this loss of equatorial symmetry
only occurs in the high Rm models, and that the dynamo obtained
with the E1e-5_hRm model and δq<sup>*</sup><sub>0</sub> = −2.9 per cent is multipolar
without showing a significantly increased equatorial antisymmetric
flow (see Table A4). On the other hand, thermal winds generated by
the +Y<sub>2,0</sub> pattern (plain upward triangles in Fig. 13) strengthen the
<figure><img src="image_15.png" /></figure>

**Figure 13.** Mean angular velocity of the flow as a function of the RMS antisymmetric zonal flow normalized by the RMS of the full velocity. Dynamos using the $Y_{2,0}$, $Y_{1,0}$ and $Y_{2,2}$ patterns are shown with the upward, rightward and leftward triangles, respectively. The reference dynamos with uniform heat flux conditions are shown as circles. The colour scale gives the average reversals period for the reversing dynamos. Dipolar stable dynamos are shown in blue, and multipolar dynamos are shown in orange. Dynamos within the reversing range defined in Fig. 12 ($10^{-5} \le \bar{M}_{Za}^* \le 10^{-3}$) are shown with grey contours. The dipolar stable dynamos outside this range in the bottom left of the figure have $\bar{M}_{Za}^* > 10^{-3}$, while the multipolar dynamos outside this range in the top right of the figure have $\bar{M}_{Za}^* < 10^{-5}$. The grey vertical dashed line for $\langle U_{\phi}^a \rangle = 0.19$ separates the range in which only multipolar dynamos are found (on the right-hand side). The orange and blue horizontal dashed lines for $\Omega_f = \pm 0.13$ delimitate the same region, shaded in grey, as the one in Fig. 10.

equatorially symmetric westward zonal flow, which is associated with non-reversing dipoles.

It has been shown in the context of giant planets that the direction of zonal flows depends on the forcing of convection (Gastine et al. 2013). In the geophysical context, westward zonal flows are found to naturally occur in numerical dynamos with homogeneous boundary conditions (Kuang & Bloxham 1997; Aubert 2005; Schaeffer et al. 2017). The competition between eastward thermal winds and westward flows could thus explain the widely observed destabilizing effect of the $-Y_{2,0}$ pattern in numerical models (Glatzmaier et al. 1999; Kutzner & Christensen 2004; Olson et al. 2010; Olson & Amit 2014). According to this mechanism, the $+Y_{2,0}$ pattern should not trigger reversals. As mentioned previously, the $+Y_{2,0}$ pattern has been nevertheless found to be destabilizing in some viscosity-dominated geodynamo models (Olson et al. 2010; Olson & Amit 2014). While certainly not the only way to destabilize a magnetic dipole, the destabilization mechanism associated with the disruption of the westward flows observed in our low viscosity simulations may be more relevant for the Earth's core.

### 4.2.3 Role of equatorial antisymmetry in the dipolar-multipolar transition

A decrease of the magnetic dipole stability has been found in geodynamo simulations when the equatorially antisymmetric flow is increased (Li et al. 2002; Olson et al. 2009, 2010; Garcia et al. 2017) as well as in the VKS experiment (Berhanu et al. 2010). Garcia et al. (2017) notably showed that the dipolar-multipolar transition in a set of simulations can be explained by a hydrodynamic transition associated with the loss of equatorial symmetry in

the flow. Pétrélis et al. (2009) suggested a reversal mechanism in which breaking the equatorial symmetry of the flow is key. Our results support a role of equatorial symmetry breaking in the dipolar-multipolar transition. Nevertheless, a transition to a multipolar state is obtained without significantly increasing the antisymmetry in the E1e-5\_SF1 model and for moderate amplitudes of the $-Y_{2,0}$ pattern in the E1e-5\_hRm model. Moreover, no systematic increase in the antisymmetric zonal flow is observed in the reversing dynamos. Thus, though the equatorial symmetry breaking seems to be often associated with a destabilization of the magnetic dipole in numerical models and in experiments, it is very likely not the only cause for such destabilization.

## 4.3 Geophysical interpretation

### 4.3.1 Destabilization mechanisms and westward flows

Geomagnetic observations support the existence of a dominant eccentric westward gyre, symmetric with respect to the equator, that dominates flow reconstructions at the core–mantle boundary (Pais & Jault 2008; Finlay et al. 2023). In our simulations, westward flows naturally occur and are well visible in the time-averaged meridional cuts (see Fig. 8 and Fig. A3). The westward zonal flow inside the core suggests that destabilizing the present-day dipole would require either an equatorially antisymmetric heat flux ($Y_{1,0}$ mechanism) or a strong equatorial cooling ($-Y_{2,0}$ mechanism). A reinforced polar cooling would, on the contrary, strengthen the westward flows and stabilize the dipole. The heat flux pattern as deduced from seismic tomography is dominated by an equatorially symmetric $Y_{2,2}$ component (Amit et al. 2015), which does not favour equatorial cooling over polar cooling on average. The effect of such pattern on the amplitude and symmetry of zonal flows needs further studies. Nevertheless, our results suggest that a tomographic pattern should be only weakly destabilizing for the dynamo.

Interestingly, the implications of these mechanisms can be, to a certain extent, tested for the Earth. The present-day symmetric westward gyre deduced from magnetic observations is compatible with a dipole-dominated dynamo. Magnetic field structures can be reconstructed from paleomagnetic data, and some persistence can be observed in the last thousand to million years (Johnson & Constable 1995; Constable et al. 2016). The recent study by Clizzie & Constable (2024) notably found that a westward drift of magnetic field structures dominates in the past 100 kyr and that the Laschamp excursion is associated with a transition from eastward drift to westward drift. This persistent westward drift could be a sign of a persistent westward flow, while the change in the drift direction suggests a relation between the dipole stability and changes in the zonal flows.

### 4.3.2 Parametrization of the dipolar-multipolar transition

We show in Fig. 11(a) that the local Rossby number, which was found to control the dipolar-multipolar transition in early dynamo models (Olson & Christensen 2006; Christensen 2010), does not discriminate between strongly dipolar and weakly dipolar dynamos within our simulation set. The inability for this parameter to predict the behaviour of the dynamo was already demonstrated (Pettidemange 2018; Menu et al. 2020). The ratio of magnetic to kinetic energy, put forward by Tassin et al. (2021), also fails to describe accurately the transition in our set of simulations as we
find multipolar dynamos with a magnetic to kinetic energy ratio as large as 6 (Fig. 11b). The dipolar–multipolar transition imposed by an equatorial cooling in the E1e-5.SF1 model is particularly interesting in this respect, as it occurs despite the magnetic energy staying significantly larger than the kinetic energy. Our results thus show that heat flux heterogeneities can be a way to destabilize the dipole while keeping a strong magnetic field, as it is expected for the Earth. Though such behaviour is not common in previously published dynamos, other choices in the model’s setup can be made to obtain unstable dipoles with large M, for example by considering larger Prandtl numbers as demonstrated by Jones & Tsang (2024).

We show that the parameter $M_{Za}^* = \frac{\Lambda_{Za}}{Rm_{Za}^2}$ (or alternatively $M_{Za}^* = \frac{E_{bZa}}{E_{uZa}} E_{\eta}$) describes more accurately the evolution of the dipolar fraction. We find a ‘reversing range’ for $10^{-5} < \overline{M}_{Za}^* < 10^{-3}$ outside of which no reversing dynamos are found. As shown in Fig. 13, this is within this range that dynamos with weak equatorially symmetric zonal flows are found (corresponding to the markers with grey contours). Dynamos with large $M_{Za}^*$ have strong equatorially symmetric westward flows, while dynamos with low $M_{Za}^*$ have equatorially antisymmetric zonal flows, suggesting that the captured transition is associated with the disruption of westward flows discussed in Section 4.2.

To obtain an estimate of this parameter for the Earth’s core, we can write $M \approx \frac{B_{RMS}^*}{\rho\mu_0 U_{RMS}^2}$. With $B_{RMS} = 4$ mT (Gillet et al. 2010), $U_{RMS} = 4 \times 10^{-4}$ m s$^{-1}$ (Finlay & Amit 2011), $\rho = 1.1 \times 10^4$ kg m$^{-3}$ (Olson 2015) and $E_{\eta} = 10^{-9}$, we have $M \sim 7 \times 10^3$. We can further assume that the zonal antisymmetric part of the kinetic energy is of about 1 per cent of the total kinetic energy, consistent with flow reconstructions (Fournier et al. 2011). Because the magnetic dipole has a zonal antisymmetric geometry and dominates the magnetic energy in the Earth’s core, the zonal antisymmetric part of the magnetic energy can be assumed to be not very different from the full magnetic energy. We will here consider that $0.1 \le E_{bZa}/E_{b} \le 1$. Using these values, we obtain $7 \times 10^{-5} \le M_{Za}^* \le 7 \times 10^{-4}$. This interval lies within the reversing range obtained for our set of simulations, consistent with a geodynamo that can reverse polarity. This suggests that the transition described by the $M_{Za}^*$ parameter could be relevant for the Earth, in contrast with all previously proposed criteria.

The correlation obtained between $\overline{f}_{\text{dip}}$ and $\overline{M}_{Za}^*$ can be understood by considering the two terms appearing in the definition of parameter $M_{Za}^*$ in eq. (28). The term $\Lambda_{Za}$ mostly quantifies the strength of the magnetic dipole that has a zonal equatorially antisymmetric geometry. It is thus expected to have an increased value of $M_{Za}^*$ when the dipolar fraction is increased. Similarly, $Rm_{Za}$ quantifies the strength of the equatorially antisymmetric zonal flows. Given the relation between the loss of equatorial symmetry and the multipolar transition, one can expect the dipolar fraction to decrease when $Rm_{Za}$ increases. It should be noted that a transition from a dipolar state towards a multipolar state is expected to be associated with a decrease of the ratio of magnetic to kinetic energy: the magnetic energy is decreased by the collapse of the dipole, and the kinetic energy can be expected to increase due to a reduced Lorentz force. We remark that all the criteria tested in Fig. 11 suffer from this inherent bias: multipolar dynamos are associated with lower magnetic field and increased velocities, so that combinations that include velocity and/or magnetic field intensity will naturally be discriminant of the transition. The $M_{Za}^*$ parameter nevertheless has the advantage of

locating the transition for similar values of $M_{Za}^*$ nearly independently of the considered reference geodynamo model. Nevertheless, all the reference models used in this study share important characteristics such as the mean heat flux, the radius of the inner core or the Prandtl number (fixed to unity in this study). The relevance of $M_{Za}^*$ when varying these parameters should be investigated. To make such a vast parametric study tractable, we plan to adopt promising rare event simulation algorithms (e.g. Bouchet et al. 2019) in future works.

## 5 CONCLUSION

We explore in this study the effect of heat flux heterogeneities in dynamo models ranging from standard parameter regimes to more Earth like regimes than previously explored. We show that heat flux heterogeneities can have different effects depending on the chosen dynamo model. We nevertheless obtain systematic behaviours, summarized below.

In our simulations, the destabilization of the magnetic dipole due to heat flux heterogeneities is associated with a disruption of the westward flows naturally occurring in the five reference models used in this study. As shown in Fig. 13, this disruption is caused either by a weakening of the westward flows through thermal winds imposed by an equatorial cooling of the core ($-Y_{2,0}$ pattern), or by a loss of equatorial symmetry in the zonal flows. This loss of equatorial symmetry can be imposed by the $Y_{1,0}$ heat flux pattern, or it can occur spontaneously when increasing the amplitude of the $Y_{2,2}$ or $-Y_{2,0}$ pattern. Conversely, a polar cooling of the core strengthens the symmetric westward flow, and stabilizes the dipole. Our results suggest that the destabilization of the magnetic dipole through heat flux heterogeneities is more difficult in dynamos with strong westward flows. However, we insist that it is not clear from our results whether the westward flow disruption is a cause or a consequence of the multipolar transition. Applied to the Earth’s core, this implies that an equatorial cooling of the Earth’s core or an equatorially antisymmetric heat flux pattern could disrupt the observed westward zonal flows and destabilize the magnetic dipole, thus increasing the reversal frequency. Conversely, a favoured polar cooling could increase the westward flows and stabilize the dipole.

Because the geodynamo is a highly nonlinear system, the effects of combinations of the simple patterns considered here are not obvious. We plan to explore the effect of complex heat flux patterns found in mantle convection simulations in a future study.

As shown in Fig. 11, we find that the local Rossby number and the magnetic to kinetic energy ratio, both previously put forward as controlling the dipolar–multipolar transition, do not capture the dynamo behaviours observed in our simulations. Instead, the dipolar fraction of the magnetic field is found to be correlated with $M_{Za}^* = \frac{\Lambda_{Za}}{Rm_{Za}^2}$ (see Fig. 12). Dynamos with high $\overline{M}_{Za}^*$ have strong and stable dipoles, while dynamos with low $\overline{M}_{Za}^*$ are multipolar. A transition is found for the range $10^{-5} < \overline{M}_{Za}^* < 10^{-3}$ in which all kinds of dynamos are found, including reversing ones. Estimates for the Earth fall into this range too, supporting the geophysical relevance of this criterion. This parameter is related to the strength and equatorial symmetry of the westward flows, suggesting that geodynamo reversals are associated to the disruption of westward flows. According to the relation we obtain between $M_{Za}^*$ and the dynamo behaviour, a superchron-type behaviour would occur for larger values of $M_{Za}^*$. Such increase in $M_{Za}^*$ could arise from a
decrease of the antisymmetric zonal flow component, or from an
increase in the magnetic field intensity.

ACKNOWLEDGMENTS

This research has been supported by the Agence Nationale de la
Recherche (grant no. ANR-19-CE31-0019). We thank Julien Aubert
and Johannes Wicht for thorough reviews on the first version of the
manuscript. We also thank Hagay Amit and Alexandre Fournier for
their comments. Numerical simulations have been performed us-
ing HPC resources of IDRIS, TGCC and CINES under allocations
A0100407382, A0120407382, A0140407382 and A0160407382
attributed by GENCI (Grand Equipement National de Calcul In-
tensif). The XSHELLS code received funding from the European
High-Performance Computing Joint Undertaking (JU) under Grant
Agreement No. 101093038 (ChEESE-2P)

DATA AVAILABILITY

The XSHELLS code used to run the geodynamo simulations is freely
available at https://nschaeff.bitbucket.io/xshells/.

REFERENCES

Amit, H., Christensen, U. R. & Langlais, B., 2011. The influence of degree-1 mantle heterogeneity on the past dynamo of Mars, *Phys. Earth planet. Inter.*, **189**(1-2), 63-79.

Amit, H., Choblet, G., Olson, P., Monteux, J., Deschamps, F., Langlais, B. & Tobie, G., 2015. Towards more realistic core-mantle boundary heat flux patterns: a source of diversity in planetary dynamos, *Prog. Earth planet. Sci.*, **2**(1), 26.

Aubert, J., 2005. Steady zonal flows in spherical shell dynamos, *J. Fluid Mech.* **542**, 53–67.

Aubert, J., Finlay, C.C. & Fournier, A., 2013. Bottom-up control of geo-
magnetic secular variation by the earth's inner core, *Nature*, **502**(7470),
219–223.

Aubert, J., Gastine, T. & Fournier, A., 2017. Spherical convective dynamos in the rapidly rotating asymptotic regime, *J. Fluid Mech.*, **813**, 558–593.

Berhanu, M. et al., 2009. Bistability between a stationary and an oscillatory dynamo in a turbulent flow of liquid sodium, *J. Fluid Mech.*, **641**, 217–226.

Berhanu, M. et al., 2010. Dynamo regimes and transitions in the VKS experiment, *Eur. Phys. J. B*, **77**, 459–468.

Biggin, A.J., Steinberger, B., Aubert, J., Suttie, N., Holme, R., Torsvik, T.H., van der Meer, D.G. & van Hinsbergen, D. J.J., 2012. Possible links between long-term geomagnetic variations and whole-mantle convection processes, *Nat. Geosci.*, **5**(8), 526–533.

Bouchet, F., Rolland, J. & Simonnet, E., 2019. Rare event algorithm links transitions in turbulent flows with activated nucleations, *Phys. Rev. Lett.*, **122**, 074502, doi: 10.1103/PhysRevLett.122.074502.

Braginsky, S.I. & Roberts, P.H., 1995. Equations governing convection in earth's core and the geodynamo, *Geophys. Astrophys. Fluid Dyn.*, **79**(1-4), 1-97.

Buffett, B.A., Ziegler, L. & Constable, C.G., 2013. A stochastic model for palaeomagnetic field variations, *Geophys. J. Int.*, **195**(1), 86-97.

Burke, K., Steinberger, B., Torsvik, T.H. & Smethurst, M.A., 2008. Plume generation zones at the margins of large low shear velocity provinces on the core-mantle boundary, *Earth planet. Sci. Lett.*, **265**(1-2), 49-60.

Busse, F.H. & Simitev, R.D., 2011. Remarks on some typical assumptions in dynamo theory, *Geophys. Astrophys. Fluid Dyn.*, **105**(2-3), 234-247.

Cao, H., Aurnou, J.M., Wicht, J., Dietrich, W., Soderlund, K.M. & Russell, C.T., 2014. A dynamo explanation for mercury's anomalous magnetic field, *Geophys. Res. Lett.*, **41**(12), 4127-4134.

Christensen, U.R., 2010. Dynamo scaling laws and applications to the planets, *Space Sci. Rev.* **152**(1-4), 565-590.

Christensen, U.R. & Aubert, J., 2006. Scaling properties of convection-
driven dynamos in rotating spherical shells and application to planetary
magnetic fields, *Geophys. J. Int.*, **166**(1), 97–114.

Clizzie, N. & Constable, C., 2024. Reversal of drift direction during the laschamp geomagnetic excursion, *Phys. Earth planet. Inter.*, **347**, 107143, doi: 10.1016/j.pepi.2024.107143.

Constable, C., Korte, M. & Panovska, S., 2016. Persistent high paleosecular variation activity in southern hemisphere for at least 10 000 years, *Earth planet. Sci. Lett.*, **453**, 78–86.

Dannberg, J., Gassmöller, R., Thallner, D., LaCombe, F. & Sprain, C., 2024. Changes in core-mantle boundary heat flux patterns throughout the supercontinent cycle, *Geophys. J. Int.*, **237**(3), 1251-1274.

Davies, C.J., Bono, R.K., Meduri, D.G., Aubert, J., Greenwood, S. & Biggin, A.J., 2022. Dynamo constraints on the long-term evolution of earth's magnetic field strength, *Geophys. J. Int.*, **228**(1), 316–336.

Dietrich, W., Hori, K. & Wicht, J., 2016. Core flows and heat transfer induced by inhomogeneous cooling with sub-and supercritical convection, *Phys. Earth planet. Inter.*, **251**, 36–51.

Dietrich, W., Gastine, T. & Wicht, J., 2017. Reversal and amplification of zonal flows by boundary enforced thermal wind, *Icarus*, **282**, 380–392.

Dormy, E., 2016. Strong-field spherical dynamos, *J. Fluid Mech.*, **789**, 500-513.

Driscoll, P. & Olson, P., 2009a. Effects of buoyancy and rotation on the polarity reversal frequency of gravitationally driven numerical dynamos, *Geophys. J. Int.*, **178**(3), 1337-1350.

Driscoll, P. & Olson, P., 2009b. Polarity reversals in geodynamo models with core evolution, *Earth planet. Sci. Lett.*, **282**(1-4), 24-33.

Dziewonski, A.M., Lekic, V. & Romanowicz, B.A., 2010. Mantle anchor structure: an argument for bottom up tectonics, *Earth planet. Sci. Lett.*, **299**(1-2), 69-79.

Finlay, C.C. & Amit, H., 2011. On flow magnitude and field-flow alignment at earth's core surface, *Geophys. J. Int.*, **186**(1), 175-192.

Finlay, C.C., Gillet, N., Aubert, J., Livermore, P.W. & Jault, D., 2023. Gyres, jets and waves in the earth's core, *Nat. Rev. Earth Environ.*, **4**(6), 377-392.

Flament, N., Bodur, Ö.F., Williams, S.E. & Merdith, A.S., 2022. Assembly of the basal mantle structure beneath Africa, *Nature*, **603**(7903), 846-851.

Fournier, A., Aubert, J. & Thébault, E., 2011. Inference on core surface flow from observations and 3-d dynamo modelling, *Geophys. J. Int.*, **186**(1), 118-136.

Frasson, T., Labrosse, S., Nataf, H.-C., Coltice, N. & Flament, N., 2024. On the impact of true polar wander on heat flux patterns at the core-mantle boundary, *Solid Earth*, **15**(5), 617-637.

Garcia, F., Oruba, L. & Dormy, E., 2017. Equatorial symmetry breaking and the loss of dipolarity in rapidly rotating dynamos, *Geophys. Astrophys. Fluid Dyn.*, **111**(5), 380-393.

Garnero, E.J. & McNamara, A.K., 2008. Structure and dynamics of earth's lower mantle, *Science*, **320**(5876), 626-628.

Gastine, T., Duarte, L. & Wicht, J., 2012. Dipolar versus multipolar dynamos: the influence of the background density stratification, *Astron. Astrophys.*, **546**, A19, doi: 10.1051/0004-6361/201219799.

Gastine, T., Wicht, J. & Aurnou, J., 2013. Zonal flow regimes in rotating anelastic spherical shells: An application to giant planets, *Icarus*, **225**(1), 156-172.

Gillet, N., Jault, D., Canet, E. & Fournier, A., 2010. Fast torsional waves and strong magnetic field within the earth's core, *Nature*, **465**(7294), 74-77.

Glatzmaier, G.A. & Roberts, P.H., 1995. A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle, *Phys. Earth planet. Inter.*, **91**(1-3), 63-75.

Glatzmaier, G.A., Coe, R.S., Hongre, L. & Roberts, P.H., 1999. The role of the Earth's mantle in controlling the frequency of geomagnetic reversals, *Nature*, **401**(6756), 885-890.

Guervilly, C., Cardin, P. & Schaeffer, N., 2019. Turbulent convective length scale in planetary cores, *Nature*, **570**(7761), 368-371.

Holdenried-Chernoff, D. & Buffett, B., 2022. Evidence for turbulent magnetic diffusion in earth's core, *Geochem. Geophys. Geosyst.*, **23**(12), e2022GC010672, doi: 10.1029/2022GC010672.
Hounslow, M.W., Domeier, M. & Biggin, A.J., 2018. Subduction flux modulates the geomagnetic polarity reversal rate, *Tectonophysics*, **742–743**, 34–49.

Johnson, C.L. & Constable, C.G., 1995. The time-averaged geomagnetic field as recorded by lava flows over the past 5 myr, *Geophys. J. Int.*, **122**(2), 489–519.

Jones, C. & Tsang, Y.-K., 2024. Low inertia reversing geodynamos, preprint(arXiv:2408.07420).

Kuang, W. & Bloxham, J., 1997. An earth-like numerical dynamo model, *Nature*, **389**(6649), 371–374.

Kutzner, C. & Christensen, U., 2002. From stable dipolar towards reversing numerical dynamos, *Phys. Earth planet. Inter.*, **131**(1), 29–45.

Kutzner, C. & Christensen, U.R., 2004. Simulated geomagnetic reversals and preferred virtual geomagnetic pole paths, *Geophys. J. Int.*, **157**(3), 1105–1118.

Landeau, M. & Aubert, J., 2011. Equatorially asymmetric convection inducing a hemispherical magnetic field in rotating spheres and implications for the past Martian dynamo, *Phys. Earth planet. Inter.*, **185**(3–4), 61–73.

Li, J., Sato, T. & Kageyama, A., 2002. Repeated and sudden reversals of the dipole field generated by a spherical dynamo action, *Science*, **295**(5561), 1887–1890.

Loper, D.E., 1978. The gravitationally powered dynamo, *Geophys. J. Int.*, **54**(2), 389–404.

Lowrie, W. & Kent, D.V., 2013. Geomagnetic polarity timescales and reversal frequency regimes, in *Geophysical Monograph Series*, pp. 117–129, eds, Channell, J., Kent, D., Lowrie, W. & Meert, J., American Geophysical Union.

Matsui, H. et al., 2016. Performance benchmarks for a next generation numerical dynamo model, *Geochem. Geophys. Geosyst.*, **17**(5), 1586–1607.

Meduri, D., Biggin, A., Davies, C., Bono, R., Sprain, C. & Wicht, J., 2021. Numerical dynamo simulations reproduce paleomagnetic field behavior, *Geophys. Res. Lett.*, **48**(5), e2020GL090544, doi: 10.1029/2020GL090544.

Menu, M.D., Petitdemange, L. & Galtier, S., 2020. Magnetic effects on fields morphologies and reversals in geodynamo simulations, *Phys. Earth planet. Inter.*, **307**, 106542, doi: 10.1016/j.pepi.2020.106542.

Monville, R., Vidal, J., Cébron, D. & Schaeffer, N., 2019. Rotating double-diffusive convection in stably stratified planetary cores, *Geophys. J. Int.*, **219**(S1), S195–S218.

Mosca, I., Cobden, L., Deuss, A., Ritsema, J. & Trampert, J., 2012. Seismic and mineralogical structures of the lower mantle from probabilistic tomography, *J. geophys. Res.*, **117**(B6). doi:

Mound, J.E. & Davies, C.J., 2023. Longitudinal structure of Earth's magnetic field controlled by lower mantle heat flow, *Nat. Geosci.*, **16**(4), 380–385.

Nakagawa, T. & Davies, C.J., 2022. Combined dynamical and morphological characterisation of geodynamo simulations, *Earth planet. Sci. Lett.*, **594**, 117752, doi: 10.1016/j.epsl.2022.117752.

Nakagawa, T. & Tackley, P.J., 2008. Lateral variations in CMB heat flux and deep mantle seismic velocity caused by a thermal–chemical-phase boundary layer in 3D spherical convection, *Earth planet. Sci. Lett.*, **271**(1–4), 348–358.

Nataf, H.-C. & Schaeffer, N., 2015. Turbulence in the core, in *Treatise on Geophysics (Second Edition)*, pp. 161–181, ed. Schubert, G., Elsevier.

Nataf, H.-C. & Schaeffer, N., 2024. Dynamic regimes in planetary cores: τ–l diagrams, *C. R. Géosci.*, **356**(G1), 1–30. doi:

Olson, P., 2015. 8.01 - Core dynamics: an introduction and overview, in *Treatise on Geophysics (Second Edition)*, pp. 1–25, ed. Schubert, G., Elsevier, Oxford.

Olson, P. & Amit, H., 2014. Magnetic reversal frequency scaling in dynamos with thermochemical convection, *Phys. Earth planet. Inter.*, **229**, 122–133.

Olson, P. & Christensen, U.R., 2002. The time-averaged magnetic field in numerical dynamos with non-uniform boundary heat flow, *Geophys. J. Int.*, **151**(3), 809–823.

Olson, P. & Christensen, U.R., 2006. Dipole moment scaling for convection-driven planetary dynamos, *Earth planet. Sci. Lett.*, **250**(3–4), 561–571.

Olson, P., Driscoll, P. & Amit, H., 2009. Dipole collapse and reversal precursors in a numerical dynamo, *Phys. Earth planet. Inter.*, **173**(1–2), 121–140.

Olson, P.L., Coe, R.S., Driscoll, P.E., Glatzmaier, G.A. & Roberts, P.H., 2010. Geodynamo reversal frequency and heterogeneous core–mantle boundary heat flow, *Phys. Earth planet. Inter.*, **180**(1–2), 66–79.

Olson, P.L., Glatzmaier, G.A. & Coe, R.S., 2011. Complex polarity reversals in a geodynamo model, *Earth planet. Sci. Lett.*, **304**(1–2), 168–179.

Pais, M. & Jault, D., 2008. Quasi-geostrophic flows responsible for the secular variation of the Earth's magnetic field, *Geophys. J. Int.*, **173**(2), 421–443.

Petitdemange, L., 2018. Systematic parameter study of dynamo bifurcations in geodynamo simulations, *Phys. Earth planet. Inter.*, **277**, 113–132.

Pétrélis, F., Fauve, S., Dormy, E. & Valet, J.-P., 2009. Simple mechanism for reversals of earth's magnetic field, *Phys. Rev. Lett.*, **102**(14), 144503, doi: 10.1103/PhysRevLett.102.144503.

Pozzo, M., Davies, C., Gubbins, D. & Alfè, D., 2014. Thermal and electrical conductivity of solid iron and iron–silicon mixtures at earth’s core conditions, *Earth planet. Sci. Lett.*, **393**, 159–164.

Roberts, P.H. & Glatzmaier, G.A., 2000. Geodynamo theory and simulations, *Revi. Mod. Phys.*, **72**(4), 1081.

Sahoo, S., Sreenivasan, B. & Amit, H., 2016. Dynamos driven by weak thermal convection and heterogeneous outer boundary heat flux, *Phys. Earth planet. Inter.*, **250**, 35–45.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations, *Geochem. Geophys. Geosyst.*, **14**(3), 751–758.

Schaeffer, N., Jault, D., Nataf, H.-C. & Fournier, A., 2017. Turbulent geodynamo simulations: a leap towards Earth's core, *Geophys. J. Int.*, **211**(1), 1–29.

Schrinner, M., Petitdemange, L. & Dormy, E., 2012. Dipole collapse and dynamo waves in global direct numerical simulations, *Astrophys. J.*, **752**(2), 121, doi: 10.1088/0004-637X/752/2/121.

Schwaiger, T., Gastine, T. & Aubert, J., 2021. Relating force balances and flow length scales in geodynamo simulations, *Geophys. J. Int.*, **224**(3), 1890–1904.

Simitev, R.D. & Busse, F.H., 2009. Bistability and hysteresis of dipolar dynamos generated by turbulent convection in rotating spherical shells, *Europhys. Lett.*, **85**(1), 19001, doi: 10.1209/0295-5075/85/19001.

Soderlund, K.M., King, E.M. & Aurnou, J.M., 2012. The influence of magnetic fields in planetary dynamo models, *Earth planet. Sci. Lett.*, **333**, 9–20.

Sreenivasan, B. & Jones, C.A., 2006. Azimuthal winds, convection and dynamo action in the polar regions of planetary cores, *Geophys. Astrophys. Fluid Dyn.*, **100**(4–5), 319–339.

Stanley, S., 2010. A dynamo model for axisymmetrizing saturn's magnetic field, *Geophys. Res. Lett.*, **37**(5), doi: 10.1029/2009GL041752.

Tassin, T., Gastine, T. & Fournier, A., 2021. Geomagnetic semblance and dipolar–multipolar transition in top-heavy double-diffusive geodynamo models, *Geophys. J. Int.*, **226**(3), 1897–1919.

Terra-Nova, F. & Amit, H., 2024. Regionally-triggered geomagnetic reversals, *Sci. Rep.*, **14**(1), 9639, doi: 10.1038/s41598-024-59849-z.

Torsvik, T.H., Burke, K., Steinberger, B., Webb, S.J. & Ashwal, L.D., 2010. Diamonds sampled by plumes from the core–mantle boundary, *Nature*, **466**(7304), 352–355.

Trampert, J., 2004. Probabilistic tomography maps chemical heterogeneities throughout the lower mantle, *Science*, **306**(5697), 853–856.

Wicht, J. & Meduri, D.G., 2016. A Gaussian model for simulated geomagnetic field reversals, *Phys. Earth planet. Inter.*, **259**, 45–60.

Wicht, J. & Tilgner, A., 2010. Theory and modeling of planetary dynamos, *Space Sci. Rev.*, **152**(1), 501–542.

Yadav, R.K., Gastine, T., Christensen, U.R., Wolk, S.J. & Poppenhaeger, K., 2016. Approaching a realistic force balance in geodynamo simulations, *Proc. Natl. Acad. Sci. U.S.A.*, **113**(43), 12 065–12 070.

Zhang, K. & Gubbins, D., 1992. On convection in the earth's core driven by lateral temperature variations in the lower mantle, *Geophys. J. Int.*, **108**(1), 247–255.
Zhang, N. & Zhong, S., 2011. Heat fluxes at the Earth's surface and core-mantle boundary since Pangea formation and their implications for the geomagnetic superchrons, *Earth planet. Sci. Lett.*, **306**(3–4), 205–216.

## APPENDIX A: SIMULATION DETAILS

We give in the following tables some details on the simulations used in this study. The spatial discretization of the numerical simulation is described in the radial direction by the number of radial points $N_r = N_r^{OC} + N_r^{IC}$ with $N_r^{OC}$ the number of points in the outer core and $N_r^{IC}$ the number of points in the inner core. The scalar fields are expanded in spherical harmonics in the horizontal direction up to the maximum degree $l_{max}$ and the maximum order $m_{max}$. Hyperdiffusion is used in two of the five reference dynamo models on the flow, the magnetic field and the codensity field. This hyperdiffusion follows the formula suggested by Nataf & Schaeffer (2015). For the flow, the viscosity is given as a function of the spherical harmonic degree $l$ by

$$v(l) = \begin{cases} v_0 & \text{for } l \le l_c \\ v_0 q_v^{l-l_c} & \text{for } l > l_c \end{cases} \qquad (A1)$$

with $q_v = \left(\frac{v_{max}}{v_0}\right) \frac{1}{l_{max} - l_c}$. Hyperdiffusion is parametrized by the ratio $\frac{v_{max}}{v_0}$ and the cut-off degree $l_c$, which are given in the tables

for the simulations using hyperdiffusion. Similarly, the magnetic diffusion is given by

$$\eta(l) = \begin{cases} \eta_0 & \text{for } l \le l_c \\ \eta_0 q_\eta^{l-l_c} & \text{for } l > l_c \end{cases} \qquad (A2)$$

with $q_\eta = \left(\frac{\eta_{max}}{\eta_0}\right) \frac{1}{l_{max} - l_c}$. Finally, the thermal diffusion follows

$$\kappa(l) = \begin{cases} \kappa_0 & \text{for } l \le l_c \\ \kappa_0 q_\kappa^{l-l_c} & \text{for } l > l_c \end{cases} \qquad (A3)$$

with $q_\kappa = \left(\frac{\kappa_{max}}{\kappa_0}\right) \frac{1}{l_{max} - l_c}$. The cut-off degree $l_c$ is the same for the magnetic field and the velocity field.

Table A1 gives the results for the E1e-4.hRm model, Table A2 for the E1e-5.hRm model, Table A3 for the E1e-6.hRm model, Table A4 for the E1e-5.SF1 model and Table A5 for the E1e-5.SF2 model.

<table><caption>**Table A1.** Characteristics of the simulations using the E1e-4.hRm reference case with $E = 10^{-4}$, $Ra = 1.5 \times 10^8$, $Pm = 3$ and $Pr = 1$. The spatial discretization is set using $N_r^{OC} = 100$, $N_r^{IC} = 44$, $l_{max} = 100$ and $m_{max} = 85$. Hyperdiffusivity has been used for these simulations choosing $\frac{v_{max}}{v_0} = 15$, $\frac{\eta_{max}}{\eta_0} = 45$, $\frac{\kappa_{max}}{\kappa_0} = 15$ and $l_c = 70$. The column $\delta q_0^*$ gives the amplitude of the pattern, $\Delta t_{sim}$ gives the duration of the simulation (after removal of the initial transition time), $f_{dip}$ gives the time-averaged value of the dipolar fraction at the CMB, $\tilde{\sigma}(f_{dip})$ gives the standard deviation of the dipolar fraction normalized by the mean, $\bar{M}$ gives the time-averaged energy ratio, $\bar{M}_{Za}^*$ gives the parameter define in eq. (28), $N_{cross}$ gives the number of times the magnetic dipole crosses the equator, $\bar{Rm}$ gives the time-averaged magnetic Reynolds number, $\bar{Ro}_i$ gives the time-averaged local Rossby number, $\bar{\Omega}_f$ gives the time-averaged mean angular velocity and $\bar{\langle U_\phi^a \rangle}$ gives the time-averaged antisymmetric azimuthal part of the flow. The values of $\bar{\Omega}_f$ that have not been computed are marked as 'NC'. The behaviours are DS for dipolar stable, R for reversing, B for bistable and M for multipolar.</caption><thead><tr><th>Pattern</th><th>δq<sub>0</sub><sup>*</sup><br>(per cent)</th><th>Δt<sub>sim</sub><br>(t<sub>η</sub>)</th><th>f<sub>dip</sub></th><th>δ̃(f<sub>dip</sub>)<br>(per cent)</th><th>M̄</th><th>M̄<sub>Za</sub><br>(×10<sup>-4</sup>)</th><th>N<sub>cross</sub></th><th>Rm̄</th><th>Ro<sub>i</sub><br>(×10<sup>-2</sup>)</th><th>Ω̃<sub>f</sub><br>(per cent)</th><th>⟨U<sub>φ</sub><sup>a</sup>⟩<br>(per cent)</th><th>Behaviour</th></tr></thead><tbody><tr><td>—</td><td>0.00</td><td>134.3</td><td>0.45</td><td>19.0</td><td>0.95</td><td>1.36</td><td>34</td><td>649</td><td>9.41</td><td>-1.2</td><td>10.4</td><td>R</td></tr><tr><td>Y10</td><td>0.72</td><td>35.8</td><td>0.39</td><td>25.7</td><td>0.91</td><td>0.55</td><td>8</td><td>649</td><td>9.70</td><td>-0.6</td><td>15.0</td><td>R</td></tr><tr><td>Y10</td><td>1.45</td><td>34.1</td><td>0.20</td><td>53.1</td><td>0.57</td><td>0.06</td><td>148</td><td>703</td><td>10.40</td><td>1.3</td><td>28.4</td><td>M</td></tr><tr><td>Y10</td><td>2.90</td><td>33.4</td><td>0.13</td><td>42.0</td><td>0.29</td><td>0.01</td><td>285</td><td>766</td><td>10.74</td><td>1.8</td><td>35.9</td><td>M</td></tr><tr><td>Y10</td><td>7.25</td><td>33.6</td><td>0.11</td><td>46.8</td><td>0.22</td><td>0.01</td><td>256</td><td>848</td><td>9.69</td><td>-0.9</td><td>47.3</td><td>M</td></tr><tr><td>Y11</td><td>0.72</td><td>38.3</td><td>0.47</td><td>11.1</td><td>1.01</td><td>1.94</td><td>0</td><td>651</td><td>9.06</td><td>NC</td><td>9.3</td><td>DS</td></tr><tr><td>Y11</td><td>1.45</td><td>39.3</td><td>0.49</td><td>9.3</td><td>1.05</td><td>2.24</td><td>0</td><td>653</td><td>8.84</td><td>NC</td><td>9.2</td><td>DS</td></tr><tr><td>Y11</td><td>2.90</td><td>35.7</td><td>0.49</td><td>10.0</td><td>0.97</td><td>1.97</td><td>0</td><td>667</td><td>8.74</td><td>NC</td><td>9.6</td><td>DS</td></tr><tr><td>Y11</td><td>7.25</td><td>39.0</td><td>0.30</td><td>38.7</td><td>0.65</td><td>0.3</td><td>81</td><td>721</td><td>9.14</td><td>NC</td><td>17.2</td><td>B</td></tr><tr><td>Y20</td><td>0.72</td><td>37.8</td><td>0.48</td><td>8.8</td><td>1.06</td><td>2.27</td><td>0</td><td>635</td><td>9.34</td><td>NC</td><td>8.8</td><td>DS</td></tr><tr><td>Y20</td><td>1.45</td><td>38.0</td><td>0.49</td><td>7.1</td><td>1.17</td><td>2.79</td><td>0</td><td>622</td><td>9.35</td><td>NC</td><td>8.9</td><td>DS</td></tr><tr><td>Y20</td><td>2.90</td><td>37.0</td><td>0.49</td><td>6.0</td><td>1.26</td><td>3.16</td><td>0</td><td>605</td><td>9.58</td><td>NC</td><td>9.2</td><td>DS</td></tr><tr><td>Y20</td><td>7.25</td><td>35.6</td><td>0.40</td><td>6.4</td><td>1.42</td><td>4.5</td><td>0</td><td>600</td><td>9.84</td><td>-14.8</td><td>9.7</td><td>DS</td></tr><tr><td>Y20</td><td>-0.72</td><td>35.5</td><td>0.21</td><td>68.2</td><td>0.53</td><td>0.09</td><td>265</td><td>706</td><td>10.31</td><td>NC</td><td>23.8</td><td>M</td></tr><tr><td>Y20</td><td>-1.45</td><td>35.2</td><td>0.13</td><td>48.9</td><td>0.37</td><td>0.03</td><td>354</td><td>741</td><td>10.56</td><td>NC</td><td>28.6</td><td>M</td></tr><tr><td>Y20</td><td>-2.90</td><td>35.7</td><td>0.11</td><td>42.5</td><td>0.32</td><td>0.02</td><td>399</td><td>771</td><td>10.45</td><td>8.1</td><td>30.8</td><td>M</td></tr><tr><td>Y20</td><td>-7.25</td><td>37.4</td><td>0.10</td><td>43.8</td><td>0.24</td><td>0.01</td><td>611</td><td>846</td><td>10.04</td><td>13.3</td><td>32.2</td><td>M</td></tr><tr><td>Y22</td><td>0.72</td><td>36.3</td><td>0.46</td><td>12.4</td><td>0.99</td><td>1.78</td><td>0</td><td>646</td><td>9.33</td><td>-1.4</td><td>9.3</td><td>DS</td></tr><tr><td>Y22</td><td>1.45</td><td>35.2</td><td>0.33</td><td>49.1</td><td>0.77</td><td>0.36</td><td>119</td><td>674</td><td>9.83</td><td>-0.2</td><td>17.5</td><td>B</td></tr><tr><td>Y22</td><td>2.90</td><td>36.2</td><td>0.37</td><td>37.5</td><td>0.81</td><td>0.58</td><td>58</td><td>670</td><td>9.54</td><td>-1.4</td><td>14.4</td><td>B</td></tr><tr><td>Y22</td><td>7.25</td><td>34.6</td><td>0.18</td><td>62.3</td><td>0.51</td><td>0.1</td><td>208</td><td>710</td><td>10.28</td><td>-0.0</td><td>20.9</td><td>M</td></tr></tbody></table>
<table><caption>**Table A2.** Characteristics of the simulations using the E1e-5.hRm reference case with $E = 10^{-5}$, $Ra = 5 \times 10^{9}$, $Pm = 2$ and $Pr = 1$. The spatial discretization is set using $N_r^{OC} = 280$, $N_r^{IC} = 80$, $l_{max} = 213$ and $m_{max} = 170$. The behaviours are DS for dipolar stable, R for reversing, B for bistable and M for multipolar.</caption><thead><tr><th>Pattern</th><th>δq<sup>*</sup><sub>o</sub><br>(per cent)</th><th>Δt<sub>sim</sub><br>(t<sub>η</sub>)</th><th>∫<sub>dip</sub></th><th>δ̂(f<sub>dip</sub>)<br>(per cent)</th><th>M̄<br>(×10<sup>-4</sup>)</th><th>N<sub>cross</sub></th><th>Rm̄</th><th>Ro<sub>l</sub><br>(×10<sup>-2</sup>)</th><th>Ω<sub>f</sub><br>(per cent)</th><th>⟨U<sup>a</sup><sub>φ</sub>⟩<br>(per cent)</th><th>Behaviour</th></tr></thead><tbody><tr><td>–</td><td>0.00</td><td>96.0</td><td>0.38</td><td>25.2</td><td>3.49</td><td>0.94</td><td>17</td><td>987</td><td>4.36</td><td>-5.9</td><td>8.8</td><td>R</td></tr><tr><td>Y20</td><td>2.90</td><td>1.9</td><td>0.47</td><td>8.4</td><td>5.67</td><td>2.3</td><td>0</td><td>913</td><td>3.74</td><td>-8.7</td><td>9.5</td><td>DS</td></tr><tr><td>Y20</td><td>14.50</td><td>0.8</td><td>0.35</td><td>6.6</td><td>13.3</td><td>14.78</td><td>0</td><td>880</td><td>3.18</td><td>-36.3</td><td>9.5</td><td>DS</td></tr><tr><td>Y20</td><td>-2.90</td><td>4.7</td><td>0.18</td><td>51.3</td><td>2.27</td><td>0.2</td><td>33</td><td>1058</td><td>4.87</td><td>-1.5</td><td>10.9</td><td>M</td></tr><tr><td>Y20</td><td>-14.50</td><td>5.3</td><td>0.08</td><td>40.2</td><td>1.43</td><td>0.03</td><td>125</td><td>1314</td><td>4.57</td><td>17.5</td><td>22.6</td><td>M</td></tr></tbody></table>

<table><caption>**Table A3.** Characteristics of the simulations using the E1e-6.hRm reference case with $E = 10^{-6}$, $Ra = 9.6 \times 10^{11}$, $Pm = 0.2$ and $Pr = 1$. Hyperdiffusivity has been used for these simulations. The spatial discretization is set using $N_r^{OC} = 512$, $N_r^{IC} = 0$ (insulating inner core), $l_{max} = 341$ and $m_{max} = 341$. Hyperdiffusivity has been used for these simulations choosing $\frac{V_{max}}{V_0} = 5$, $\frac{\eta_{max}}{\eta_0} = 1$ (no hyperdiffusion on the magnetic field), $\frac{\kappa_{max}}{\eta_0} = 5$ and $l_c = 307$. The values of $\Omega_f$ that have not been computed are marked as 'NC'. The behaviours are DS for dipolar stable, R for reversing, B for bistable and M for multipolar.</caption><thead><tr><th>Pattern</th><th>δq<sup>*</sup><sub>o</sub><br>(per cent)</th><th>Δt<sub>sim</sub><br>(t<sub>η</sub>)</th><th>∫<sub>dip</sub></th><th>δ̂(f<sub>dip</sub>)<br>(per cent)</th><th>M̄<br>(×10<sup>-4</sup>)</th><th>N<sub>cross</sub></th><th>Rm̄</th><th>Ro<sub>l</sub><br>(×10<sup>-2</sup>)</th><th>Ω<sub>f</sub><br>(per cent)</th><th>⟨U<sup>a</sup><sub>φ</sub>⟩<br>(per cent)</th><th>Behaviour</th></tr></thead><tbody><tr><td>–</td><td>0.00</td><td>3.7</td><td>0.72</td><td>7.4</td><td>3.53</td><td>1.06</td><td>0</td><td>575</td><td>3.93</td><td>1.3</td><td>14.4</td><td>DS</td></tr><tr><td>Y20</td><td>9.06</td><td>4.2</td><td>0.60</td><td>10.5</td><td>9.68</td><td>11.85</td><td>0</td><td>476</td><td>3.20</td><td>-27.0</td><td>8.7</td><td>DS</td></tr><tr><td>Y20</td><td>-1.81</td><td>15.3</td><td>0.41</td><td>44.1</td><td>2.26</td><td>0.33</td><td>51</td><td>622</td><td>4.79</td><td>NC</td><td>17.1</td><td>B</td></tr><tr><td>Y20</td><td>-2.54</td><td>5.8</td><td>0.22</td><td>44.3</td><td>1.18</td><td>0.03</td><td>74</td><td>678</td><td>5.51</td><td>NC</td><td>24.4</td><td>M</td></tr><tr><td>Y20</td><td>-2.90</td><td>10.8</td><td>0.33</td><td>51.8</td><td>2.06</td><td>0.17</td><td>74</td><td>630</td><td>4.92</td><td>NC</td><td>19.1</td><td>B</td></tr><tr><td>Y20</td><td>-9.06</td><td>8.8</td><td>0.19</td><td>39.7</td><td>1.61</td><td>0.04</td><td>135</td><td>688</td><td>4.66</td><td>18.8</td><td>26.9</td><td>M</td></tr></tbody></table>

<table><caption>**Table A4.** Characteristics of the simulations using the E1e-5.SF1 reference case with $E = 10^{-5}$, $Ra = 5 \times 10^{8}$, $Pm = 3$ and $Pr = 1$. The spatial discretization is set using $N_r^{OC} = 280$, $N_r^{IC} = 80$, $l_{max} = 213$ and $m_{max} = 170$. The behaviours are DS for dipolar stable, R for reversing, B for bistable and M for multipolar.</caption><thead><tr><th>Pattern</th><th>δq<sup>*</sup><sub>o</sub><br>(per cent)</th><th>Δt<sub>sim</sub><br>(t<sub>η</sub>)</th><th>∫<sub>dip</sub></th><th>δ̂(f<sub>dip</sub>)<br>(per cent)</th><th>M̄<br>(×10<sup>-4</sup>)</th><th>N<sub>cross</sub></th><th>Rm̄</th><th>Ro<sub>l</sub><br>(×10<sup>-2</sup>)</th><th>Ω<sub>f</sub><br>(per cent)</th><th>⟨U<sup>a</sup><sub>φ</sub>⟩<br>(per cent)</th><th>Behaviour</th></tr></thead><tbody><tr><td>–</td><td>0.00</td><td>5.1</td><td>0.61</td><td>7.6</td><td>32.12</td><td>19.92</td><td>0</td><td>399</td><td>0.76</td><td>-13.3</td><td>9.4</td><td>DS</td></tr><tr><td>Y20</td><td>5.80</td><td>3.1</td><td>0.53</td><td>9.8</td><td>21.24</td><td>12.7</td><td>0</td><td>444</td><td>0.83</td><td>-13.2</td><td>9.4</td><td>DS</td></tr><tr><td>Y20</td><td>14.50</td><td>1.8</td><td>0.33</td><td>8.7</td><td>18.03</td><td>10.7</td><td>0</td><td>489</td><td>0.82</td><td>-38.1</td><td>9.2</td><td>DS</td></tr><tr><td>Y20</td><td>-4.35</td><td>6.1</td><td>0.66</td><td>13.3</td><td>17.46</td><td>8.62</td><td>0</td><td>448</td><td>0.96</td><td>-9.5</td><td>8.6</td><td>DS</td></tr><tr><td>Y20</td><td>-8.70</td><td>9.1</td><td>0.33</td><td>40.5</td><td>7.28</td><td>1.25</td><td>15</td><td>525</td><td>1.27</td><td>-2.0</td><td>9.3</td><td>B</td></tr><tr><td>Y20</td><td>-14.50</td><td>2.4</td><td>0.18</td><td>43.8</td><td>6.09</td><td>0.48</td><td>17</td><td>558</td><td>1.38</td><td>5.2</td><td>11.5</td><td>M</td></tr></tbody></table>

<table><caption>**Table A5.** Characteristics of the simulations using the E1e-5.SF2 reference case with $E = 10^{-5}$, $Ra = 1.5 \times 10^{8}$, $Pm = 3$ and $Pr = 1$. The spatial discretization is set using $N_r^{OC} = 280$, $N_r^{IC} = 80$, $l_{max} = 213$, $m_{max} = 170$. The behaviours are DS for dipolar stable, R for reversing, B for bistable and M for multipolar.</caption><thead><tr><th>Pattern</th><th>δq<sup>*</sup><sub>o</sub><br>(per cent)</th><th>Δt<sub>sim</sub><br>(t<sub>η</sub>)</th><th>∫<sub>dip</sub></th><th>δ̂(f<sub>dip</sub>)<br>(per cent)</th><th>M̄<br>(×10<sup>-4</sup>)</th><th>N<sub>cross</sub></th><th>Rm̄</th><th>Ro<sub>l</sub><br>(×10<sup>-2</sup>)</th><th>Ω<sub>f</sub><br>(per cent)</th><th>⟨U<sup>a</sup><sub>φ</sub>⟩<br>(per cent)</th><th>Behaviour</th></tr></thead><tbody><tr><td>–</td><td>0.00</td><td>5.0</td><td>0.79</td><td>5.3</td><td>53.5</td><td>74.47</td><td>0</td><td>214</td><td>0.29</td><td>-17.5</td><td>7.4</td><td>DS</td></tr><tr><td>Y20</td><td>5.80</td><td>4.3</td><td>0.74</td><td>5.8</td><td>38.01</td><td>50.91</td><td>0</td><td>236</td><td>0.31</td><td>NC</td><td>7.6</td><td>DS</td></tr><tr><td>Y20</td><td>14.50</td><td>5.6</td><td>0.56</td><td>10.0</td><td>24.56</td><td>31.62</td><td>0</td><td>283</td><td>0.32</td><td>-31.4</td><td>7.6</td><td>DS</td></tr><tr><td>Y20</td><td>-2.90</td><td>7.1</td><td>0.89</td><td>3.1</td><td>34.81</td><td>69.47</td><td>0</td><td>238</td><td>0.34</td><td>NC</td><td>6.1</td><td>DS</td></tr><tr><td>Y20</td><td>-5.80</td><td>4.8</td><td>0.84</td><td>5.7</td><td>21.33</td><td>31.07</td><td>0</td><td>252</td><td>0.43</td><td>-10.5</td><td>6.0</td><td>DS</td></tr><tr><td>Y20</td><td>-8.70</td><td>10.2</td><td>0.54</td><td>35.1</td><td>11.74</td><td>3.34</td><td>3</td><td>279</td><td>0.55</td><td>-2.7</td><td>9.4</td><td>R</td></tr><tr><td>Y20</td><td>-10.15</td><td>15.5</td><td>0.65</td><td>23.9</td><td>12.96</td><td>7.24</td><td>1</td><td>278</td><td>0.53</td><td>-2.5</td><td>7.5</td><td>R</td></tr><tr><td>Y20</td><td>-11.02</td><td>9.3</td><td>0.54</td><td>31.7</td><td>11.55</td><td>2.88</td><td>4</td><td>287</td><td>0.58</td><td>0.2</td><td>10.2</td><td>R</td></tr><tr><td>Y20</td><td>-11.60</td><td>27.6</td><td>0.60</td><td>31.0</td><td>12.86</td><td>5.63</td><td>12</td><td>284</td><td>0.56</td><td>-0.5</td><td>8.1</td><td>R</td></tr><tr><td>Y20</td><td>-12.18</td><td>12.0</td><td>0.64</td><td>25.2</td><td>15.14</td><td>8.94</td><td>1</td><td>276</td><td>0.52</td><td>-1.0</td><td>7.5</td><td>R</td></tr><tr><td>Y20</td><td>-13.05</td><td>11.1</td><td>0.54</td><td>35.8</td><td>10.05</td><td>3.47</td><td>7</td><td>298</td><td>0.61</td><td>1.9</td><td>8.5</td><td>R</td></tr><tr><td>Y20</td><td>-14.50</td><td>13.9</td><td>0.49</td><td>42.9</td><td>10.63</td><td>3.36</td><td>17</td><td>304</td><td>0.64</td><td>2.7</td><td>9.0</td><td>R</td></tr><tr><td>Y20</td><td>-17.40</td><td>18.2</td><td>0.50</td><td>38.0</td><td>10.38</td><td>3.79</td><td>21</td><td>311</td><td>0.64</td><td>5.6</td><td>8.5</td><td>R</td></tr></tbody></table>
<figure><img src="image_16.png" /><figcaption>**Figure A1.** Time-series of the dipolar fraction and of the $M_{Z\alpha}^*$ parameter in a bistable dynamo. The dynamo is obtained using the E1e-4.hRm reference case and the $Y_{2,2}$ pattern with an amplitude $\delta q_o^* = 1.4$ per cent. The time-series are shown in blue when the dipolar fraction is higher than 0.25 and in orange when the dipolar fraction is lower than 0.25.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure A2.** Values of the magnetic Reynolds number (top) and of the Elsasser number (bottom) as a function of the heat flux heterogeneity amplitudes for the different patterns applied on the E1e-4.hRm reference case.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure A3.** Meridional cuts showing the azimuthal component of the magnetic field and of the velocity averaged in time and in the azimuthal direction relative to the root mean square of the fields for dynamo simulations using the ±Y₂.₀ pattern. See Fig. 9 for more details. The amplitudes of equatorial and polar cooling are given below for each reference dynamo model. E1e-4_hRm: δq₀* = ±7.2 per cent; E1e-5_hRm: δq₀* = ±14.5 per cent; E1e-6_hRm: δq₀* = ±9.1 per cent; E1e-5_SF1: δq₀* = ±14.5 per cent; E1e-5_SF2: δq₀* = ±14.5 per cent.</figcaption></figure>

## APPENDIX B: DYNAMO BISTABILITY

We classified six of our dynamo simulations as bistable. We recall that we define a dynamo as bistable when $f_{dip} > 0.25$ and $f_{dip} < 0.25$ for more than 20 per cent of the simulation time. Two of the simulations categorized as bistable have a relatively low dipolar fraction on average with large fluctuations which imply a significant amount of time spent in a multipolar state. The four other bistable cases, nonetheless, show two clear stable modes and the dynamo switches from one mode to another. In Fig. A1 are shown the time-series of the dipolar fraction and of $M_{Za}^*$ for one of the bistable case. This bistable dynamo is obtained using the E1e-4_hRm reference case and the Y₂.₂ pattern with an amplitude $δq_0^* = 1.4$ per cent (see Fig. 5). In this figure, the periods when the dynamo is dipolar ($f_{dip} > 0.25$) are shown in blue while the periods when the dynamo is multipolar ($f_{dip} < 0.25$) are shown in orange. Over the 37 diffusive time simulated, the dynamo switches four times between a dipolar and a multipolar behaviour. The dynamo stays for a long time in both modes, with up to ~13 $t_η$ (~520 kyr) spent continuously in a dipolar mode and ~5 $t_η$ (~200 kyr) spent continuously in a multipolar mode. Bistability has been previously documented in dynamo models (Simitev & Busse 2009; Busse & Simitev 2011; Schrinner et al. 2012; Petitdemange 2018) as well as in the context of the VKS dynamo experiment (Berhanu et al. 2009). However, this term usually refers to a hysteresis behaviour

in which the dynamo can reach different regimes depending on the initial conditions. Here, we obtain bistable dynamos that oscillate between two significantly different regime in the duration of one single run.

The bistable dynamos are obtained from originally dipolar simulations when destabilizing heat flux geometries are applied at the top of the core. Importantly, long simulation runs are necessary to identify such a behaviour, as bistable dynamos stay for a long time in either a dipolar or a multipolar regime. As a strongly dipolar magnetic field is used for the initial conditions of all the simulations, shorter simulation runs could have, for example, falsely led to the interpretation of these bistable cases as dipolar dynamos.

## APPENDIX C: EFFECT OF HEAT FLUX HETEROGENEITIES ON THE MAGNETIC FIELD AND FLOW AMPLITUDES

The evolutions of the magnetic Reynolds number and the Elsasser number as a function of the pattern amplitudes are shown in Fig. A2. In the main text, we saw that all the Y₁.₀ and -Y₂.₀ patterns tend to decrease the magnetic to kinetic energy ratio. As it can be seen in Fig. A2, this translates into both an increase of the strength of the flow, quantified by the magnetic Reynolds number, and a decrease of the magnetic field amplitude, quantified by the Elsasser
number. The −*Y*₂,₀ pattern is the only pattern to induce significant variations in the magnetic Reynolds number and the Elsasser number for the lowest amplitude. The +*Y*₂,₀ has an opposite effect, by decreasing *Rm* and increasing Λ. The opposite effects of the +*Y*₂,₀ pattern and the −*Y*₂,₀ pattern is particularly striking. As shown in Fig. 6, the +*Y*₂,₀ pattern preserves the dipole-dominated structure of the magnetic field, even for large amplitudes, and slightly increases the magnetic energy relatively to the kinetic energy. On the contrary, the −*Y*₂,₀ significantly decreases the energy ratio and trig-

gers a transition towards a multipolar behaviour even for the lowest heterogeneity amplitudes.

APPENDIX D: EFFECT OF THE Y
PATTERN ON THE MEAN FLOW AND
MAGNETIC FIELD

We show in Fig. A3 an extension of Fig. 9 for the five dynamo
models.