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

HAL
open science

# Dynamo Models With a Mercury-Like Magnetic Offset Dipole

P. Kolhey, D. Heyner, J. Wicht, T. Gastine, K. -H. Glassmeier, F. Plaschke

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

To cite this version:

P. Kolhey, D. Heyner, J. Wicht, T. Gastine, K. -H. Glassmeier, et al., Dynamo Models With a Mercury-Like Magnetic Offset Dipole. Journal of Geophysical Research. Planets, 2025, 130, pp.2024JE008660. 10.1029/2024JE008660. insu-05080919

HAL Id: insu-05080919

https://insu.hal.science/insu-05080919v1

Submitted on 23 May 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Planets

## RESEARCH ARTICLE
10.1029/2024JE008660

**Key Points:**

* We analyze a double-diffusive and single-diffusive numerical model for Mercury's dynamo, both including a thick thermally stratified layer

* Both models are able to explain Mercury's field characteristics, when dynamo action sets in with a fundamental equatorially symmetric mode

* The double-diffusive model needs to operate close to the onset of convection to generate Mercury-like fields, which seems unrealistic

**Correspondence to:**

P. Kolhey, p.kolhey@tu-braunschweig.de

**Citation:**

Kolhey, P., Heyner, D., Wicht, J., Gastine, T., Glassmeier, K.-H., & Plaschke, F. (2025). Dynamo models with a Mercury-like magnetic offset dipole. *Journal of Geophysical Research: Planets*, 130, e2024JE008660. https://doi.org/10.1029/2024JE008660

Received 29 JUL 2024
Accepted 3 FEB 2025

**Author Contributions:**

**Conceptualization:** P. Kolhey, D. Heyner, J. Wicht, T. Gastine, K.-H. Glassmeier

**Data curation:** P. Kolhey, J. Wicht, T. Gastine

**Formal analysis:** P. Kolhey

**Funding acquisition:** F. Plaschke

**Investigation:** P. Kolhey

**Project administration:** D. Heyner, J. Wicht

**Software:** P. Kolhey, J. Wicht, T. Gastine

**Supervision:** D. Heyner, J. Wicht, K.-H. Glassmeier, F. Plaschke

**Validation:** P. Kolhey, J. Wicht, T. Gastine

**Visualization:** P. Kolhey

**Writing – original draft:** P. Kolhey

**Writing – review & editing:** D. Heyner, J. Wicht, T. Gastine, K.-H. Glassmeier

© 2025. The Author(s).

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

## Dynamo Models With a Mercury-Like Magnetic Offset Dipole

P. Kolhey¹, D. Heyner¹, J. Wicht², T. Gastine³, K.-H. Glassmeier¹,², and F. Plaschke¹

¹Institut für Geophysik und Extraterrestrische Physik, Technische Universität Braunschweig, Braunschweig, Germany,

²Max Planck Institute for Solar System Research, Göttingen, Germany, ³Université Paris Cité, Institut de Physique du Globe de Paris, UMR 7154 CNRS, Paris, France

**Abstract** Mercury has a global magnetic field that is unique in our solar system. It is dominated by a weak axial dipole and a substantial axial quadrupole contribution. The field has a dipole tilt which is smaller than one degree and an offset of the magnetic equator toward north by about 20% of the planetary radius. The only dynamo model that succeeds in continuously reproducing these features without having to rely on an unrealistic heat-flux pattern through the core-mantle boundary is a double-diffusive model by Takahashi et al. (2019), https://doi.org/10.1038/s41467-018-08213-7, where thermal effects cause a thick, stably stratified layer in the outer parts of the core. In this study we further show the sensitivity of this model to the choice of parameters. In the explored parameter set, the model yields solutions which are unrealistic for Mercury when convection becomes vigorous. We present a new single-diffusive model, that captures Mercury's field characteristics during long periods that are interrupted by dipole field reversals causing the magnetic equator to switch between the northern and southern hemisphere. The Mercury-like solutions are stable over a broader range of parameters. An important ingredient for the success is the fact that dynamo action sets in with the equatorially symmetric magnetic field family.

**Plain Language Summary** Mercury's magnetic field is quite unique when compared to other planetary magnetic fields in our solar system. It is very weak and has a very simple dipole structure. However, the magnetic equator is shifted northward. Numerical models which try to explain the origin of the magnetic field from the dynamics of the electrically conducting fluid outer core struggle in explaining the latter feature. So far only one model has succeeded in continuously reproducing these properties. The model includes a thick thermally stably stratified region, where the convective motions in the core are strongly suppressed and convection below this region is driven by differences in temperature and chemical composition. In this study we explore this model further, showing its sensitivity to the choice of model parameters. By exploring the influence of parameter changes on dynamo action, we show that this model yields unrealistic solutions for Mercury when convection becomes vigorous. We also present a new model where core convection is solely driven by differences in temperature, which also captures all properties of Mercury's magnetic field, but in a wider range of parameters.

## 1. Introduction

Mercury's interior structure and dynamics remain enigmatic. However, the peculiar magnetic field produced in the planet's iron core provides important constraints on the internal structure of the planet. After the orbital MESSENGER (MErcury Surface, Space ENvironment, GEochemistry, and Ranging) mission (Solomon et al., 2018) it seems that the exterior far field is best described by a dipole field which is strictly aligned along the planet's rotation axis but significantly shifted northward by about 20% of Mercury's radius. Various detailed examination of the magnetic field structure by Anderson et al. (2012), Thébault et al. (2018), Wardinski et al. (2019), Toepfer et al. (2021), Wardinski et al. (2021), and Pump et al. (2024) (see Genova et al. (2021) for a review) reveal that the dipole field strength is about -215 to -190 nT, which is rather weak compared to Earth's magnetic dipole field strength of around -29405 nT (Alken et al., 2021).

Planetary dynamos are expected to operate in a regime where the Coriolis force and the Lorentz force are of comparable magnitude in the Navier-Stokes equation. Their ratio can be characterized by the Elsasser number

$$ \Lambda = \frac{B^2 \sigma}{\rho \Omega}, \qquad (1) $$
<table><caption>**Table 1**
Summary of Dynamo Models to Explain Mercury's Magnetic Field Properties From Various Studies</caption><thead><tr><th colspan="3">Model setup</th><th colspan="5">Magnetic field properties</th></tr><tr><th>Convective driving</th><th>Inner core size</th><th>Stable layer? (Thickness)</th><th>Therm. Boundary condition at CMB</th><th>Imposed external magnetic field?</th><th>Weak magnetic field?</th><th>Axisymmetric field?</th><th>Persistent offset dipole field?</th><th>Reference</th></tr></thead><tbody><tr><td>Single-diffusive/bottom driving</td><td>a = 0.35–0.60</td><td>Yes (d<sub>s</sub> = 0.56–0.7)</td><td>Homogeneous (isothermal)</td><td>No</td><td>Yes</td><td>Yes</td><td>No</td><td>Christensen (2006), Christensen and Wicht (2008)</td></tr><tr><td>Double-diffusive/bottom driving</td><td>a = 0.35</td><td>Yes (d<sub>s</sub> = 0.36–0.44)</td><td>Homogeneous (fixed heat flux)</td><td>No</td><td>No</td><td>No</td><td>No</td><td>Manglik et al. (2010)</td></tr><tr><td>Single-diffusive/bottom driving</td><td>a = 0.34</td><td>Yes (double snow layer)</td><td>Homogeneous (fixed heat flux)</td><td>No</td><td>Yes</td><td>Yes</td><td>No</td><td>Vilim et al. (2010)</td></tr><tr><td>Single-diffusive/bottom driving</td><td>a = 0.35</td><td>No</td><td>Homogeneous (isothermal)</td><td>Yes</td><td>Yes</td><td>Yes</td><td>No</td><td>Glassmeier et al. (2007), Heyner, Wicht, et al. (2011)</td></tr><tr><td>Single-diffusive/volumetric driving</td><td>a = 0.2</td><td>No</td><td>SH degree-2 heat flux pattern</td><td>No</td><td>No</td><td>Yes</td><td>Yes</td><td>Cao et al. (2014)</td></tr><tr><td>Single-diffusive/bottom driving</td><td>a = 0.05–0.50</td><td>Yes (d<sub>s</sub> = 0.31–0.54)</td><td>SH degree-1 heat flux pattern</td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Tian et al. (2015)</td></tr><tr><td>Double-diffusive/therm. volumetric driving/comp. bottom driving</td><td>a = 0.2</td><td>Yes (d<sub>s</sub> = 0.625)</td><td>Homogeneous (fixed heat flux)</td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Takahashi et al. (2019), this study</td></tr><tr><td>Single-diffusive/bottom driving</td><td>a = 0.2</td><td>Yes (d<sub>s</sub> = 0.625)</td><td>Homogeneous (fixed heat flux)</td><td>No</td><td>Yes</td><td>Yes</td><td>Yes</td><td>This study</td></tr></tbody></table>

Note. The inner core size is expressed by the ratio $a = r_i/r_o$, where $r_i$ and $r_o$ are the inner core radius and outer core radius, respectively. The stable layer thickness is given by $d_s = (r_o - r_i)/(r_o - r_i)$, where $r_i$ is the radius of the bottom of the stable layer.

where B is the magnetic field strength, σ the electrical conductivity, ρ the mass density, and Ω the planetary rotation rate. While we expect it to be $Λ ≈ 1$ it is only about $Λ ≈ 10^{-5}$ at Mercury's surface (Wicht & Heyner, 2014).

Over the past decades, several dynamo studies proposed different model solutions in order to explain the observed magnetic field strength and structure. Table 1 therefore provides an overview of important model setups and their resulting field characteristics, which are explained in detail below.

Owing to its close proximity to the Sun, Mercury is exposed to a strong solar wind which interacts with the planet's internally generated magnetic field, leading to a small magnetosphere. At the magnetopause, the boundary of this interaction region, so-called Chapman-Ferraro currents can generate a secondary magnetic field oppositely aligned to the dynamo field at the polar regions. As pointed out by Glassmeier et al. (2007), it is plausible that the long-period, large-scale components of the external magnetic field diffuse down to the core region where they would interact with the internal field produced by the dynamo. This feedback mechanism might account for Mercury's weak magnetic field as demonstrated by simplified mean-field dynamo models (Heyner et al., 2010; Heyner, Schmitt, et al., 2011). Depending on the dipole orientation and field strength, fully nonlinear dynamo simulations showed that a homogeneous external magnetic field can weaken the internal field generated by the dynamo (Gómez-Pérez & Solomon, 2010; Gómez-Pérez & Wicht, 2010). Furthermore, Heyner, Wicht et al. (2011) have been able to demonstrate that an external field depending on the internal dynamo field strength is able to maintain a weak internal magnetic field if the initially imposed internal field is also weak. For instance, the onset of dynamo action or polarity reversals of the dipole field would allow for a weak internal field strength.
Alternative explanations for Mercury's magnetic field arise from the planet's unique internal structure and dynamics. In the context of an external weakening magnetic field, Vilim et al. (2010) proposed a dynamo model wherein two convection zones in the outer core are separated by a stably stratified iron snow region, leading to two independent dynamo regions that generate magnetic fields with mutually negative effects.

Given the uncertain size of Mercury's solid inner core, early dynamo models of the pre-MESSENGER era assumed both, a very small solid inner core (Heimpel et al., 2005) or a very large inner core (Stanley et al., 2005; Takahashi & Matsushima, 2006). These models were able to produce a weak poloidal magnetic field which is observed outside the fluid core. However, none of the above mentioned models reproduce the field morphology observed by MESSENGER. Especially the large quadrupole component, which causes the offset of the magnetic equator in the far field, is challenging to explain.

In order to achieve such a dipole offset, Cao et al. (2014) introduced a heat flux pattern of spherical harmonic degree 2, allowing more heat to escape at low latitudes than in higher latitudes. In combination with a volumetric buoyancy source, the resulting convection was able to sustain a persistent dipole offset in the magnetic field. Although this dynamo model offers a possible explanation for Mercury's unusual magnetic field morphology, it does not account for the weak field strength. Additionally, the origin of the heat flux pattern necessary for this solution remains dubious.

Already before MESSENGER, a stably stratified layer underneath the core-mantle boundary (CMB) was suggested, since the temperature gradient at the CMB could be sub-adiabatic (see Christensen, 2006; Genova et al., 2021, and reference therein). Hence, Christensen (2006) and Christensen and Wicht (2008) were the first who performed dynamo simulations for Mercury that included a stably stratified layer in the outer part of the iron core. This layer acts as an electromagnetic low-pass filter, damping small-scale, fast-varying magnetic field components more efficiently due to a magnetic skin effect. In consequence, a weak and steady dipole-dominated field remains at the CMB. Differential rotation in the stable layer can further attenuate non-axisymmetric field components. While the dipole tilt and the quadrupole-to-dipole ratio of this dynamo model matches Mercury values at some instances in time, these properties exhibit significant temporal variations and also unrealistically larger non-axisymmetric field contributions, as later analysis revealed (Wicht & Heyner, 2014).

Tian et al. (2015) demonstrated that it is possible to produce a Mercury-like magnetic field when the model includes a stably stratified layer and a degree-1 heat flux pattern at the CMB where more heat escapes in the northern hemisphere. Just as for the study by Cao et al. (2014), the assumed CMB heat flux patterns of both studies seem unlikely (e.g., Tosi et al., 2013, 2015). Mantle convection, if still active, would likely consist of small scale convection cells that do not give rise to large-scale temperature variations at the CMB. If the mantle is in a conductive state (no convection anymore), long-term insolation patterns of degree and order 2 or long-wavelength variations in topography (Smith et al., 2012; Zuber et al., 2012) could induce temperature variations at the CMB but the resulting heat flux would be opposite to the imposed heat flux proposed by Cao et al. (2014). Furthermore it does not create a north-south asymmetry as considered by Tian et al. (2015).

Manglik et al. (2010) suggest a double-diffusive approach, treating thermal and compositional convection separately. While thermal effects would yield a thick stratified outer layer like in the models by Christensen (2006) and Christensen and Wicht (2008), already a mild contribution of compositional convection is sufficient to destroy this layer. Thus the resulting magnetic field is too strong for Mercury. In contrast, the recent model by Takahashi et al. (2019) uses a double-diffusive approach that not only maintains a thick enough thermally stable outer layer but also shows a persistent dipole offset without assuming any exceptional heat flux patterns at the CMB. Different to Manglik et al. (2010) they assumed a smaller solid inner core and a stronger thermal stratification. Interestingly, the magnetic field, and thereby also the Lorentz force, appears to play a crucial role in the symmetry-breaking of the convective motions. While a simulation omitting the action of Lorentz force only generates an axial quadrupolar field, self-consistent dynamo simulations reproduce a Mercury-like offset dipole field. A detailed explanation of the underlying mechanism how the symmetry-breaking is achieved is still missing, but the results from Takahashi et al. (2019) suggest that the magnetic field plays a key role.

Based on the results of currently debated dynamo simulations and thermal evolution models, it seems reasonable to assume that the upper part of Mercury's core is stably stratified. This allows to explain the weakness and axisymmetry of the magnetic field. Different thicknesses of the stable layer have been suggested for the various
dynamo models and the inner core size varies among these, too. While up to about two-thirds of the remaining fluid layer depth are assumed to be stably stratified, the inner core radii vary from 5% to 60% of the total core. Recent geodetic models favor a large inner core with 50% of the total core size (Genova et al., 2019), however smaller inner cores are not fully excluded (Knibbe et al., 2021; Steinbrügge et al., 2021).

Another critical aspect in the dynamo models is the convective driving mechanism, which is crucially tied to the chemical composition of the planetary core. To keep Mercury's core partially liquid, relative large contributions of lighter elements are required. One potential candidate is sulfur (Hauck et al., 2013; Nittler et al., 2018). In sub-eutectic compositions, sulfur does not partition into the solid phase. It is thus released at the growing inner-core front and hence drives compositional convection. However, large amounts of sulfur found at Mercury's surface suggest that the mantle has solidified under reducing conditions. Little to no sulfur would then partition into the core. However, it seems also possible that the detected high sulfur concentrations (Nittler et al., 2018) only represent very local concentrations, opening the possibility for more sulfur in the core (Renggli et al., 2022). For the other light element candidate, silicon (Davies et al., 2024; Edmund et al., 2022; Knibbe & van Westrenen, 2018), this does not seem possible since it almost equally partitions into solid and liquid phase. However, an Fe-Si core may release a sufficient amount of latent heat during core crystallization to power the dynamo (Davies et al., 2024; Knibbe & van Westrenen, 2018). Therefore, if a substantial amount of silicon is present in the core, latent heat release at the inner core boundary (ICB) cannot be neglected, as was done in previous studies (Cao et al., 2014; Takahashi et al., 2019). However, diamond-anvil-cell experiments by Edmund et al. (2022) suggest that depending on the silicon concentration a lighter solid Fe-Si phase may form and rise to the surface as a form of compositional convection. Although volumetric heating by radioactive elements is a plausible driver for convection, given the surface abundance of potassium, thorium, and uranium (Peplowski et al., 2011), abundances in the core are likely too small to significantly contribute to powering convection.

In this study, we first revisit the double-diffusive dynamo model by Takahashi et al. (2019) and systematically explore several open questions. A parameter study reveals how typical the published solution is. We then carefully analyze the reason for the dipole offset. Finally, we present simulations of a new purely thermal (co-density) model including a stably stratified layer at the top of the core. This differs from the model used by Christensen (2006) and Christensen and Wicht (2008) in terms of inner core size, as well as in the implementation of the stably stratified layer. Our new model is able to reproduce all of Mercury's magnetic field characteristics very well.

The paper is structured as follows: In Section 2 we introduce the mathematical and numerical foundation of the dynamo models presented in this paper. Afterward we present our results of the simulations in Section 3, beginning with the double-diffusive model in Section 3.1, which is then contrasted to the results of the single-diffusive model (purely thermal/co-density) in Section 3.2. Finally, we discuss the outcomes in Section 4 before drawing our conclusions in Section 5.

## 2. The Dynamo Models

### 2.1. Model Equations

In this study we consider an electrically conducting fluid with magnetic diffusivity $\lambda$ in a rotating spherical shell with an inner core radius $r_i$ and an outer core radius $r_o$. The outer and inner boundary represent the CMB and solid ICB, respectively. The system rotates with a constant rotation frequency $\Omega$ about the vertical z-axis. We use a dimensionless form of the magnetohydrodynamic set of equations in the Boussinesq limit. We assume that convection is driven by thermal and by compositional gradients. Thus, density variations are given by the simplified equation of state

$$ \rho = \rho_0 [1 - \alpha(T - T_0) - \beta(\xi - \xi_0)], \quad (2) $$

where $T$ is the temperature, $\xi$ the concentration of light elements, $\rho_0$, $T_0$ and $\xi_0$ denote the reference states of the density, temperature and chemical composition in the core, respectively. Furthermore, $\alpha$ is the thermal expansion coefficient, whereas $\beta$ is an analogue quantity for the chemical composition. We use the shell thickness $d = r_o - r_i$ as a length scale, the magnetic diffusion time $t_\lambda = d^2/\lambda$ as a time scale and consequently $\lambda/d$ as a velocity scale. The magnetic field is scaled by $\sqrt{\rho_0\mu_0\lambda\Omega}$ with $\mu_0$ denoting the magnetic permeability in vacuum.
Depending on the boundary conditions used in the model, temperature is scaled by either the temperature gradient at the outer boundary $F_{T,o} = |\partial T/\partial r|_{r_o}$ or the temperature gradient at the inner boundary $F_{T,i} = |\partial T/\partial r|_{r_i}$, both multiplied with the length scale $d$. The chemical composition is scaled by the compositional gradient at the ICB $F_{\xi,i} = |\partial \xi/\partial r|_{r_i}$ multiplied with the length scale $d$. The dynamo process is then described by the conservation of mass, the condition of a divergence-free magnetic field, the Navier-Stokes equation, the induction equation and the transport equations for temperature, as well as the chemical composition:

$$ \nabla \cdot \mathbf{u} = 0, \quad \nabla \cdot \mathbf{B} = 0, \qquad (3) $$

$$ \frac{E}{Pm} \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) = -\nabla p - 2 \mathbf{e}_z \times \mathbf{u} + \frac{Ra_\theta PmE}{Pr} \frac{r}{r_o} \mathbf{\theta} \mathbf{e}_r + \frac{Ra_\xi PmE}{Sc} \frac{r}{r_o} \mathbf{\xi} \mathbf{e}_r + (\nabla \times \mathbf{B}) \times \mathbf{B} + E \nabla^2 \mathbf{u}, \qquad (4) $$

$$ \frac{\partial \mathbf{B}}{\partial t} = \nabla \times (\mathbf{u} \times \mathbf{B}) + \nabla^2 \mathbf{B}, \qquad (5) $$

$$ \frac{\partial \theta}{\partial t} + (\mathbf{u} \cdot \nabla)\theta + u_r \frac{dT_0}{dr} = \frac{Pm}{Pr} \nabla^2 \theta, \qquad (6) $$

$$ \frac{\partial \xi}{\partial t} + (\mathbf{u} \cdot \nabla)\xi = \frac{Pm}{Sc} \nabla^2 \xi + \epsilon_{\xi}, \qquad (7) $$

Where $\mathbf{u}$ is the velocity field, $\mathbf{B}$ the magnetic field, $p$ the non-hydrostatic pressure; $\mathbf{e}_r$ denotes the unit vector in the radial direction, $\theta = T - T_o$ the temperature fluctuation and $\xi$ the chemical composition, respectively. $\epsilon_\xi$ describes a uniform compositional source/sink and gravity is assumed to vary linearly with radius. The system described by Equations 3–7 is controlled by six non-dimensional parameters which are the Ekman number $E$, the thermal Rayleigh number $Ra_\theta$, the compositional Rayleigh number $Ra_\xi$, the Prandtl number $Pr$, the magnetic Prandtl number $Pm$ and the Schmidt number $Sc$

$$ \begin{aligned} E &= \frac{\nu}{\Omega d^2}, \quad Ra_\theta = \frac{\alpha g_o F_T d^4}{\nu \kappa_\theta}, \quad Ra_\xi = \frac{\beta g_o F_{\xi,i} d^4}{\nu \kappa_\xi}, \\ Pr &= \frac{\nu}{\kappa_\theta}, \quad Pm = \frac{\nu}{\lambda}, \quad Sc = \frac{\nu}{\kappa_\xi} \end{aligned} \qquad (8) $$

where $\nu$ is the kinematic viscosity and $g_o$ is the gravitational acceleration at the outer boundary. The Schmidt number $Sc$ is equivalent to the Prandtl number, replacing the thermal diffusivity $\kappa_\theta$ by the compositional diffusivity $\kappa_\xi$. In planetary interiors, the compositional diffusivity is much smaller than the thermal: $Pr \ll Sc$.

## 2.2. The Double-Diffusive Model

Following Takahashi et al. (2019) we assume mechanical no-slip boundary conditions, which corresponds to $\mathbf{u} = 0$ at $r \in \{r_i; r_o\}$. Furthermore, the region above the CMB ($r > r_o$) is assumed to be electrically insulating, such that the magnetic field in this region can be described by a potential field and the field inside the fluid shell has to match a potential field at $r = r_o$. At the inner boundary, the region $r < r_i$ is assumed to be electrically insulating as well, which is a strong simplification motivated by the small inner core size. A fixed compositional flux is assumed through the ICB, whereas it is zero at the CMB. For the temperature boundary conditions it is assumed that the heat flux at the ICB is zero, while the heat flux through the CMB has non-zero value. Hence, temperature is scaled by the temperature gradient at the CMB $F_{T,o}$ without stably stratified layer, which defines $Ra_\theta$. We note that assuming a vanishing thermal gradient at the lower boundary seems unrealistic, since latent heat and secular cooling of the inner core will certainly yield a significant heat flux. Following Takahashi et al. (2019) a thermally stratified region at the top of the fluid core is included by adding a sub-adiabatic region to the conductive solution of the background temperature gradient, corresponding to
<figure><img src="image_5.png" /><figcaption>**Figure 1.** (a) Temperature and compositional background gradient profiles for the double-diffusive reference case DD2 and the purely thermal reference case TH11. The thermal background state includes a stably stratified region defined by Equation 9 with $\Gamma_{T,o} = 10$, $\tau_s = 0.05 r_o$ and $r_s = 0.5 r_o$ in the double-diffusive case, while the thermal background state for the thermal model is defined by Equation 12 with $\Gamma_{T,o} = 3$, $\tau_s = 0.05 r_o$ and $r_s = 0.5 r_o$. The vertical dotted line marks the radius $r_s$. (b) Resulting ratio between the Brunt-Väisälä frequency and planetary rotation frequency defined in Equation 10 for the double-diffusive reference case DD2 and the single-diffusive reference case TH11.</figcaption></figure>

$$ \frac{dT_0}{dr} = \Gamma_{T,o} - \frac{1}{2} \left\{ 1 - \tanh \left( \frac{r-r_s}{\tau_s} \right) \right\} \left\{ \Gamma_{T,o} + \frac{r^3 - r_i^3}{r_o^3 - r_i^3} \left( \frac{r_o}{r} \right)^2 \right\}, \quad (9) $$

where $\Gamma_{T,o}$ is the amplitude of the temperature gradient at the outer boundary, $r_s$ the radius of the bottom of the stably stratified layer and $\tau_s$ is the thickness of the transition zone between the thermally stably stratified region and the convective region, respectively. In Figure 1 the background temperature gradient is shown for the reference case DD2, which will be presented below. We define the convective part as the region where the combined buoyancy power is positive. The thermal component is driving in the lower part of this region but stabilizing in the upper part.

In rapidly rotating systems a measure for the stability of the stratified layer is the ratio $\tilde{N}$ between the Brunt-Väisälä frequency $N$ and the planetary rotation frequency $\Omega$ (Monville et al., 2019):

$$ \tilde{N}^2 = \frac{N^2}{\Omega^2} = \frac{Ra_\theta E^2}{Pr} \frac{r}{r_o} \frac{dT_0}{dr} + \frac{Ra_\xi E^2}{Sc} \frac{r}{r_o} \frac{d\xi_0}{dr}. \quad (10) $$

Since the maximum degree of stratification is reached at the CMB (cf. Figure 1), where the destabilizing effect of the compositional component is negligible, we use

$$ \tilde{N}_o = \frac{N_o}{\Omega} \approx \sqrt{\frac{Ra_\theta E^2}{Pr} \Gamma_{T,o}}, \quad (11) $$

as a measure for the degree of stratification with $N_o$ as the Brunt-Väisälä frequency at the CMB.

Our simulation run DD2 (short for double-diffusive and called BU1 in the study by Takahashi et al. (2019)) will serve as reference for this study, because it produces a Mercury-like magnetic field structure (see Table A1). Similarly to Takahashi et al. (2019), we set the Prandtl number to $Pr = 0.1$ and the Schmidt number is $Sc = 1.0$,
corresponding to a Lewis number of $Le = Sc/Pr = 10$. Higher diffusivity contrasts are numerically demanding and will not be investigated in this study. The choice of the Ekman number is very moderate with $E = 2.0 \times 10^{-4}$ in the reference model, but a few simulations are performed at $E = 10^{-4}$. The thermal Rayleigh number is set to $Ra = 5.0 \times 10^{6}$, whereas the compositional Rayleigh number is $Ra_{s} = 10^{8}$, but we will vary both Rayleigh numbers systematically to investigate the systems behavior on the convective state of the core. The magnetic Prandtl number is $Pm = 3$ for the reference model and will be varied in the range $2 \le Pm \le 5$ in this study. With the moderate choice of Ekman number, dynamo action is not achieved for $Pm < 2$. The aspect ratio is fixed to $a = r_{i}/r_{o} = 0.2$. For the thermally stratified layer defined in Equation 9 the amplitude of the background temperature gradient at the CMB is $\Gamma_{T,o} = 10$ for the reference model and the thickness of the transition zone between the convective and stably stratified region is $\tau_{s} = 0.05 r_{o}$. If necessary, the amplitude $\Gamma_{T,o}$ will be changed to keep the degree of stratification of the stable layer $\tilde{N}_{o}$ in the same range as the reference case, while the thickness of the transition zone will be held fixed. The bottom radius of the stably stratified layer is set to $r_{s} = 0.5 r_{o}$.

## 2.3. The Single-Diffusive Model

In our purely thermally driven model (single-diffusive) we neglect the convective driving by the release of light elements such that Equation 7 does not need to be solved. Similar to the double-diffusive model we assume mechanical no-slip boundary conditions and electrically insulating regions above the CMB and the below the ICB. The heat fluxes at the boundaries are assumed to be constant corresponding to the background state of the system. In this case we use the temperature gradient at the ICB $\tilde{F}_{T,i}$ to scale the temperature and to define $Ra_{\theta}$. With a constant heat flux at the inner core and a stratified layer at the top of the core the background profile is given by

$$ \frac{dT_0}{dr} = \Gamma_{T,o} - \frac{1}{2} \left\{ 1 - \tanh \left( \frac{r - r_s}{\tau_s} \right) \right\} (\Gamma_{T,o} + 1), \quad (12) $$

which is shown in Figure 1 for the reference case TH11 defined below. Temperature and/or compositional profiles in the potential stable layer underneath Mercury's CMB remain unclear. We chose simple profiles that succeed in suppressing radial motions.

The simulation run TH11 (short for purely thermal) will serve as a reference case (see Table A2). In order to allow for a comparison to the double-diffusive model, we choose $Pr = 1.0$, $E = 1 \times 10^{-4} - 2 \times 10^{-4}$ and we explore magnetic Prandtl numbers ranging from $Pm = 2$ to $Pm = 5$. The Rayleigh number is varied between $14 Ra_{c} \le Ra \le 228 Ra_{c}$, where $Ra_{c} \approx 3.5 \times 10^{6}$ for $E = 10^{-4}$ and $Ra_{c} \approx 1.8 \times 10^{6}$ for $E = 2 \times 10^{-4}$. Since convection will essentially develop in the lower part of the fluid domain, the Rayleigh number defined on the shell gap $d$ actually overestimates the effective Rayleigh number of the convective part alone, explaining the highly overcritical values needed for $Ra$ to ensure dynamo action. Again, we fix the aspect ratio to $a = 0.2$. For the stably stratified layer we choose $\tau_{s} = 0.05 r_{o}$ and $r_{s} = 0.5 r_{o}$, similar to Takahashi et al. (2019), while $\Gamma_{T,o}$ will be varied to keep the degree of stratification $\tilde{N}_{o}$ in the same range as the reference case TH11 with $\Gamma_{T,o} = 3$ presented below. Note that the lower value of $\Gamma_{T,o}$ for the thermal model is a consequence of the different background state (bottom-up driving) and the single-diffusive approach, where convection is now solely driven by the temperature perturbations but simultaneously stratified in the upper part of the core. Hence, higher thermal Rayleigh numbers are necessary to drive convection compared to the double-diffusive case, such that the amplitude of the thermal background gradient at the CMB can be reduced to obtain a similar degree of stratification (see Equation 11).

## 2.4. Numerical Method

Equations 3–7 are solved numerically using the well-established pseudo-spectral code MagIC (freely available under https://github.com/magic-sph/magic, see Wicht, 2002). The code makes use of the property that the solenoidal velocity and magnetic field can be separated into their scalar poloidal and toroidal potentials,

$$ u = \nabla \times (\nabla \times W e_r) + \nabla \times Z e_r, \quad (13) $$
$$B = \nabla \times (\nabla \times G e_r) + \nabla \times H e_r. \quad (14)$$

The remaining variables in the transport equations are expanded in spherical harmonic functions in $(\theta, \phi)$-direction up to degree $l_{max}$ and the transformations between spectral and grid representations are carried out using the SHTns library (Schaeffer, 2013). In the radial direction, the variables are expanded in Chebyshev polynomials up to degree $N_r$.

Typically, Chebyshev polynomials are defined on a Gauss-Lobatto grid, which has a dense spacing of nodal points at the boundaries with the advantage that boundary layers are resolved properly. In dynamo simulations, however, this can cause severe time step limitations due to the propagation of Alfvén waves near the boundaries when the magnetic field is of noticeable strength (Christensen et al., 1999). Therefore, we use a mapping function defined by Kosloff and Tal-Ezer (1993) to stretch the grid to be more uniformly spaced and avoid the severe time step limitations (see Gastine et al., 2020, for more details). Additional time step limitations can arise by introducing a strongly stratified thermal layer since gravity waves can be excited. Hence, the advection term $u_r$ $dT_0/dr$ is treated implicitly (Gastine et al., 2020). An implicit-explicit third-order Runge-Kutta scheme (BPR353) (Boscarino et al., 2013) is used for time integration, where the Coriolis and nonlinear terms are treated explicitly (see Gopinath et al., 2022, for more details on the time stepping scheme). Typically, we start a dynamo simulation from a purely hydrodynamic simulation run, where we include a randomly perturbed seed magnetic field. Some simulations, where we change one of the control parameters, however, are started from a fully evolved dynamo case.

## 2.5. Simulation Diagnostics

In our analysis of the simulations we use different output diagnostics. For the diagnostic parameters we use overbars for a temporal average and angle brackets for spatial average of an arbitrary variable $f$:

$$\bar{f}(\mathbf{r}) = \frac{1}{t_a} \int_{t_0}^{t_0+t_a} f(\mathbf{r}, t) \, dt, \quad \langle f \rangle_V(t) = \frac{1}{V} \int_V f(\mathbf{r}, t) \, dV,$$

where $t_a$ is the total averaging time and $V$ the volume of the part of the core we are interested in. Note that the diagnostic quantities will also be evaluated in different radial zones of the fluid region, mainly in the convection zone. Hence, the volume integral will be evaluated accordingly.

As a proxy for the flow amplitude we use the Reynolds number and magnetic Reynolds number which are defined as

$$\mathrm{Re} = \overline{\langle \mathbf{u}^2 \rangle_V^{1/2}} / Pm, \quad \mathrm{Rm} = \overline{\langle \mathbf{u}^2 \rangle_V^{1/2}}. \quad (15)$$

The kinetic energy $E_k$ and magnetic energy $E_m$ are given by

$$E_k + E_m = \frac{1}{2} \overline{\langle \mathbf{u}^2 \rangle_V} + \frac{Pm}{2E} \overline{\langle \mathbf{B}^2 \rangle_V}. \quad (16)$$

In the following we will analyze different contributions of the kinetic and magnetic energy, which are axisymmetric (ax), equatorial symmetric (es) and equatorial anti-symmetric (eas).

The mean magnetic field amplitude is measured by the Elsasser number

$$\Lambda = \overline{\langle \mathbf{B}^2 \rangle_V}. \quad (17)$$

In order to quantify the dipole offset we assume that such a magnetic field configuration can be described by an internal axial dipole which is shifted by an equatorial offset Z. It is possible to expand the field of an offset axial dipole by an infinite series of centered axisymmetric Gauss coefficients that obey the relation
$$g_{\ell}^{0} = \ell g_{1}^{0} Z^{\ell-1}, \qquad (18)$$

where $\ell$ is the spherical harmonic degree and $Z = Z/R_M$ is the relative offset (Alexeev et al., 2010; Wicht & Heyner, 2014), with $R_M$ as the radius of Mercury's surface. Further away from the planet only dipole and quadrupole contributions remain sizable and the relative offset simply becomes

$$Z \approx \frac{g_2^0}{2g_1^0}. \qquad (19)$$

We adopt this simplified definition in the following, keeping in mind that the offset may actually change closer to the planet unless relation (18) is also fulfilled for further harmonics. See Appendix B for the calculation of the Gauss coefficients from the poloidal magnetic field coefficients of the dynamo model. Also the dipole tilt $\Delta_{dip}$ will be investigated, which is defined as the smallest angle between rotation axis and dipole axis.

# 3. Results

## 3.1. The Double-Diffusive Model

Figure 2a shows how the equatorially symmetric and anti-symmetric contributions of kinetic and magnetic energy evolve over time. We start with a random initialization of all flow contributions. The equatorially anti-symmetric flow harmonics quickly decay and a (nearly) perfectly equatorially symmetric flow remains. Such a pure symmetry is a characteristic of solutions close to onset of convection.

At about $t = 0.04$ we initialize strong, randomly distributed magnetic modes which all first decay. However, after about $t = 0.1$ the equatorially symmetric flow starts to amplify the equatorially symmetric magnetic field. When this has reached sufficient strength at about $t = 0.4$, equatorially anti-symmetric flow and magnetic field modes also start to grow. The system finally settles into a solution where the r.m.s. equatorially anti-symmetric flow field is about a factor three weaker than its equatorially symmetric counterpart and the equatorially anti-symmetric magnetic field is about a factor 1.5 weaker than the symmetric one. Equatorially anti-symmetric flows are required to produce equatorially anti-symmetric magnetic field from equatorially symmetric field (and vice versa) (Gubbins & Zhang, 1993). This explains why the anti-symmetric field contributions grow together with the anti-symmetric flow (Figure 2a).

Figure 2b shows the evolution of the dominant Gauss coefficients defining the surface field. The axial dipole $g_1^0$ and axial quadrupole $g_2^0$ clearly dominate, and the equatorially anti-symmetric $g_1^0$ is about a factor three larger than the equatorially symmetric $g_2^0$. The ratio yields an offset $Z$ in the range of the one observed at Mercury. Axial dipole and axial quadrupole are also relatively stable in time while all other components keep on changing their sign. The time variability of the axial dipole is smaller than that of the axial quadrupole. These Gauss coefficients are given in non-dimensional form and their square would yield an Elsasser number. Since Mercury's Elsasser number is about $\Lambda \approx 10^{-5}$, the simulation seems to run at a realistic surface field strength. However, the internal field in the convective region reaches $\Lambda \approx 1$, as is expected for planetary dynamos. The Elsasser number at the CMB and Elsasser number in the convective region are listed in Table A1.

The snapshots in Figure 3 show the flow field at the equator and the radial magnetic field at the CMB and below the stratified layer. The left column illustrates the solutions during the kinematic growth phase of the magnetic field when Lorentz forces are still too small to modify the flow, while the right column illustrates the fully developed saturated solution. As seen in the equatorial slices in Figures 3a and 3b, the stable stratification effectively inhibits the flow from penetrating the outer region. The convective flow is very similar in both the growth and the saturated phase. The magnetic field, on the other hand, shows strong differences. During the kinematic growth phase, the CMB field is dominated by an axial quadrupole component (cf. Figure 3c), while in the saturated regime an axial dipole dominates with a sufficient axial quadrupole contribution to shift the magnetic equator significantly northward (cf. Figure 3d). Obviously, the initial convective flow preferably produces an equatorially symmetric magnetic field of the quadrupolar symmetry family. Once strong enough, the Lorentz force changes the flow in such a way that equatorially anti-symmetric flow and magnetic contributions
<figure><img src="image_6.png" /><figcaption>**Figure 2.** (a) Time evolution of the equatorial symmetric (es) and anti-symmetric (eas) contributions of the kinetic and magnetic energy at the beginning of the simulation DD2. (b) Time evolution of dipole offset $Z$ and the corresponding non-dimensional Gauss coefficients scaled by $\sqrt{\rho_0\mu_0\lambda^2\Omega}$ at the planetary surface up to the axial octupole $g_3^0$ for the DD2 model. For the computation of the coefficients at the planetary surface we assumed a ratio $r_o/R_M = 0.82$, which is a rough estimate from recent interior models (Genova et al., 2019; Knibbe et al., 2021; Steinbrügge et al., 2021). Non-axisymmetric Gauss coefficients ($m>0$) above degree $l=2$ are at least two orders of magnitude smaller and are not shown.</figcaption></figure>

can grow. A similar behavior has been described by Aubert and Wicht (2004), who studied axial and equatorial dipolar dynamos close to the onset of convection.

In the convective region, however, the magnetic fields are roughly similar and show a patchy structure somewhat concentrated at low latitudes (see Figures 3e and 3f). The axial dipole and quadrupole contributions only become dominant at the outer boundary of the saturated solution. The reason becomes apparent when analyzing the radial dependence of the different field contributions.

Figure 4 illustrates the variations of different kinetic and magnetic energy contributions with radius. Panels (a) and (c) show that the poloidal kinetic energy is damped very efficiently in the stable layer, while the toroidal kinetic energy remains relatively strong. The zonal flows, that is, the axisymmetric toroidal contributions (panel (a)), account for a significant part of the kinetic energy in the stable layer, but non-axisymmetric contributions remain strong. The equatorially anti-symmetric flows, which were allowed to grow by the presence of the magnetic field, remain subdominant at all radii (panel (c)). The exception are the toroidal flows in the very outer part of the stable layer.

Zonal flows produce an axisymmetric toroidal magnetic field by shearing the poloidal field by a so-called $\Omega$-effect (Parker, 1955). We attribute the significant rise in axisymmetric toroidal magnetic energy starting at $r = 0.4 r_o$ and its large amplitude in the stable layer to this effect (see Figure 4b). In the outer parts of the convective layer we thus expect a so called $\alpha\Omega$-dynamo at work. This expression refers to a mechanism where the poloidal field is produced by the small-scale convection while the toroidal field is predominantly produced by an $\Omega$-effect.
<figure><img src="image_7.png" /><figcaption>(b) Saturated Phase</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c) $B_r$ at $r/r_o = 1.00$ and $t_\lambda = 0.44$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) $B_r$ at $r/r_o = 1.00$ and $t_\lambda = 6.12$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(e) $B_r$ at $r/r_o = 0.30$ and $t_\lambda = 0.44$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(f) $B_r$ at $r/r_o = 0.30$ and $t_\lambda = 6.12$</figcaption></figure>

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

**Figure 3.** Equatorial slice of the radial velocity $u_r$ (a, b) and Hammer projection of the radial magnetic field $B_r$ at different radial levels (c–f) during the kinematic growth phase ($t_\lambda = 0.44$) and the saturated regime ($t_\lambda = 6.12$) of the magnetic field for the DD2 model. The dashed line in the equatorial slices indicates the bottom of the stably stratified layer $r_s$.

Panels (b) and (d) of Figure 4 show that the poloidal magnetic field is strongly attenuated over the stable layer where we expect negligible poloidal field production because of the weak radial flows. Two attenuation mechanisms are a magnetic skin effect due to the time dependence (Christensen, 2006; Christensen et al., 2018; Christensen & Wicht, 2008; Gastine & Wicht, 2021) and the action of zonal winds in the stable layer (Christensen & Wicht, 2008; Stanley et al., 2005; Stevenson, 1980; Wicht & Christensen, 2024). The former mainly affects non-axisymmetric field components that have high temporal variations. Skin effect and zonal wind action explain the steep decay in the lower part of the stable layer between $r = 0.5 r_o$ and $r = 0.65 r_o$.

The third reason for the attenuation is the geometric decay of a potential field like $(r_s/r)^{1+2}$. Over the stable layer from $r_s = 0.5 r_o$ to $r = r_o$ the dipole energy decays by a factor 0.0156 and the quadrupole energy by a factor 0.0039. The decrease in quadrupole energy is thus four times higher. Panel (d) of Figure 4 shows that in the convective region the equatorially symmetric magnetic field dominates its anti-symmetric counterpart. This changes in the stable layer since the symmetric contributions decay more severely likely due to the stronger decay of quadrupole than dipole contributions. Obviously, the combined attenuation effects strongly change the field produced in the convective region and lead to a surface field that is dominated by the equatorially anti-symmetric axisymmetric dipole $g_1^0$ and the equatorially symmetric axisymmetric quadrupole $g_2^0$.

Previous simulations have shown that the Ω-effect tends to promote more complex field geometries with stronger equatorially symmetric field contributions (Gastine et al., 2012; Schrinner et al., 2012). In these models, stress-free flow boundary condition are required to allow for fast zonal flows. While we use rigid boundary conditions,
<figure><img src="image_13.png" /><figcaption>**Figure 4.** Top: Time-averaged radial profiles of the (axisymmetric) poloidal and toroidal parts of the kinetic (a) and magnetic energy (b) for the model DD2. Bottom: Equatorially (anti-)symmetric poloidal and toroidal parts of the kinetic (c) and magnetic energy (d) are shown. The vertical dotted line in the plots indicates the bottom of the stably stratified layer.</figcaption></figure>

the stable layer largely decouples the convective region from the outer boundary and therefore allows faster zonal winds to develop. The strong Ω-effect is thus a possible reason for the initial production of an equatorially symmetric magnetic field in the growth phase.

To better understand the zonal winds responsible for the Ω-effect we analyze the zonal force balance, following Gastine and Wicht (2021), Wulff et al. (2022, 2024). The time- and azimuthally averaged φ-component of the Navier-Stokes Equation 4 in cylindrical coordinates (s, φ, z) is given by

$$0 = \overline{\langle F_{Ad} \rangle}_{\phi} + \overline{\langle F_{R} \rangle}_{\phi} + \overline{\langle F_{C} \rangle}_{\phi} + \overline{\langle F_{L} \rangle}_{\phi} + \overline{\langle F_{L} \rangle}_{\phi}; \qquad (20)$$

$$F_{AD} = \frac{1}{s} \nabla \cdot (s \langle \boldsymbol{u} \rangle_{\phi} \langle u_{\phi} \rangle_{\phi}), \qquad (21)$$

$$F_{R} = \frac{1}{s} \nabla \cdot (s \langle \boldsymbol{u}' \boldsymbol{u}_{\phi} \rangle_{\phi}), \qquad (22)$$

$$F_C = \frac{2 Pm}{E} \langle u_s \rangle_{\phi}, \qquad (23)$$
<figure><img src="image_14.png" /><figcaption>Figure 5. Meridional cuts of the time-averaged zonal flow $u_\phi$ and the single contributions of the zonal force balance in Equation 20 during (a) the growth phase of the magnetic field and (b) the saturated phase of model DD2.</figcaption></figure>

$$F_{\nu} = -\frac{Pm}{s} \nabla \cdot \left[ s^2 \nabla \left( \frac{\langle u_{\phi} \rangle_{\phi}}{s} \right) \right], \qquad (24)$$

$$F_L = -\frac{Pm}{E s} \nabla \cdot (s \langle BB_\phi \rangle_\phi). \qquad (25)$$

Here the velocity field is decomposed in its axisymmetric and non-axisymmetric part $u = \langle u \rangle_\phi + u'$, resulting in an advective force $F_{AD}$ and forces due to Reynolds stresses $F_R$. The other forces are the Coriolis force $F_C$, the viscous forces $F_\nu$ and the Lorentz force $F_L$, represented by the Maxwell stress tensor. In the following analysis, we neglect the advective force $F_{AD}$, since it is rather small in comparison to the other forces.
A comparison of the winds during the growth and saturated phase in Figure 5 (left panels) shows the symmetry-breaking by the equatorial anti-symmetric contribution in the saturated case most strongly in the outer parts of the convective region (panel (b)). In both cases, the zonal winds have deeply penetrated into the stable layer at all latitudes and are decisively non-geostrophic, that is, they have strong z-gradients.

The zonal forces are dominated by Coriolis and Reynolds stresses outside of the tangent cylinder. At and inside the tangent cylinder, viscous forces also contribute surprisingly strongly, likely due to the strong z-gradients. Lorentz forces only become important in the equatorial plane and at the transition to the stable layer in the northern hemisphere. This obvious north/south asymmetry relates to the asymmetry in magnetic field, zonal flows and the other zonal forces. However, the Lorentz force seems too weak and too localized to play any important role in zonal-wind driving or the equatorial symmetry breaking.

The variation of zonal winds in the direction of the rotation axis is described by the thermal and compositional zonal wind equation:

$$ \frac{\partial \langle u_{\phi} \rangle_{\phi}}{\partial z} = \frac{Ra_{\theta} E Pm}{2 Pr r_o} \frac{\partial \langle T \rangle_{\phi}}{\partial \theta} + \frac{Ra_{\xi} E Pm}{2 Sc r_o} \frac{\partial \langle \xi \rangle_{\phi}}{\partial \theta}. \quad (26) $$

Figure 6 compares left and right-hand sides of Equation 26. Inside the convective region, latitudinal gradients in the compositional field dominate. Here the light elements rise from the ICB in the form of confined plumes inside the tangent cylinder and around the equatorial plane. Latitudinal temperature gradients become important at mid to high latitudes at the transition to and in the stable layer and determine the respective wind structure. Low latitude winds in the stable layer are mostly driven by latitudinal compositional gradients. The convection seems vigorous enough to equilibrate the temperature deeper inside the convective region but not closer to the stable layer. The equatorial symmetry-breaking of the flow in the saturated phase is also apparent in the lower left panel of Figure 6.

We now investigate to which degree the results depend on the convective forcing by changing $Ra_{\theta}$ and $Ra_{\xi}$. We also vary $Pm$ to change the impact of the Lorentz force. The parameters for the different simulations along with key characteristics are listed in Table A1.

Figure 7 illustrates how the time-averaged quadrupole-to-dipole ratio $g_2^0/g_1^0$ and dipole tilt $\Delta_{dip}$ (color coded) depend on the Reynolds number $Re_{conv} = Rm_{conv}/Pm$ (left panel) and on the magnetic Reynolds number $Rm_{conv}$ (right panel) where the subscripts “conv” accounts for the fact that those diagnostics have been computed in the convection zone only. Time variations in $g_2^0/g_1^0$ and Reynolds number are indicated by error bars with the width of one standard variation. In panels (a) and (b) we keep thermal Rayleigh number $Ra_{\theta}$ fixed but vary $Ra_{\xi}$. In panel (c) $Ra_{\xi}$ is kept fixed while varying $Ra_{\theta}$.

Generally, the cases can be divided into two groups: Mercury-like solutions similar to DD2 where the mean $g_2^0/g_1^0$ and dipole tilt $\Delta_{dip}$ are relatively small, that is, $g_2^0/g_1^0 < 1.0$ and $\Delta_{dip} < 5^{\circ}$, and quadrupole-dominated solutions with $g_2^0/g_1^0 > 1.0$ and $\Delta_{dip} > 15^{\circ}$. The latter cases are characterized by a highly time variable reversing axial dipole while the axial quadrupole remains comparatively stable. This explains the strong time dependence of dipolar tilt, $g_2^0/g_1^0 > 1.0$, and Reynolds number. Table A1 lists mean and standard variation. The cases in the group of high $g_2^0/g_1^0$ and high $\Delta_{dip}$ show a better collapse when plotted against $Re_{conv}$ than against $Rm_{conv}$. This suggests a marginal influence of $Pm$ in this regime.

Increasing $Ra_{\xi}$ and thus the compositional convective driving causes a transition from the first group to the second, for example, for $Ra_{\theta} = 5 \times 10^6$ (panel (a)) when increasing $Ra_{\xi}$ from $2 \times 10^8$ (DD5) to $2.5 \times 10^8$ (DD6) at $Re_{conv} \approx 140$. Increasing $Ra_{\theta}$ to $1.0 \times 10^7$ (panel (b)) shifts this transition, which now happens at $Re_{conv} \approx 170$ for $Ra_{\xi}$ between $2.5 \times 10^8$ (DD15) and $Ra_{\xi} = 3 \times 10^8$ (DD16). However, panel (c) of Figure 7 shows that an increase of $Ra_{\theta}$ for fixed $Ra_{\xi}$ will also ultimately lead from Mercury-like to reversing cases. This somewhat complex behavior may have to do with the fact that the thermal Rayleigh number scales the stabilizing convective effects but also the thermal driving closer to the inner boundary, where the thermal buoyancy power highly increases with $Ra_{\theta}$ due to the small Prandtl number $Pr = 0.1$ (Simitev & Busse, 2005).
<figure><img src="image_15.png" /><figcaption>(b) Saturated Phase</figcaption></figure>

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

**Figure 6.** Thermo-compositional wind balance of the reference model DD2 during the kinematic growth phase (a) and the saturated phase (b). Displayed are the left and right hand side of the thermo-compositional wind Equation 26, respectively.

The number of reversals in each simulation, $N_{\text{rev}}$, and the mean time between reversals, $T_{\text{rev}}$, are listed in Table A1. The reversal frequency increases with either Rayleigh number. The higher time dependence yields a stronger skin effect, leading to a smaller relative amplitude of the $g_1^0$ coefficient and thus a larger mean ratio $g_2^0/g_1^0$.
<figure><img src="image_17.png" /><figcaption>**Figure 7.** The quadrupole-to-dipole ratio $|g_2^0/g_1^0|$ against the Reynolds number $Re_{conv}$ (left) and the magnetic Reynolds number $Rm_{conv}$ (right) in the convection zone for different values of constant $Ra_\theta$ (a, b) or $Ra_\xi$ (c). The color represents the dipole tilt $\Delta_{dip}$, whereas the error bars represent the standard deviation of $Re_{conv}$ and $|g_2^0/g_1^0|$. The hatched area indicates the range of possible values of $|g_2^0/g_1^0|$ at the present-day Mercury resulting from various internal magnetic field models (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019, 2021).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**Figure 8.** Snapshots of the radial magnetic field at the core-mantle boundary in the Hammer projection (upper left panel in each case), the radial magnetic field in the convection zone (lower left panel in each case) and a meridional cut of the azimuthally averaged zonal flow (right panel) for (a) model DD19 where thermal and compositional Rayleigh numbers were increased and (b) model DD3 where the magnetic Prandtl number for the reference model DD2 was increased to *Pm* = 5.</figcaption></figure>

Figure 8a illustrates DD19, a case with a strong time dependence and thus a large mean tilt and a large mean $g_2^0/g_1^0$. In the convective region, both magnetic symmetries are produced (see radial field in Figure 8a) with the equatorial symmetric contribution dominating, which is similar to the saturated phase of DD2. As already
discussed, from the simple geometric decay we would expect that $g_2^0$ decays four times more than $g_1^0$ over the stable layer. However, the skin effect yields a stronger decay of the axial dipole and results in a quadrupole-dominated field at the outer boundary.

Like in case DD2, the zonal winds, shown in Figure 8a, penetrate deeply into the stable layer but show a more complex and equatorially asymmetric pattern due to the larger Rayleigh number. Purely hydrodynamical simulations without a magnetic field have a similar structure but are more equatorially symmetric and geostrophic (not shown).

A special case is DD3, which has the same parameters as DD2 but $Pm = 5$ instead of $Pm = 3$. This amplifies the Lorentz force, enabling a classical dipole-dominated solution. Case DD3 appears as the left-most point in panels (a) and (c) of Figure 7 with practically zero $g_2^0/g_1^0$ value. Figure 8b illustrates the solution and shows that the equatorially symmetric field is already much smaller in the convective region than for DD2. Since the $g_2^0$ component also becomes highly time dependent, irregularly changing its sign, it is strongly attenuated by the skin effect and the surface field is by far dominated by the axial dipole. Additionally, the zonal winds are more altered by the stronger Lorentz force. While for DD2 prograde winds were present in mid-latitude at the bottom of the stable layer, they are now strongly attenuated in the strong dipolar case.

Generally, increasing $Pm$ yields stronger magnetic fields and thus likely stronger Lorentz forces, but this effect does not help to re-establish a Mercury-like solution at higher Rayleigh numbers.

Overall Mercury-like solutions for the double-diffusive model are rather limited in a narrow parameter regime close to the onset of convection. In the next section we show that a single-diffusive model approach can serve as a suitable alternative in order to achieve a persistent Mercury-like offset dipole solution.

## 3.2. The Single-Diffusive Model

Some authors argue that Mercury's dynamo is solely driven by thermal convection (Knibbe & van Westrenen, 2018). The different roles of thermal and compositional convection may therefore be only an unnecessary complication and a co-density approach (Braginsky & Roberts, 1995), or rather a purely thermal approach, may be more realistic. In this section we present cases that keep the inner core size and stable layer thickness from the double-diffusive model but are purely driven by thermal convection. All respective cases are listed in Table A2.

We start with discussing our reference single-diffusive model TH11 which successfully reproduces a Mercury-like field structure. The convection is driven by a heat flux released at the ICB, and the stably stratified layer is introduced by the background temperature gradient profile given by Equation 12.

After initiating the simulation with a random magnetic seed field, the dominant structure at the CMB during the kinematic growth phase is once more an axial quadrupole (not shown). This changes to a Mercury-like configuration in the saturated phase due to the impact of Lorentz forces. Figure 9 shows snapshots of the saturated magnetic field and azimuthal flow for the single-diffusive reference case TH11. Since the axial dipole switches polarity, we show snapshots for both dipole signs.

Because of the larger Rayleigh number, both equatorial symmetric and anti-symmetric flow modes are always present. Compared to case DD2, the convective flow and magnetic field inside the convective region are considerably more complex (see Figures 9c and 9d). Furthermore, the zonal flow is much more geostrophic though an equatorial anti-symmetric component is also clearly apparent. This component switches its sign with the axial dipole: Panels (e) and (f) demonstrate that the prograde flows at mid latitudes are stronger in the northern or southern hemisphere respectively.

Figure 10 shows the evolution of Gauss coefficients, of the dipole offset $Z$ and of the dipole tilt $\Delta_{dip}$. Like in the large Rayleigh number model DD19, the axial quadrupole always keeps its sign, which causes the magnetic equator to switch hemispheres. Other coefficients are rather small due to their decay over the stable layer. The offset matches the expected value for the present-day Mercury field quite well during epochs when the axial dipole remains roughly stable. The average dipole tilt is only 0.55°, and the tilt typically remains below 1° and thus within the range suggested for Mercury (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019, 2021). Since the reversals are rare and last much shorter than the inter-reversal periods, both the Mercury-like offset of the magnetic equator and small dipole tilts can be found during about 80% of the time. In the models
<figure><img src="image_19.png" /></figure>

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

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

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

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

<figure><img src="image_24.png" /><figcaption>**Figure 9.** Snapshots of the radial magnetic field $B_r$ in the Hammer projection at different radial levels and meridional cuts of the azimuthal flow component for the reference model TH11. (a, c) The radial magnetic field during the saturated phase of the field at the core-mantle boundary (CMB) and in the convective zone, respectively. (e) Meridional cut of the azimuthally averaged zonal flow component during the saturated phase. (b, d) Radial magnetic field during the saturated phase after a dipole reversal at the CMB and in the convective zone, respectively. (f) Meridional cut of the azimuthally averaged zonal flow component after a dipole reversal.</figcaption></figure>

presented by Christensen and Wicht (2008) all components constantly changed sign and the Mercury-like offset was only present for some instances (Wicht & Heyner, 2014). The square of the Gauss coefficients yields a realistic surface Elsasser number of about $\Lambda = 10^{-6}$, while the Elsasser number in the convective region is of order one. Table A2 lists both numbers for all our models.

The radial profiles of different magnetic and kinetic energy contributions are shown in Figure 11. They reveal a similar radial dependence as in the double-diffusive reference case DD2. The poloidal flow is damped efficiently by the stably stratified layer. While the non-axisymmetric toroidal flow also decreases significantly, the zonal flows seem largely unaffected and penetrate deep into the stable layer. Compared to the case DD2, the kinetic energy in the convection zone is almost one order of magnitude higher due to the larger Rayleigh number (larger super-criticality). In the stable layer, strong zonal winds are present and sustain an $\Omega$-effect.

Once more, the poloidal magnetic field is effectively attenuated by the magnetic skin effect and by the zonal wind action. These effects are more pronounced than in DD2 because the poloidal field is of smaller scale, more time-
<figure><img src="image_25.png" /><figcaption>**Figure 10.** (a) Time evolution of the non-dimensional Gauss coefficients scaled by $\sqrt{\rho_0\mu_0\lambda\Omega}$ up to the axial octupole $g_1^0$. For the computation of the coefficients at the planetary surface we assumed a ratio $r_o/R_M = 0.82$, which is a rough estimate from recent interior models (Genova et al., 2019; Knibbe et al., 2021; Steinbrügge et al., 2021). Non-axisymmetric Gauss coefficients ($m>0$) above degree $\ell = 2$ are at least two orders of magnitude smaller and are not shown. (b) Time evolution of the dipole offset $\mathcal{Z}$ and the dipole tilt $\Delta_{\text{dip}}$ for the purely thermal model TH11. The hatched areas indicate the range for the dipole offset from present-day Mercury magnetic field models (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019, 2021).</figcaption></figure>

dependent, and more non-axisymmetric. It nevertheless becomes nearly perfectly axisymmetric in the stable layer. At $r \approx 0.8$, the dipole component starts to dominate over the quadrupole component.

Figure 12 shows the axial quadrupole-to-dipole ratio $g_2^0/g_1^0$ against $Re_{\text{conv}}$ (panel (a)) and $Rm_{\text{conv}}$ (panel (b)) for all single-diffusive simulation runs. Like in Figure 7, the symbol color indicates the dipole tilt and the hatched area highlights the possible values for Mercury. Choosing $Rm_{\text{conv}}$ rather than $Re_{\text{conv}}$ yields a better collapse of the cases. This indicates a stronger impact of the magnetic Prandtl number and thus of magnetic effects. Nearly all the cases have reversing dipoles.

Focusing on $Pm = 3$ (circles), we observe that for $Rm_{\text{conv}} < 900$ the quadrupole-to-dipole ratio becomes too high. Close to the onset of dynamo action at $Rm_{\text{conv}} \approx 700$, the fundamental equatorially symmetric magnetic solution clearly dominates and the CMB field is close to a pure axial quadrupole. However, as the Rayleigh number and thus $Rm_{\text{conv}}$ increase, the equatorial anti-symmetric flow components (indicated by an increase of $E_{k,\text{eas}}/E_{k,\text{es}}$ in
<figure><img src="image_26.png" /><figcaption>**Figure 11.** Time-averaged radial profiles of the (axisymmetric) poloidal and toroidal parts of the kinetic (a) and magnetic energy (b) for the model TH11. Below the equatorially (anti-) symmetric poloidal and toroidal parts of the kinetic (c) and magnetic energy (d) are shown. The vertical dotted line in the plots indicates the bottom of the stably stratified layer.</figcaption></figure>

Table A2) and equatorially anti-symmetric magnetic field components grow. A present-day Mercury-like quadrupole-to-dipole ratio is reached in the range $880 \le Rm_{\text{conv}} \le 1050$.

When $Ra$ is increased even further, the quadrupole-to-dipole ratio becomes too small. The dependence of the ratio on $Ra$ is mostly explained by the fact that stronger convective driving promotes a larger axial dipole while the axial quadrupole increases more mildly. Variations in the skin effect only play a small role since the reversal rate remains comparable over the Rayleigh number range explored here while the axial quadrupole never reverses in any of the cases.

Figure 12 also shows an outlier where for $Pm = 5$ (squares) and $Rm_{\text{conv}} \approx 940$ the field is still strongly quadrupolar. When increasing the Rayleigh number and thus $Rm$ further the $Pm = 5$ cases roughly collapse with the $Pm = 2$ (triangles) and $Pm = 3$ (circles) solutions. Solutions become slightly more dipolar when increasing $Pm$ (cf. Table A2). This effect is less drastic than close to the onset of dynamo action in the double-diffusive case (cf. Figure 7). The $Pm$ dependence needs to be explored further. For the lower magnetic Prandtl numbers and hence weaker Lorentz forces the parameter space for a Mercury-like offset dipole is wider.

Overall the mean dipole tilt is rather small, especially for the cases where the quadrupole-to-dipole ratio is in a range expected for Mercury. Mean tilts smaller than 1°, consistent with Mercury field models, in combination
<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /><figcaption>**Figure 12.** Quadrupole-to-dipole ratio $|g_2^0/g_1^0|$ against (a) the Reynolds number $Re_{conv}$ in the convection zone and (b) the magnetic Reynolds number $Rm_{conv}$. The color bar shows the mean dipole tilt. The hatched area indicates the range of possible values of $|g_2^0/g_1^0|$ at the present-day Mercury resulting from various internal magnetic field models (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019, 2021).</figcaption></figure>

with a Mercury-like offset of the magnetic equator are found for several cases. The tilt standard deviation amounts to a few degree while the relative standard deviation in the offset is about 50% or smaller. Examples are TH9, TH11, or TH14. See Table A2.

All single-diffusive cases discussed so far are driven by a heat flux from the inner boundary while heat sources in the double-diffusive models are purely volumetric. We also performed a small set of single-diffusive simulations where convection is driven purely volumetrically (see cases marked with a V at the end of the ID in Table A2) with the thermal background state represented by Equation 9. Particularly large magnetic Reynolds numbers of $Rm_{conv} > 1000$ are required to initiate dynamo action in these cases. For all these simulations the resulting magnetic field at the CMB is dominated by an axial dipole with a quadrupole-to-dipole ratio of about $g_2^0/g_1^0 = 0.1$. The simple reason is that the quadrupole is not any more the fundamental dynamo mode for volumetrically driven single-diffusive cases.

### 3.3. Ekman Number Dependence

All double-diffusive and single-diffusive models discussed so far have been performed with different Ekman numbers, that is, $E = 2 \times 10^{-4}$ for the double-diffusive cases and $E = 10^{-4}$ for the single-diffusive cases. For a better comparison of the different setups, we performed a small set of double-diffusive simulations at $E = 10^{-4}$ and single-diffusive simulations at $E = 2 \times 10^{-4}$ and varied $Ra$ with a fixed $Pm = 3$. In the double-diffusive configuration, we fixed the thermal Rayleigh number at $Ra_g = 1.25 \times 10^7$ and changed $Ra_\xi$ to explore the effect of the convective forcing. These simulations are marked with a trailing “E” in their ID provided in Tables A1 and A2. Figure 13a shows the quadrupole-to-dipole ratio $|g_2^0/g_1^0|$ against the magnetic Reynolds number $Rm_{conv}$ in the convection zone, where the symbol color indicates the dipole tilt $\Delta_{dip}$. The results reveal that the decrease of the Ekman number by a factor of two has already a strong impact on the model solution. In the double-diffusive configuration, the dynamo is dipole-dominated close to the onset of convection but the quadrupole-to-dipole ratio is too large as compared to Mercury. As $Ra_\xi$ increases, the quadrupole-to-dipole ratio decreases rapidly to values around $|g_2^0/g_1^0| \le 0.2$, which is below the observed ratio for Mercury (cf. simulations DD2E-DD4E in Table A1). The dipole tilt is in these cases below $1^\circ$. As $Ra_\xi$ is increased from $4 \times 10^8$ (DD4E) to $5 \times 10^8$ (DD5E), we observe a strong increase of the dipole energy, resulting in an almost vanishing quadrupole-to-dipole ratio. In this case, the Lorentz force has a

stronger influence, which yields a nearly perfect dipole field at the surface with a lower $Rm_{conv}$ for case DD5 than for DD4. The strong increase in magnetic energy is indicated by the increase of the Elsasser number in the convection zone and the rise of the equatorial antisymmetric magnetic energy (cf. Table A1). None of these extra double-diffusive simulations show any sign changes of the axial dipole and quadrupole, explaining the small variance of $|g_2^0/g_1^0|$. Compared to the double-diffusive cases with $E = 2 \times 10^{-4}$, the influence of the equatorial symmetric mode becomes nearly negligible and cannot be controlled by the Rayleigh number anymore. On the contrary, the influence of the quadrupole becomes even smaller with increasing Rayleigh number, which narrows the parameter regime for Mercury-like magnetic fields even further. Figure 13b shows that for the single-diffusive configuration the increase of the Ekman number $E = 2 \times 10^{-4}$ has clearly the opposite effect. Close to the onset of dynamo action at $Rm_{conv} \approx 500$ the magnetic field at the surface is dominated by axial quadrupolar field, indicated by $|g_2^0/g_1^0| > 1$. Although the ratio decreases with increasing Rayleigh number and the magnetic field becomes dipole-dominated at the surface, it settles at $|g_2^0/g_1^0| \approx 0.6$. All single-diffusive simulations feature
<figure><img src="image_29.png" /><figcaption>**Figure 13.** Quadrupole-to-dipole ratio against the magnetic Reynolds number $Rm_{\text{conv}}$ in the convection zone for (a) the double-diffusive model with lower Ekman number $E = 10^{-4}$ for constant $Ra_{\theta} = 1.25 \times 10^7$ and (b) for the single-diffusive model with a higher Ekman number $E = 2 \times 10^{-4}$. The color indicates the dipole tilt $\Delta_{\text{dip}}$, whereas the error bars represent the standard deviation of $Rm_{\text{conv}}$ and $|g_2^0/g_1^0|$. The hatched area indicates the range of possible values of $|g_2^0/g_1^0|$ at the present-day Mercury resulting from various internal magnetic field models (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019, 2021).</figcaption></figure>

dipole tilt greater than 1° and a high variance in $g_2^0/g_1^0$ due to frequent reversals of the axial dipole. Therefore, these simulations are not suitable for explaining the observed surface field at Mercury.

## 4. Discussion

We have considered two different dynamo configurations that yield Mercury-like quadrupole-to-dipole ratios and small dipole tilts: a double-diffusive setup identical to the dynamo model suggested by Takahashi et al. (2019) and a single-diffusive (thermal/co-density) setup. The latter one is an attempt to capture the essence of Mercury's magnetic field in a simpler model. Both setups feature a thick stably stratified layer that occupies 50% of the fluid outer core volume. Consequently, the convective layer occupies the innermost 30% of the total core radius. The effective Ekman number defined on the thickness of the convective layer hence reads $E_{\text{eff}} = E/0.375^2$, rising the effective values by a factor 7.1 to $7.1 \times 10^{-4}$ and $1.4 \times 10^{-3}$. A small Rossby number, $Ro = U/\Omega d \ll 1$, indicates that the dynamics is rotationally constrained. The typical Rossby numbers in our simulations never exceed $Ro_{\text{eff}} = Ro/0.375 \approx 0.1$. Despite the large effective Ekman numbers, our simulations are hence still dominated by rotation, which is also confirmed by the columnar structure of the solution. Effective magnetic Reynolds numbers for the onset of dynamo action in our setups amount to $Rm_{\text{eff}} \approx 77$ for the double-diffusive reference case and $Rm_{\text{eff}} \approx 230$ for the single-diffusive case. This is more in line with previously reported values for simulations without stable layers (Petitdemange, 2018). Flow velocity estimates for Mercury based on scaling laws predict magnetic Reynolds number of a few hundreds in agreement with our effective values (Christensen & Wicht, 2008; Knibbe & van Westrenen, 2018).

An essential ingredient for the Mercury-like field characteristics is the fact that dynamo action sets in with the equatorially symmetric quadrupole family in both cases. However, the equatorially anti-symmetric dipole family sets in through different mechanisms. In the double-diffusive case, the Lorentz force associated to the quadrupole family magnetic fields triggers the onset of equatorially anti-symmetric flow fields. These in turn produce the dipole family magnetic fields. This is reminiscent of a scenario described by Aubert and Wicht (2004) using a single-diffusive setup without stable layer. In our single-diffusive setup, both symmetry families are present in flow and magnetic field right at the beginning of dynamo action because these dynamos require higher Rayleigh numbers. Since the equatorially anti-symmetric family grows with increasing Rayleigh number, we can adjust this parameter to reach the desired quadrupole-to-dipole ratio.

With this study, it becomes evident that the small inner core size plays a crucial role in featuring the quadrupole symmetry family so prominent. Dynamo simulations by Cao et al. (2014) and Tian et al. (2015) successfully generated Mercury-like dipole-to-quadrupole ratios with a small inner core ($a \le 0.2$), while models with larger inner cores ($a \ge 0.35$) could not produce a

persistent Mercury-like ratio (Christensen & Wicht, 2008; Wicht & Heyner, 2014). Christensen and Wicht (2008) indeed report that a smaller inner core ($a = 0.35$) seems to promote the quadrupole family. Strong zonal flows may also be important. Several authors report that the $\Omega$-effect, powered by zonal flows, tends to produce more complex magnetic fields with strong or dominant quadrupole-family contributions (Aubert, 2005; Gastine et al., 2012; Schrinner et al., 2012). Fast zonal flows typically develop in simulations with stress-free boundaries. While we used rigid boundaries here, the stable layer decouples convective region and outer boundary and allows for relatively prominent zonal flows. Consequently, the $\Omega$-effect is quite active.
A second essential ingredient for a Mercury-like magnetic field is the modification over the stable layer. Three mechanisms are at play: a magnetic skin effect damps faster varying smaller scale contributions; second, the zonal flows in the stable layer attenuate non-axisymmetric features; and third, high harmonic contributions decay faster with distance to their source region. The right combination results in a Mercury-like surface field that is dominated by the axial dipole and the axial quadrupole in the right proportion. We could at least adjust the time dependence and thus the first mechanism with the Rayleigh number(s). The thickness of the stable layer and the speed of the zonal winds are two other factors which control the field modification but have not been explored here. The octupole-to-dipole ratio in our Mercury-like models is by almost a factor 4 too low when compared to observations. Using a thinner stable layer would boost the octupole but would also likely require to adjust the Rayleigh number to keep a favorable quadrupole-to-dipole ratio.

Christensen and Wicht (2008) had to run their simulations at a relatively high Rayleigh number to promote a sufficiently high contribution of the equatorially symmetry family. Since this results in a generally strong time dependence, the magnetic field is Mercury-like only during a few instances in time (Wicht & Heyner, 2014). This is in difference to our model where we start with dominantly equatorially symmetric field. At a Rayleigh number sufficient to yield the required axial dipole contribution, the axial quadrupole remains relatively stable while the axial dipole undergoes only rare reversals.

While a realistic Ekman number for Mercury would be $E = \mathcal{O}(10^{-13})$ (Seuren et al., 2023) we use $E = \mathcal{O}(10^{-4})$. Numerical expenses significantly increase with lowering $E$. The moderate values considered here allowed exploring the other parameter dependencies. A small factor two change in the Ekman number however reveals that low Ekman numbers tend to increase the influence of the equatorial anti-symmetric modes of the dipole family, while high Ekman numbers favor the equatorial symmetric modes. The change in Ekman number resulted in a complete different Rayleigh number dependence for the double-diffusive and single-diffusive model, respectively. These results are in line with Christensen and Wicht (2008) who observed that a lower Ekman number of $E = 3 \times 10^{-5}$ tends to promote the axial dipole at the surface of the planet. Lower Ekman numbers also allow stronger zonal winds to develop, which will yield stronger $\Omega$-effect and possibly more attenuation of non-axisymmetric features in the stable layer. Lower Ekman number simulation are definitely required to explore the stability of the Mercury-like regimes reported here.

How realistic are our Mercury-like setups? The single-diffusive approach would be strictly realistic when only thermal effects drive convection. This could be the case for a sulfur-poor Fe-Si core discussed by Knibbe and van Westrenen (2018) and Davies et al. (2024) because Si fractionates to similar degrees into solid and liquid phase. The latent heat released at Si crystallization would then be a main driver of convection in this scenario. If there is an appreciable amount of sulfur (and carbon) in the core—at least 2 wt% in order to explain the small inner core (Knibbe et al., 2021)—a double-diffusive approach seems more realistic. Our double-diffusive simulations show complex solutions close to onset of convection with surprisingly non-geostrophic zonal flows and a high sensitivity to any parameter change. But, a zero heat flux from the inner core is assumed. Takahashi et al. (2019) report that a realistic non-zero heat flux leads to a solution that is not Mercury-like any more. The interplay of both drivers of convection in this setup is not really understood. Some of the complication may only exist close to onset and would vanish for higher Rayleigh numbers. However, we have shown that the magnetic fields lose their Mercury-like character for stronger convective driving.

The single-diffusive model allows for a heat flow through the inner boundary and seems more realistic in this respect. Also, flows way beyond the onset of convection are needed ($Re_{conv} \gg 1$) in order to obtain Mercury-like solutions. The small inner core of $a = 0.2$, required for a Mercury-like magnetic field in this setup, is still allowed by the newest interior models though larger cores are preferred (Genova et al., 2019; Knibbe et al., 2021; Steinbrügge et al., 2021).

Currently, the thickness of the stably stratified layer is not well constrained by thermal evolution models. The core cooling controlled by mantle dynamics and the inner core growth play a role (Davies et al., 2024; Knibbe & Van Hoolst, 2021; Knibbe & van Westrenen, 2018). These mechanisms depend on core and mantle compositions. Thermal evolution models by Knibbe and Van Hoolst (2021) for an Fe-Si-Ni core suggest a stable layer thickness of around $d_s = 0.25 r_o$, which is half of the thickness assumed in our models. Additionally, the inner core can grow up to $r_i = 0.6 r_o$, which is larger compared to Mercury-like dynamo models. Most recently, the models by Davies et al. (2024), where some of the models include a FeS-layer, result in a stable layer thickness of roughly
$d_s = 0.25 - 0.35 r_o$ and $r_i = 0.5 - 0.6 r_o$, similar to Knibbe et al. (2021). Clearly, future thermal evolution and dynamo models will have to gradually converge to obtain a consistent idea of Mercury's core structure.

**5. Conclusion**

Mercury's surface magnetic field is dominated by an axial dipole and a relatively strong axial quadrupole contribution. This combination causes the characteristic northern offset of the magnetic equator (20% of the planetary radius) and the very small dipole tilt (below one degree) (Anderson et al., 2012; Thébault et al., 2018; Wardinski et al., 2019, 2021). Moreover, the field is rather weak, corresponding to an Elsasser number of only $\Lambda = 10^{-5}$ at the CMB.

While the dynamo model of Christensen and Wicht (2008) only captures these characteristics for a short period of time (below 10% of the magnetic diffusion time, corresponding to roughly 3 kyr), those of Cao et al. (2014) and Tian et al. (2015) require an unrealistic heat-flux pattern through the CMB. The double-diffusive dynamo model by Takahashi et al. (2019) is the only one that features a permanently Mercury-like field without relying on such heat-flux pattern. All these models incorporate a thick stably stratified layer at the top of the mantle.

Here we explored the model by Takahashi et al. (2019) in more detail. We find it to be very sensitive to the parameter choices. It requires a delicate balance of thermal and compositional driving close to onset of convection to yield a Mercury-like field. Furthermore, no heat flux through the inner boundary is allowed, a rather low convective driving is required, and the zonal flow is strongly non-geostrophic. In conclusion, the model seems to be an unlikely candidate for Mercury.

Therefore we present an alternative single-diffusive model, which represents a purely thermal driving or a simplified combination of thermal and compositional effects. In this model's reference case, long periods of roughly 1.4 magnetic diffusion times (corresponding to 114 kyr when the smaller aspect ratio compared to Christensen and Wicht (2008) is considered) with Mercury-like magnetic field are interrupted by reversals of the axial dipole component. With each reversal, the offset of the magnetic equator switches hemispheres. Driving from the inner core now is a prerequisite for a Mercury-like field. A smaller inner core also seems to be advantageous.

Dynamo action sets in with an equatorially symmetric magnetic field of the quadrupole family. Since the dipole family contribution increases with Rayleigh number, we can use this parameter to reach the desired quadrupole-to-dipole ratio. The small mean dipole tilt below one degree and an Elsasser number between $10^{-5}$ and $10^{-4}$ prove how well these models can explain Mercury's dynamo. However, more studies are required to explore the parameters dependence, in particular at lower Ekman numbers.

The BepiColombo mission (Benkhoff et al., 2021), currently en-route to Mercury, will provide new insights onto the planet's interior and yield a more complete understanding of the magnetic environment. Geophysical and geodetic measurements will help to further constrain the interior (Genova et al., 2021). The magnetometer experiment on board the planetary orbiter (Glassmeier et al., 2010; Heyner et al., 2021) will provide an updated map of Mercury's internal field and improved constraints on its secular variation that will allow additional validations of dynamo models.
Appendix A: Numerical Simulations

<caption>Table A1</caption>

Summary of Model Parameters and Simulation Diagnostics for the Double-Diffusive Model

<table><thead><tr><th>ID</th><th>ℓ<sub>max</sub></th><th>N<sub>r</sub></th><th>Ra<sub>θ</sub> (×10<sup>7</sup>)</th><th>Ra<sub>ζ</sub> (×10<sup>8</sup>)</th><th>Pm</th><th>Ñ<sub>o</sub></th><th>Rm<sub>conv</sub></th><th>Λ<sub>conv</sub> (×10<sup>-4</sup>)</th><th>Λ<sub>CMB</sub> (×10<sup>-4</sup>)</th><th>Z</th><th>σ(Z)</th><th>Δ in [degree]</th><th>σ(Δ) in [degree]</th><th>N<sub>rev</sub></th><th>T<sub>rev</sub></th><th>(E<sub>k,conv</sub>/E<sub>L,es</sub>)<sub>conv</sub></th><th>(E<sub>m,conv</sub>/E<sub>L,es</sub>)<sub>conv</sub></th></tr></thead><tbody><tr><td>DD1</td><td>85</td><td>101</td><td>0.5</td><td>0.75</td><td>3</td><td>4.47</td><td>225</td><td>1.83</td><td>1.92</td><td>0.09</td><td>0.02</td><td>1.29</td><td>0.75</td><td>0</td><td>∞</td><td>0.08</td><td>0.81</td></tr><tr><td>DD2<sup>†</sup></td><td>128</td><td>129</td><td>0.5</td><td>1</td><td>3</td><td>4.47</td><td>276</td><td>1.53</td><td>1.26</td><td>0.14</td><td>0.02</td><td>1.09</td><td>0.62</td><td>0</td><td>∞</td><td>0.09</td><td>0.53</td></tr><tr><td>DD3</td><td>85</td><td>101</td><td>0.5</td><td>1</td><td>5</td><td>4.47</td><td>373</td><td>11.35</td><td>38.21</td><td>0</td><td>0</td><td>1.02</td><td>0.52</td><td>0</td><td>∞</td><td>0.18</td><td>2.22</td></tr><tr><td>DD4</td><td>85</td><td>101</td><td>0.5</td><td>1.5</td><td>3</td><td>4.47</td><td>343</td><td>1.59</td><td>2.27</td><td>0.12</td><td>0.02</td><td>1.04</td><td>0.6</td><td>0</td><td>∞</td><td>0.12</td><td>0.71</td></tr><tr><td>DD5</td><td>85</td><td>101</td><td>0.5</td><td>2</td><td>3</td><td>4.47</td><td>392</td><td>1.57</td><td>1.44</td><td>0.16</td><td>0.04</td><td>1.89</td><td>1.89</td><td>0</td><td>∞</td><td>0.16</td><td>0.59</td></tr><tr><td>DD6</td><td>85</td><td>101</td><td>0.5</td><td>2.5</td><td>3</td><td>4.47</td><td>431</td><td>1.63</td><td>0.34</td><td>0.88</td><td>0.69</td><td>20.99</td><td>19.69</td><td>12</td><td>0.66</td><td>0.19</td><td>0.49</td></tr><tr><td>DD7</td><td>85</td><td>101</td><td>0.5</td><td>3</td><td>3</td><td>4.47</td><td>469</td><td>1.79</td><td>0.32</td><td>1.24</td><td>0.9</td><td>31.15</td><td>23.69</td><td>24</td><td>0.38</td><td>0.23</td><td>0.54</td></tr><tr><td>DD8</td><td>85</td><td>101</td><td>0.5</td><td>3</td><td>5</td><td>4.47</td><td>738</td><td>5.92</td><td>0.47</td><td>0.97</td><td>0.7</td><td>25.47</td><td>22.05</td><td>18</td><td>0.24</td><td>0.28</td><td>0.71</td></tr><tr><td>DD9</td><td>85</td><td>101</td><td>0.5</td><td>3</td><td>2</td><td>4.47</td><td>322</td><td>0.39</td><td>0.09</td><td>2.6</td><td>1.89</td><td>52.4</td><td>22.1</td><td>33</td><td>0.33</td><td>0.2</td><td>0.43</td></tr><tr><td>DD10</td><td>128</td><td>129</td><td>0.5</td><td>6</td><td>3</td><td>4.47</td><td>611</td><td>4.5</td><td>0.76</td><td>2.41</td><td>1.74</td><td>35.4</td><td>23.03</td><td>20</td><td>0.22</td><td>0.37</td><td>0.73</td></tr><tr><td>DD11<sup>†</sup></td><td>128</td><td>129</td><td>1</td><td>1</td><td>3</td><td>6.32</td><td>307</td><td>3.16</td><td>2.59</td><td>0.16</td><td>0.02</td><td>0.97</td><td>0.47</td><td>0</td><td>∞</td><td>0.08</td><td>0.49</td></tr><tr><td>DD12</td><td>85</td><td>101</td><td>1</td><td>1.5</td><td>3</td><td>4.47</td><td>399</td><td>5.25</td><td>3.56</td><td>0.22</td><td>0.03</td><td>1.68</td><td>0.91</td><td>0</td><td>∞</td><td>0.12</td><td>0.48</td></tr><tr><td>DD13</td><td>128</td><td>129</td><td>1</td><td>2</td><td>3</td><td>4.47</td><td>454</td><td>5.15</td><td>3.12</td><td>0.25</td><td>0.04</td><td>2.11</td><td>1.22</td><td>0</td><td>∞</td><td>0.14</td><td>0.5</td></tr><tr><td>DD14</td><td>85</td><td>101</td><td>1</td><td>2.5</td><td>3</td><td>4.47</td><td>500</td><td>4.93</td><td>1.98</td><td>0.35</td><td>0.15</td><td>5.32</td><td>8.11</td><td>3</td><td>2.55</td><td>0.16</td><td>0.48</td></tr><tr><td>DD15</td><td>85</td><td>101</td><td>1</td><td>2.5</td><td>5</td><td>4.47</td><td>782</td><td>10.93</td><td>4.11</td><td>0.2</td><td>0.04</td><td>2.62</td><td>1.49</td><td>0</td><td>∞</td><td>0.2</td><td>0.74</td></tr><tr><td>DD16</td><td>85</td><td>101</td><td>1</td><td>3</td><td>3</td><td>4.47</td><td>536</td><td>5.35</td><td>1.33</td><td>0.58</td><td>0.36</td><td>15.25</td><td>18.08</td><td>14</td><td>0.64</td><td>0.19</td><td>0.51</td></tr><tr><td>DD17</td><td>85</td><td>101</td><td>1</td><td>3.5</td><td>3</td><td>4.47</td><td>560</td><td>5.28</td><td>0.93</td><td>1.16</td><td>0.93</td><td>27.4</td><td>22.97</td><td>18</td><td>0.4</td><td>0.22</td><td>0.54</td></tr><tr><td>DD18</td><td>85</td><td>101</td><td>1</td><td>3.5</td><td>2</td><td>4.47</td><td>405</td><td>2.21</td><td>0.35</td><td>2.29</td><td>1.88</td><td>47.71</td><td>23.33</td><td>18</td><td>0.31</td><td>0.17</td><td>0.41</td></tr><tr><td>DD19</td><td>85</td><td>101</td><td>1</td><td>4</td><td>3</td><td>4.47</td><td>591</td><td>5.26</td><td>0.68</td><td>1.74</td><td>1.42</td><td>38.52</td><td>22.83</td><td>20</td><td>0.32</td><td>0.25</td><td>0.59</td></tr><tr><td>DD20</td><td>85</td><td>101</td><td>1</td><td>4</td><td>2</td><td>4.47</td><td>421</td><td>2.16</td><td>0.35</td><td>2.04</td><td>1.58</td><td>48.88</td><td>24.07</td><td>24</td><td>0.23</td><td>0.21</td><td>0.47</td></tr><tr><td>DD21</td><td>128</td><td>129</td><td>1</td><td>6</td><td>3</td><td>4.47</td><td>681</td><td>6.5</td><td>0.73</td><td>2.93</td><td>2.2</td><td>50.62</td><td>23.42</td><td>22</td><td>0.16</td><td>0.32</td><td>0.69</td></tr><tr><td>DD22</td><td>85</td><td>101</td><td>1.5</td><td>1</td><td>3</td><td>4.47</td><td>502</td><td>6.72</td><td>1.38</td><td>1.17</td><td>1.3</td><td>26.96</td><td>21.82</td><td>22</td><td>0.36</td><td>0.03</td><td>0.2</td></tr><tr><td>DD23</td><td>128</td><td>129</td><td>2</td><td>1</td><td>3</td><td>4.47</td><td>709</td><td>8.82</td><td>1.02</td><td>2.11</td><td>1.66</td><td>38.87</td><td>22.28</td><td>20</td><td>0.26</td><td>0.03</td><td>0.36</td></tr><tr><td>DD24</td><td>128</td><td>129</td><td>2</td><td>2</td><td>3</td><td>4.9</td><td>695</td><td>10.85</td><td>2.13</td><td>1.77</td><td>1.49</td><td>25.64</td><td>21.77</td><td>22</td><td>0.31</td><td>0.09</td><td>0.43</td></tr><tr><td>DD25</td><td>128</td><td>129</td><td>2</td><td>4</td><td>3</td><td>4.9</td><td>786</td><td>12.63</td><td>1.95</td><td>1.68</td><td>1.14</td><td>27.92</td><td>22.23</td><td>15</td><td>0.33</td><td>0.19</td><td>0.59</td></tr><tr><td>DD26</td><td>128</td><td>129</td><td>2</td><td>4</td><td>5</td><td>4.9</td><td>1242</td><td>25.08</td><td>2.27</td><td>1.36</td><td>0.98</td><td>26.18</td><td>19.94</td><td>8</td><td>0.29</td><td>0.21</td><td>0.72</td></tr><tr><td>DD27</td><td>128</td><td>129</td><td>2.5</td><td>1</td><td>3</td><td>4.47</td><td>930</td><td>10.08</td><td>0.47</td><td>1.01</td><td>0.77</td><td>28.12</td><td>20.07</td><td>13</td><td>0.31</td><td>0.06</td><td>0.59</td></tr><tr><td>DD28</td><td>128</td><td>129</td><td>3</td><td>1</td><td>3</td><td>4.46</td><td>1078</td><td>12.49</td><td>0.28</td><td>0.95</td><td>0.76</td><td>30.29</td><td>23.35</td><td>18</td><td>0.21</td><td>0.1</td><td>0.79</td></tr><tr><td>DD29</td><td>128</td><td>129</td><td>3</td><td>6</td><td>3</td><td>4.58</td><td>1114</td><td>19.61</td><td>1.04</td><td>1.69</td><td>1.36</td><td>43.61</td><td>24</td><td>16</td><td>0.17</td><td>0.26</td><td>0.78</td></tr><tr><td>DD30</td><td>128</td><td>129</td><td>3</td><td>9</td><td>3</td><td>4.58</td><td>1221</td><td>17.88</td><td>0.7</td><td>1.95</td><td>1.64</td><td>58.17</td><td>21.46</td><td>21</td><td>0.13</td><td>0.27</td><td>0.85</td></tr><tr><td>DD1E</td><td>85</td><td>97</td><td>1.25</td><td>1.5</td><td>3</td><td>4.47</td><td>297</td><td>0.36</td><td>0.03</td><td>0.32</td><td>0.03</td><td>2.84</td><td>1.66</td><td>0</td><td>∞</td><td>0.01</td><td>0.06</td></tr><tr><td>DD2E</td><td>85</td><td>97</td><td>1.25</td><td>2.5</td><td>3</td><td>4.47</td><td>336</td><td>2.62</td><td>2.50</td><td>0.08</td><td>0.02</td><td>0.84</td><td>0.49</td><td>0</td><td>∞</td><td>0.13</td><td>0.91</td></tr><tr><td>DD3E</td><td>85</td><td>97</td><td>1.25</td><td>3</td><td>3</td><td>4.47</td><td>374</td><td>2.58</td><td>2.25</td><td>0.09</td><td>0.01</td><td>0.87</td><td>0.50</td><td>0</td><td>∞</td><td>0.14</td><td>0.88</td></tr><tr><td>DD4E</td><td>128</td><td>97</td><td>1.25</td><td>4</td><td>3</td><td>4.47</td><td>446</td><td>2.19</td><td>2.03</td><td>0.09</td><td>0.02</td><td>0.98</td><td>0.62</td><td>0</td><td>∞</td><td>0.15</td><td>0.98</td></tr><tr><td>DD5E</td><td>128</td><td>97</td><td>1.25</td><td>5</td><td>3</td><td>4.47</td><td>447</td><td>5.09</td><td>10.24</td><td>0.00</td><td>0.00</td><td>0.73</td><td>0.38</td><td>0</td><td>∞</td><td>0.20</td><td>1.96</td></tr></tbody></table>

Note. All simulations were performed with $E = 2 \times 10^{-4}$, $Pr = 0.1$, $Sc = 1.0$ and $a = 0.2$. Simulations with $E = 10^{-4}$ are marked with an “E” at the end of their ID. Diagnostics with subscript conv are evaluated in the convection zone and with subscript CMB at the core-mantle-boundary. $σ$ stands for the standard deviation of the variable, whereas $N_{\text{rev}}$ is the number of dipole reversals and $T_{\text{rev}}$ in magnetic diffusion times the average period between two reversals. Simulations already investigated by Takahashi et al. (2019) are marked with a dagger.
<table><caption>**Table A2**
Summary of Model Parameters and Simulation Diagnostics for the Purely Thermal Model (See Table A1 for Details)</caption><thead><tr><th>ID</th><th>ℓ<sub>max</sub></th><th>N<sub>r</sub></th><th>Ra<br>(×10<sup>8</sup>)</th><th>Pm</th><th>Ñ<sub>o</sub></th><th>Rm<sub>conv</sub></th><th>Λ<sub>conv</sub></th><th>Λ<sub>CMB</sub><br>(×10<sup>-5</sup>)</th><th>Z</th><th>σ(Z)</th><th>Δ<br>in [degree]</th><th>σ(Δ)<br>in [degree]</th><th>N<sub>rev</sub></th><th>T<sub>rev</sub></th><th>(E<sub>L,es</sub>/E<sub>L,es</sub>)<sub>conv</sub></th><th>(E<sub>m,es</sub>/E<sub>m,es</sub>)<sub>conv</sub></th></tr></thead><tbody><tr><td>TH1*</td><td>85</td><td>101</td><td>0.5</td><td>3</td><td>1.58</td><td>298</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.00</td><td>—</td></tr><tr><td>TH2*</td><td>85</td><td>101</td><td>1</td><td>3</td><td>1.73</td><td>453</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.00</td><td>—</td></tr><tr><td>TH3*</td><td>85</td><td>101</td><td>1.5</td><td>3</td><td>3.46</td><td>528</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.02</td><td>—</td></tr><tr><td>TH4*</td><td>85</td><td>101</td><td>2</td><td>3</td><td>3.16</td><td>658</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.07</td><td>—</td></tr><tr><td>TH5</td><td>85</td><td>101</td><td>2</td><td>5</td><td>3.74</td><td>972</td><td>1.40</td><td>0.33</td><td>0.59</td><td>0.46</td><td>3.64</td><td>9.24</td><td>19</td><td>0.36</td><td>0.07</td><td>0.82</td></tr><tr><td>TH6</td><td>85</td><td>101</td><td>2.5</td><td>3</td><td>3.54</td><td>714</td><td>0.50</td><td>0.44</td><td>0.55</td><td>0.37</td><td>2.29</td><td>6.01</td><td>14</td><td>0.75</td><td>0.10</td><td>0.68</td></tr><tr><td>TH7</td><td>85</td><td>101</td><td>3</td><td>3</td><td>3.87</td><td>759</td><td>1.07</td><td>0.80</td><td>0.49</td><td>0.36</td><td>2.84</td><td>7.52</td><td>13</td><td>0.65</td><td>0.12</td><td>0.74</td></tr><tr><td>TH8</td><td>85</td><td>101</td><td>3</td><td>5</td><td>3.87</td><td>1148</td><td>5.06</td><td>6.61</td><td>0.07</td><td>0.03</td><td>0.14</td><td>0.09</td><td>0</td><td>∞</td><td>0.19</td><td>0.94</td></tr><tr><td>TH9</td><td>85</td><td>101</td><td>3.5</td><td>3</td><td>3.74</td><td>820</td><td>1.60</td><td>1.57</td><td>0.26</td><td>0.14</td><td>1.04</td><td>4.17</td><td>9</td><td>0.73</td><td>0.16</td><td>0.83</td></tr><tr><td>TH10</td><td>85</td><td>101</td><td>3.5</td><td>5</td><td>4.39</td><td>1222</td><td>6.73</td><td>8.47</td><td>0.08</td><td>0.02</td><td>0.10</td><td>0.06</td><td>0</td><td>∞</td><td>0.24</td><td>0.93</td></tr><tr><td>TH11</td><td>85</td><td>101</td><td>4</td><td>3</td><td>3.46</td><td>888</td><td>2.49</td><td>2.85</td><td>0.19</td><td>0.07</td><td>0.68</td><td>3.04</td><td>7</td><td>1.44</td><td>0.17</td><td>0.86</td></tr><tr><td>TH12</td><td>85</td><td>101</td><td>4</td><td>5</td><td>3.46</td><td>1341</td><td>9.35</td><td>12.76</td><td>0.06</td><td>0.01</td><td>0.10</td><td>0.05</td><td>0</td><td>∞</td><td>0.25</td><td>0.93</td></tr><tr><td>TH13</td><td>85</td><td>101</td><td>4</td><td>2</td><td>3.46</td><td>625</td><td>0.72</td><td>1.00</td><td>0.36</td><td>0.19</td><td>1.71</td><td>5.18</td><td>11</td><td>0.82</td><td>0.15</td><td>0.67</td></tr><tr><td>TH14</td><td>85</td><td>101</td><td>4.5</td><td>3</td><td>3.67</td><td>942</td><td>2.89</td><td>2.85</td><td>0.20</td><td>0.10</td><td>0.90</td><td>4.22</td><td>14</td><td>0.74</td><td>0.19</td><td>0.87</td></tr><tr><td>TH15</td><td>85</td><td>101</td><td>5</td><td>3</td><td>3.87</td><td>957</td><td>3.45</td><td>3.90</td><td>0.16</td><td>0.06</td><td>0.39</td><td>1.67</td><td>5</td><td>1.96</td><td>0.21</td><td>0.87</td></tr><tr><td>TH16</td><td>85</td><td>101</td><td>5.5</td><td>3</td><td>3.71</td><td>1030</td><td>4.28</td><td>5.48</td><td>0.14</td><td>0.05</td><td>0.40</td><td>2.13</td><td>4</td><td>1.71</td><td>0.21</td><td>0.90</td></tr><tr><td>TH17</td><td>85</td><td>101</td><td>6</td><td>3</td><td>3.87</td><td>1065</td><td>4.53</td><td>5.56</td><td>0.14</td><td>0.06</td><td>0.42</td><td>2.38</td><td>7</td><td>1.07</td><td>0.23</td><td>0.90</td></tr><tr><td>TH18</td><td>128</td><td>129</td><td>6.5</td><td>3</td><td>3.61</td><td>1116</td><td>5.23</td><td>6.65</td><td>0.12</td><td>0.05</td><td>0.28</td><td>1.41</td><td>2</td><td>2.58</td><td>0.22</td><td>0.93</td></tr><tr><td>TH19</td><td>128</td><td>129</td><td>7</td><td>2</td><td>3.50</td><td>821</td><td>2.37</td><td>3.37</td><td>0.23</td><td>0.12</td><td>1.05</td><td>4.00</td><td>7</td><td>0.70</td><td>0.20</td><td>0.81</td></tr><tr><td>TH20</td><td>128</td><td>129</td><td>7</td><td>3</td><td>3.50</td><td>1149</td><td>5.95</td><td>8.37</td><td>0.10</td><td>0.04</td><td>0.28</td><td>1.66</td><td>3</td><td>1.54</td><td>0.24</td><td>0.93</td></tr><tr><td>TH21</td><td>128</td><td>129</td><td>7.5</td><td>3</td><td>3.62</td><td>1203</td><td>6.02</td><td>7.47</td><td>0.10</td><td>0.05</td><td>0.26</td><td>1.37</td><td>4</td><td>0.88</td><td>0.24</td><td>0.93</td></tr><tr><td>TH22</td><td>128</td><td>129</td><td>7.5</td><td>2</td><td>3.62</td><td>838</td><td>2.63</td><td>3.28</td><td>0.23</td><td>0.12</td><td>0.99</td><td>4.01</td><td>9</td><td>0.71</td><td>0.21</td><td>0.82</td></tr><tr><td>TH23</td><td>128</td><td>129</td><td>8</td><td>2</td><td>3.46</td><td>875</td><td>2.88</td><td>3.85</td><td>0.23</td><td>0.13</td><td>1.08</td><td>4.12</td><td>10</td><td>0.61</td><td>0.20</td><td>0.83</td></tr><tr><td>TH1V*</td><td>85</td><td>101</td><td>10</td><td>3</td><td>3.54</td><td>934</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.07</td><td>—</td></tr><tr><td>TH2V</td><td>85</td><td>101</td><td>15</td><td>3</td><td>4.33</td><td>1065</td><td>1.09</td><td>0.25</td><td>0.05</td><td>0.04</td><td>0.64</td><td>2.88</td><td>8</td><td>1.07</td><td>0.17</td><td>1.07</td></tr><tr><td>TH3V</td><td>85</td><td>101</td><td>20</td><td>3</td><td>4.47</td><td>1169</td><td>3.18</td><td>1.59</td><td>0.04</td><td>0.02</td><td>0.12</td><td>0.07</td><td>0</td><td>∞</td><td>0.25</td><td>1.07</td></tr><tr><td>TH4V*</td><td>85</td><td>101</td><td>20</td><td>2</td><td>4.47</td><td>893</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>0.17</td><td>—</td></tr><tr><td>TH5V</td><td>128</td><td>129</td><td>25</td><td>3</td><td>4.33</td><td>1290</td><td>4.66</td><td>4.42</td><td>0.04</td><td>0.02</td><td>0.07</td><td>0.04</td><td>0</td><td>∞</td><td>0.27</td><td>1.16</td></tr><tr><td>TH6V</td><td>128</td><td>129</td><td>30</td><td>3</td><td>4.58</td><td>1384</td><td>6.45</td><td>9.46</td><td>0.04</td><td>0.01</td><td>0.05</td><td>0.03</td><td>0</td><td>∞</td><td>0.28</td><td>1.20</td></tr><tr><td>TH7V</td><td>128</td><td>129</td><td>30</td><td>2</td><td>4.58</td><td>1104</td><td>0.13</td><td>0.02</td><td>0.04</td><td>0.05</td><td>1.51</td><td>6.29</td><td>15</td><td>0.22</td><td>0.19</td><td>1.09</td></tr><tr><td>TH8V</td><td>128</td><td>129</td><td>35</td><td>3</td><td>4.58</td><td>1524</td><td>6.99</td><td>0.82</td><td>0.04</td><td>0.05</td><td>0.92</td><td>3.89</td><td>8</td><td>0.56</td><td>0.27</td><td>1.04</td></tr><tr><td>TH9V</td><td>128</td><td>129</td><td>40</td><td>3</td><td>4.47</td><td>1587</td><td>9.44</td><td>7.22</td><td>0.03</td><td>0.01</td><td>0.06</td><td>0.03</td><td>0</td><td>∞</td><td>0.30</td><td>1.12</td></tr><tr><td>TH1E</td><td>85</td><td>101</td><td>1</td><td>3</td><td>3.50</td><td>522</td><td>1.71</td><td>1.85</td><td>0.77</td><td>0.47</td><td>22.71</td><td>21.17</td><td>12</td><td>1.42</td><td>0.14</td><td>0.26</td></tr><tr><td>TH2E</td><td>85</td><td>101</td><td>1.5</td><td>3</td><td>3.50</td><td>644</td><td>3.25</td><td>3.80</td><td>0.33</td><td>0.17</td><td>7.96</td><td>12.10</td><td>6</td><td>1.87</td><td>0.21</td><td>0.37</td></tr><tr><td>TH3E</td><td>85</td><td>101</td><td>2.05</td><td>3</td><td>3.51</td><td>763</td><td>4.19</td><td>4.49</td><td>0.26</td><td>0.13</td><td>6.16</td><td>9.93</td><td>5</td><td>1.86</td><td>0.23</td><td>0.44</td></tr><tr><td>TH4E</td><td>85</td><td>101</td><td>2.5</td><td>3</td><td>3.50</td><td>852</td><td>5.09</td><td>3.20</td><td>0.31</td><td>0.19</td><td>8.63</td><td>14.21</td><td>4</td><td>2.59</td><td>0.22</td><td>0.47</td></tr><tr><td>TH5E</td><td>128</td><td>129</td><td>3</td><td>3</td><td>3.50</td><td>931</td><td>6.01</td><td>2.70</td><td>0.32</td><td>0.23</td><td>7.67</td><td>12.71</td><td>15</td><td>0.88</td><td>0.24</td><td>0.53</td></tr><tr><td>TH6E</td><td>128</td><td>129</td><td>4</td><td>3</td><td>3.50</td><td>1105</td><td>6.53</td><td>1.16</td><td>0.32</td><td>0.25</td><td>7.55</td><td>12.63</td><td>17</td><td>0.64</td><td>0.22</td><td>0.70</td></tr></tbody></table>

Note. All simulations were performed with $E = 10^{-4}$, $Pr = 1.0$ and $a = 0.2$. Simulations which are purely volumetrically driven are marked with a “V” at the end of the ID, whereas simulations with $E = 2 \times 10^{-4}$ are marked with an “E” at the end of the ID. Simulations were the dynamo fails to sustain a magnetic field are denoted by an asterisk.
## Appendix B: Determination of the Gauss Coefficients

We obtain an expression for the Gauss coefficients by assuming that outside the fluid shell the magnetic field is described by an electrically insulating environment, such that $\nabla \times \mathbf{B} = 0$, and therefore we find a magnetic scalar potential with $\mathbf{B} = -\nabla V$. Typically the magnetic scalar potential is expanded in Schmidt semi-normalized spherical harmonics, where

$$S_{\ell}^{m}(\cos \theta) = \left[ (2 - \delta_{0}^{m}) \frac{(\ell - m)!}{(\ell + m)!} \right]^{1/2} P_{\ell}^{m}(\cos \theta),$$

is the normalization formula of the associate Legendre polynomials $P_{\ell}^{m}$, with $\delta_{0}^{m}$ as Kronecker delta. Then, the scalar magnetic potential with reference to the planetary surface is

$$V = R_M \sum_{\ell=0}^{L_{\max}} \sum_{m=0}^{\ell} \left(\frac{R_M}{r}\right)^{\ell+1} (g_{\ell}^{m} \cos m\phi + h_{\ell}^{m} \sin m\phi) S_{\ell}^{m}(\cos \theta), \quad (B1)$$

where $g_{\ell}^{m}$ and $h_{\ell}^{m}$ are the real Gauss coefficients resulting from the expansion. The numerical code used in this study, decomposes the magnetic field into poloidal and toroidal scalars. In this case, the scalars are expanded in orthonormalized spherical harmonics, where

$$\mathcal{N}_{\ell}^{m}(\cos \theta) = \left[ \frac{2\ell + 1}{4\pi} \frac{(\ell - m)!}{(\ell + m)!} \right]^{1/2} P_{\ell}^{m}(\cos \theta) \quad (B2)$$

is the normalization formula for the associate Legendre polynomials. The resulting expansion coefficients, however, are complex. Therefore, we need to find a conversion formula to determine the Gauss coefficients from the given poloidal and toroidal expansion coefficients. Since we can only compare both potential descriptions of the magnetic field at the CMB, the insulating boundary conditions forces the toroidal expansion coefficients to be zero. Comparing for example, the radial magnetic field at CMB obtained from the magnetic potential description (B1) with $B_{r}$ obtained from the poloidal potential and expanding it in spherical harmonics, the Gauss coefficients can be calculated by

$$g_{\ell}^{0} = \ell \left( \frac{2\ell + 1}{4\pi} \right)^{1/2} \left( \frac{r_{o}^{\ell}}{R_{M}^{\ell+2}} \right) a_{\ell}^{0}, \quad (B3)$$

$$g_{\ell}^{m>0} = \ell \left( \frac{2(2\ell + 1)}{4\pi} \right)^{1/2} \left( \frac{r_{o}^{\ell}}{R_{M}^{\ell+2}} \right) a_{\ell}^{m}, \quad (B4)$$

$$h_{\ell}^{0} = 0, \quad (B5)$$

$$h_{\ell}^{m>0} = -\ell \left( \frac{2(2\ell + 1)}{4\pi} \right)^{1/2} \left( \frac{r_{o}^{\ell}}{R_{M}^{\ell+2}} \right) b_{\ell}^{m}.$$

Here $a_{\ell}^{m}$ and $b_{\ell}^{m}$ are the real and imaginary part of the poloidal expansion coefficients, respectively. A similar derivation can be found in Glatzmaier (2013), but with the usage of a Condon-Shortley phase.

## Data Availability Statement

All simulations were performed with the open-source code MagIC freely available under https://github.com/magic-sph/magic. Figures were made with Matplotlib version 3.8.2 (Caswell et al., 2023; Hunter, 2007). Data for reproducing the figures is available at Kolhey (2024).
Acknowledgments

This work used the supercomputer Phoenix and was supported by the Gauß-IT-Zentrum of the University of Braunschweig (GITZ). We are grateful to the GITZ supercomputer staff, especially Rashid Kratou and Fritz Lewall. The authors gratefully acknowledge the computing time granted by the Resource Allocation Board and provided on the supercomputer Lise and Emmy at NHR@ZIB and NHR@Göttingen as part of the NHR infrastructure. The calculations for this research were conducted with computing resources under the project nip00046. DH was supported by the German Ministerium für Wirtschaft und Klimaschutz and the German Zentrum für Luft-und Raumfahrt under contract 50QW2202. Open Access funding enabled and organized by Projekt DEAL.

References

Alexeev, I. I., Belenkaya, E. S., Slavin, J. A., Korth, H., Anderson, B. J., Baker, D. N., et al. (2010). Mercury's magnetospheric magnetic field after the first two MESSENGER flybys. *Icarus*, **209**(1), 23–39. https://doi.org/10.1016/j.icarus.2010.01.024

Alken, P., Thébault, E., Beggan, C. D., Amit, H., Aubert, J., Baerenzung, J., et al. (2021). International geomagnetic reference field: The thirteenth generation. *Earth Planets and Space*, **73**(1), 49. https://doi.org/10.1186/s40623-020-01288-x

Anderson, B. J., Johnson, C. L., Korth, H., Winslow, R. M., Borovsky, J. E., Purucker, M. E., et al. (2012). Low-degree structure in Mercury's planetary magnetic field. *Journal of Geophysical Research*, **117**(E12). https://doi.org/10.1029/2012JE004159

Aubert, J. (2005). Steady zonal flows in spherical shell dynamos. *Journal of Fluid Mechanics*, **542**(-1), 53–67. https://doi.org/10.1017/S0022112005006129

Aubert, J., & Wicht, J. (2004). Axial vs. equatorial dipolar dynamo models with implications for planetary magnetic fields. *Earth and Planetary Science Letters*, **221**(1), 409–419. https://doi.org/10.1016/S0012-821X(04)00102-5

Benkhoff, J., Murakami, G., Baumjohann, W., Besse, S., Bunce, E., Casale, M., et al. (2021). BepiColombo - Mission overview and science goals. *Space Science Reviews*, **217**(8), 90. https://doi.org/10.1007/s11214-021-00861-4

Boscarino, S., Pareschi, L., & Russo, G. (2013). Implicit-explicit Runge-Kutta schemes for hyperbolic systems and kinetic equations in the diffusion limit. *SIAM Journal on Scientific Computing*, **35**(1), A22–A51. https://doi.org/10.1137/110842855

Braginsky, S. I., & Roberts, P. H. (1995). Equations governing convection in Earth's core and the geodynamo. *Geophysical & Astrophysical Fluid Dynamics*, **79**(1–4), 1–97. https://doi.org/10.1080/03091929508228992

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. *Geophysical Research Letters*, **41**(12), 4127–4134. https://doi.org/10.1002/2014GL060196

Caswell, T. A., Lee, A., de Andrade, E. S., Droettboom, M., Hoffmann, T., Klymak, J., et al. (2023). Matplotlib: Visualization with Python. Zenodo. https://doi.org/10.5281/zenodo.10150955

Christensen, U. R. (2006). A deep dynamo generating Mercury's magnetic field. *Nature*, **444**(7122), 1056–1058. https://doi.org/10.1038/nature05342

Christensen, U. R., Cao, H., Dougherty, M. K., & Khurana, K. (2018). Saturn's magnetic field and dynamo. In K. H. Baines, F. M. Flasar, N. Krupp, & T. Stallard (Eds.), *Saturn in the 21st century* (pp. 69–96). Cambridge University Press. https://doi.org/10.1017/9781316227220.004

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

Christensen, U. R., & Wicht, J. (2008). Models of magnetic field generation in partly stable planetary cores: Applications to Mercury and Saturn. *Icarus*, **196**(1), 16–34. https://doi.org/10.1016/j.icarus.2008.02.013

Davies, C. J., Pommier, A., Greenwood, S., & Wilson, A. (2024). Thermal and magnetic evolution of Mercury with a layered Fe-Si(-S) core. *Earth and Planetary Science Letters*, **641**, 118812. https://doi.org/10.1016/j.epsl.2024.118812

Edmund, E., Morard, G., Baron, M. A., Rivoldini, A., Yokoo, S., Boccato, S., et al. (2022). The Fe-FeSi phase diagram at Mercury's core conditions. *Nature Communications*, **13**(1), 387. https://doi.org/10.1038/s41467-022-27991-9

Gastine, T., Aubert, J., & Fournier, A. (2020). Dynamo-based limit to the extent of a stable layer atop Earth's core. *Geophysical Journal International*, **222**(2), 1433–1448. https://doi.org/10.1093/gji/gga0250

Gastine, T., Duarte, L., & Wicht, J. (2012). Dipolar versus multipolar dynamos: The influence of the background density stratification. A&, **546**, A19. https://doi.org/10.1051/0004-6361/201219799

Gastine, T., & Wicht, J. (2021). Stable stratification promotes multiple zonal jets in a turbulent Jovian dynamo model. *Icarus*, **368**, 114514. https://doi.org/10.1016/j.icarus.2021.114514

Genova, A., Goossens, S., Mazarico, E., Lemoine, F. G., Neumann, G. A., Kuang, W., et al. (2019). Geodetic evidence that Mercury has A solid inner core. *Geophysical Research Letters*, **46**(7), 3625–3633. https://doi.org/10.1029/2018GL081135

Genova, A., Hussmann, H., Van Hoolst, T., Heyner, D., Iess, L., Santoli, F., et al. (2021). Geodesy, geophysics and fundamental physics investigations of the BepiColombo mission. *Space Science Reviews*, **217**(2), 1–62. https://doi.org/10.1007/s11214-021-00808-9

Glassmeier, K.-H., Auster, H.-U., Heyner, D., Okafka, K., Carr, C., Berghofer, G., et al. (2010). The fluxgate magnetometer of the BepiColombo Mercury planetary orbiter. *Planetary and Space Science*, **58**(1), 287–299. https://doi.org/10.1016/j.pss.2008.06.018

Glassmeier, K.-H., Auster, H.-U., & Motschmann, U. (2007). A feedback dynamo generating Mercury's magnetic field. *Geophysical Research Letters*, **34**(22). https://doi.org/10.1029/2007GL031662

Glatzmaier, G. (2013). Boundaries and geometries. In *Introduction to modeling convection in planets and stars: Magnetic field, density stratification, rotation (chap. 10)*. Princeton University Press.

Gómez-Pérez, N., & Solomon, S. C. (2010). Mercury's weak magnetic field: A result of magnetospheric feedback? *Geophysical Research Letters*, **37**(20). https://doi.org/10.1029/2010GL044533

Gómez-Pérez, N., & Wicht, J. (2010). Behavior of planetary dynamos under the influence of external magnetic fields: Application to Mercury and Ganymede. *Icarus*, **209**(1), 53–62. https://doi.org/10.1016/j.icarus.2010.04.006

Gopinath, V., Fournier, A., & Gastine, T. (2022). An assessment of implicit-explicit time integrators for the pseudo-spectral approximation of Boussinesq thermal convection in an annulus. *Journal of Computational Physics*, **460**, 110965. https://doi.org/10.1016/j.jcp.2022.110965

Gubbins, D., & Zhang, K. (1993). Symmetry properties of the dynamo equations for palaeomagnetism and geomagnetism. *Physics of the Earth and Planetary Interiors*, **75**(4), 225–241. https://doi.org/10.1016/0031-9201(93)90003-R

Hauck, S. A., Margot, J.-L., Solomon, S. C., Phillips, R. J., Johnson, C. L., Lemoine, F. G., et al. (2013). The curious case of Mercury's internal structure. *Journal of Geophysical Research: Planets*, **118**(6), 1204–1220. https://doi.org/10.1002/jgre.20091

Heimpel, M., Aurnou, J., Al-Shamali, F., & Gomez Perez, N. (2005). A numerical study of dynamo action as a function of spherical shell geometry. *Earth and Planetary Science Letters*, **236**(1), 542–557. https://doi.org/10.1016/j.epsl.2005.04.032

Heyner, D., Auster, H.-U., Fornaçon, K.-H., Carr, C., Richter, I., Mieth, J. Z. D., et al. (2021). The BepiColombo planetary magnetometer MPO-MAG: What can we learn from the Hermean magnetic field? *Space Science Reviews*, **217**(4), 52. https://doi.org/10.1007/s11214-021-00822-x

Heyner, D., Schmitt, D., Glassmeier, K.-H., & Wicht, J. (2011). Dynamo action in an ambient field. *Astronomische Nachrichten*, **332**(1), 36–42. https://doi.org/10.1002/asna.201011466

Heyner, D., Schmitt, D., Wicht, J., Glassmeier, K.-H., Korth, H., & Motschmann, U. (2010). The initial temporal evolution of a feedback dynamo for Mercury. *Geophysical & Astrophysical Fluid Dynamics*, **104**(4), 419–429. https://doi.org/10.1080/03091921003776839

Heyner, D., Wicht, J., Gómez-Pérez, N., Schmitt, D., Auster, H.-U., & Glassmeier, K.-H. (2011). Evidence from numerical experiments for a feedback dynamo generating Mercury's magnetic field. *Science*, **334**(6063), 1690–1693. https://doi.org/10.1126/science.1207290

Hunter, J. D. (2007). Matplotlib: A 2D graphics environment. *Computing in Science & Engineering*, **9**(03), 90–95. https://doi.org/10.1109/MCSE.2007.55
Knibbe, J. S., Rivoldini, A., Luginbuhl, S. M., Namur, O., Charlier, B., Mezouar, M., et al. (2021). Mercury's interior structure constrained by density and P-wave velocity measurements of liquid Fe-Si-C alloys. *Journal of Geophysical Research: Planets*, 126(1), e2020JE006651.
https://doi.org/10.1029/2020JE006651

Knibbe, J. S., & Van Hoolst, T. (2021). Modelling of thermal stratification at the top of a planetary core: Application to the cores of Earth and Mercury and the thermal coupling with their mantles. *Physics of the Earth and Planetary Interiors*, 321, 106804. https://doi.org/10.1016/j.pepi.2021.106804

Knibbe, J. S., & van Westrenen, W. (2018). The thermal evolution of Mercury's Fe-Si core. *Earth and Planetary Science Letters*, 482, 147-159.
https://doi.org/10.1016/j.epsl.2017.11.006

Kolhey, P. (2024). Dataset: Dynamo models with a Mercury-like magnetic offset dipole [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.14536330

Kosloff, D., & Tal-Ezer, H. (1993). A modified Chebyshev pseudospectral method with an O(N-1) time step restriction. *Journal of Computational Physics*, 104(2), 457-469. https://doi.org/10.1006/jcph.1993.1044

Manglik, A., Wicht, J., & Christensen, U. R. (2010). A dynamo model with double diffusive convection for Mercury's core. *Earth and Planetary Science Letters*, 289(3), 619-628. https://doi.org/10.1016/j.epsl.2009.12.007

Monville, R., Vidal, J., Cébron, D., & Schaeffer, N. (2019). Rotating double-diffusive convection in stably stratified planetary cores. *Geophysical Journal International*, 219(Supplement_1), S195-S218. https://doi.org/10.1093/gji/ggz347

Nittler, L. R., Chabot, N. L., Grove, T. L., & Peplowski, P. N. (2018). The chemical composition of Mercury. In S. C. Solomon, L. R. Nittler, & B. J. Anderson (Eds.), *Mercury: The view after MESSENGER* (pp. 30-51). Cambridge University Press. https://doi.org/10.1017/9781316650684.003

Parker, E. N. (1955). Hydromagnetic dynamo models. *The Astrophysical Journal*, 122, 293. https://doi.org/10.1086/146087

Peplowski, P. N., Evans, L. G., Hauck, S. A., McCoy, T. J., Boynton, W. V., Gillis-Davis, J. J., et al. (2011). Radioactive elements on Mercury's surface from MESSENGER: Implications for the planet's formation and evolution. *Science*, 333(6051), 1850-1852. https://doi.org/10.1126/science.1211576

Petitdemange, L. (2018). Systematic parameter study of dynamo bifurcations in geodynamo simulations. *Physics of the Earth and Planetary Interiors*, 277, 113-132. https://doi.org/10.1016/j.pepi.2018.02.001

Pump, K., Heyner, D., Schmid, D., Exner, W., & Plaschke, F. (2024). Revised magnetospheric model reveals signatures of field-aligned current systems at Mercury. *Journal of Geophysical Research: Space Physics*, 129(3), e2023JA031529. https://doi.org/10.1029/2023JA031529

Renggli, C., Klemme, S., Morlok, A., Berndt, J., Weber, I., Hiesinger, H., & King, P. (2022). Sulfides and hollows formed on Mercury's surface by reactions with reducing S-rich gases. *Earth and Planetary Science Letters*, 593, 117647. https://doi.org/10.1016/j.epsl.2022.117647

Schaeffer, N. (2013). Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations. *Geochemistry, Geophysics, Geosystems*, 14(3), 751-758. https://doi.org/10.1002/ggge.20071

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

Seuren, F., Triana, S. A., Rekier, J., Barik, A., & Van Hoolst, T. (2023). Effects of the librationally induced flow in Mercury's fluid core with an outer stably stratified layer. *The Planetary Science Journal*, 4(9), 161. https://doi.org/10.3847/PSJ/acee77

Simitov, R., & Busse, F. H. (2005). Prandtl-number dependence of convection-driven dynamos in rotating spherical fluid shells. *Journal of Fluid Mechanics*, 532, 365-388. https://doi.org/10.1017/S0022112005004398

Smith, D. E., Zuber, M. T., Phillips, R. J., Solomon, S. C., Hauck, S. A., Lemoine, F. G., et al. (2012). Gravity field and internal structure of Mercury from MESSENGER. *Science*, 336(6078), 214-217. https://doi.org/10.1126/science.1218809

Solomon, S. C., Nittler, L. R., & Anderson, B. J. (2018). *Mercury: The view after MESSENGER* (Vol. 21). Cambridge University Press.

Stanley, S., Bloxham, J., Hutchison, W. E., & Zuber, M. T. (2005). Thin shell dynamo models consistent with Mercury's weak observed magnetic field. *Earth and Planetary Science Letters*, 234(1), 27-38. https://doi.org/10.1016/j.epsl.2005.02.040

Steinbrügge, G., Dumberry, M., Rivoldini, A., Schubert, G., Cao, H., Schroeder, D. M., & Soderlund, K. M. (2021). Challenges on Mercury's interior structure posed by the new measurements of its obliquity and tides. *Geophysical Research Letters*, 48(3), e2020GL089895. https://doi.org/10.1029/2020GL089895

Stevenson, D. J. (1980). Saturn's luminosity and magnetism. *Science*, 208(4445), 746-748. https://doi.org/10.1126/science.208.4445.746

Takahashi, F., & Matsushima, M. (2006). Dipolar and non-dipolar dynamos in a thin shell geometry with implications for the magnetic field of Mercury. *Geophysical Research Letters*, 33(10). https://doi.org/10.1029/2006GL025792

Takahashi, F., Shimizu, H., & Tsunakawa, H. (2019). Mercury's anomalous magnetic field caused by a symmetry-breaking self-regulating dynamo. *Nature Communications*, 10(1), 208. https://doi.org/10.1038/s41467-018-08213-7

Thébault, E., Langlais, B., Oliveira, J., Amit, H., & Leclercq, L. (2018). A time-averaged regional model of the Hermean magnetic field. *Physics of the Earth and Planetary Interiors*, 276, 93-105. https://doi.org/10.1016/j.pepi.2017.07.001

Tian, Z., Zuber, M. T., & Stanley, S. (2015). Magnetic field modeling for Mercury using dynamo models with a stable layer and laterally variable heat flux. *Icarus*, 260, 263-268. https://doi.org/10.1016/j.icarus.2015.07.019

Toepfer, S., Narita, Y., Glassmeier, K.-H., Heyner, D., Kolhey, P., Motschmann, U., & Langlais, B. (2021). The Mie representation for Mercury's magnetic field. *Earth Planets and Space*, 73(1), 65. https://doi.org/10.1186/s40623-021-01386-4

Tosi, N., Čadek, O., Běhounková, M., Káňová, M., Plesa, A.-C., Grott, M., et al. (2015). Mercury's low-degree geoid and topography controlled by insolation-driven elastic deformation. *Geophysical Research Letters*, 42(18), 7327-7335. https://doi.org/10.1002/2015GL065314

Tosi, N., Grott, M., Plesa, A.-C., & Breuer, D. (2013). Thermochemical evolution of Mercury's interior. *Journal of Geophysical Research: Planets*, 118(12), 2474-2487. https://doi.org/10.1002/jgre.20168

Vilim, R., Stanley, S., & Hauck, S. A. (2010). Iron snow zones as a mechanism for generating Mercury's weak observed magnetic field. *Journal of Geophysical Research*, 115(E11). https://doi.org/10.1029/2009JE003528

Wardinski, I., Amit, H., Langlais, B., & Thébault, E. (2021). The internal structure of Mercury's core inferred from magnetic observations. *Journal of Geophysical Research: Planets*, 126(12), e2020JE006792. https://doi.org/10.1029/2020JE006792

Wardinski, I., Langlais, B., & Thébault, E. (2019). Correlated time-varying magnetic fields and the core size of Mercury. *Journal of Geophysical Research: Planets*, 124(8), 2178-2197. https://doi.org/10.1029/2018JE005835

Wicht, J. (2002). Inner-core conductivity in numerical dynamo simulations. *Physics of the Earth and Planetary Interiors*, 132(4), 281-302. https://doi.org/10.1016/S0031-9201(02)00078-X

Wicht, J., & Christensen, U. R. (2024). Contributions of Jupiter's deep-reaching surface winds to magnetic field structure and secular variation. *Journal of Geophysical Research: Planets*, 129(4), e2023JE007890. https://doi.org/10.1029/2023JE007890

Wicht, J., & Heyner, D. (2014). Mercury's magnetic field in the MESSENGER era. In J. Shuanggen (Ed.), *Planetary geodesy and remote sensing*. CRC Press.
Wulff, P. N., Christensen, U. R., Dietrich, W., & Wicht, J. (2024). The effects of a stably stratified region with radially varying electrical conductivity on the formation of zonal winds on gas planets. *Journal of Geophysical Research: Planets*, 129(2), e2023JE008042. https://doi.org/10.1029/2023JE008042

Wulff, P. N., Dietrich, W., Christensen, U. R., & Wicht, J. (2022). Zonal winds in the gas planets driven by convection above a stably stratified layer. *Monthly Notices of the Royal Astronomical Society*, 517(4), 5584–5593. https://doi.org/10.1093/mnras/stac3045

Zuber, M. T., Smith, D. E., Phillips, R. J., Solomon, S. C., Neumann, G. A., Hauck, S. A., et al. (2012). Topography of the northern hemisphere of Mercury from MESSENGER laser altimetry. *Science*, 336(6078), 217–220. https://doi.org/10.1126/science.1218805