Mercury’s anomalous magnetic field caused by a symmetry-breaking self-regulating dynamo

The discovery of Mercury’s unusually axisymmetric, anomalously axially offset dipolar magnetic field reveals a new regime of planetary magnetic fields. The cause of the offset dipole remains to be resolved, although some exotic models have been proposed. Deciphering why Mercury has such an anomalous field is crucial not only for understanding the internal dynamics, evolutionary history and origin of the planet, but also for establishing the general dynamo theory. Here we present numerical dynamo models, where core convection is driven as thermo-compositional, double-diffusive convection surrounded by a thermally stably stratified layer. We show that the present models produce magnetic fields similar in morphology and strength to that of Mercury. The dynamo-generated fields act on the flow to force interaction between equatorially symmetric and antisymmetric components that results in north-south asymmetric helicity. This symmetry-breaking magnetic feedback causes the flow to generate and maintain Mercury’s axially offset dipolar field. A new regime of planetary magnetic fields was revealed through the MESSENGER spacecraft mission to Mercury. Here, the authors present a numerical dynamo model that can re-produce both the axisymmetric and anomalously axially offset dipolar magnetic field of Mercury.

I n-orbit observations by the MErcury Surface, Space ENvironment, GEochemistry and Ranging (MESSENGER) spacecraft confirms that Mercury currently possesses a global magnetic field generated by convective motions in the liquid iron core through dynamo processes [1][2][3] . Mercury's field intensity is about 1% of the Earth's field intensity at the surface (~700 nT), and has shown very weak secular variation over the past 40 years 4 . Moreover, MESSENGER has revealed the unusual morphology of Mercury's magnetic field, which is unlike that of any other planetary magnetic field: strongly axisymmetric dipolar fields with a dipole tilt angle <0.8°and the magnetic equator displaced northward by 0.2 R H , where R H = 2440 km is the Hermean radius [1][2][3] . Note that the MESSENGER measurements of the magnetic field were strongly biased towards the northern hemisphere, and that a smaller displacement of 0.14 R H was reported in an analysis with a larger dataset 5 .
Based on the co-density approach of treating thermal and compositional convection together 6 , the weak field and its low secular variation can be explained by dynamo models incorporating a thermally stably stratified layer beneath the core mantle boundary (CMB), through which small-scale, high-frequency components deep inside the core are attenuated due to the skin effect 7,8 . However, previous models have hitherto had difficulty reproducing the unusual morphology without rather speculative CMB boundary conditions 9,10 .
There are two key issues limiting these earlier findings: the codensity formulation, instead of the treatment of the flow as double-diffusive convection, was used in the anticipation of turbulent diffusivity, and this formulation may be invalid in a thick stably stratified layer 11 ; and the unique core crystallization due to the pressure-temperature condition of small bodies such as Mercury was not taken into account, but could result in compositional convection in a non-Earth-like manner, depending on the unknown core sulfur concentration 12,13 .
Three representative mechanisms have been proposed as potential drivers of the compositional convection in Mercury's core: a bottom-up (BU), top-down (TD), and snow-layer (SL) mode 14 . The BU mode is powered by either ejection of an Earthlike light element from the inner core boundary (ICB) or a Ganymede-like floatation of light FeS solid from an Fe-FeS alloy on the FeS-rich side of the eutectic 12 . The TD mode corresponds to precipitation of an Fe-rich solid as iron snow from the CMB or the bottom of the stably stratified layer 15 . The SL mode represents the case where iron snow occurs at a certain depth of the core 16 . The combined effect on Mercury's dynamo of double-diffusive convection and a core crystallization regime has not yet been explored.
Here, we present numerical dynamo models driven by thermo-compositional, double-diffusive convection in a rotating spherical shell that reproduce all the characteristic features of Mercury's magnetic field (see Methods). As in our previous study 17 , the diffusivity contrast between thermal and compositional diffusivities is 10, which is rather small compared with those expected in planetary cores. Even such a small difference could yield the magnetic fields different from those in the codensity model 11,17 . In the present model, a thermally stably stratified layer is imposed in the upper half of the core (Fig. 1). The heat/compositional flux is assumed to be fixed on the ICB and CMB. In particular, the zero-heat-flux condition is applied to the ICB to minimize the effects of bottom heating or maximize those of internal heating to drive thermal convection 9 . In all, we performed 12 runs for the BU, TD, and SL modes, where thermal and compositional driving forces (Rayleigh numbers) were varied, while other parameters were fixed. The parameters used in the present study are described in detail in the Supplementary Table 1.

Results
Comparison with observations. The radial component of the model magnetic field at the planetary surface is dominantly axisymmetric, dipolar, and asymmetric about the equator, which is consistent with the observations as well as with the strength of the field (Fig. 2a, b). The average dipole offsets from the center are 0.14 R H and 0.2 R H in the cases of BU1 and BU2, respectively (Supplementary Table 1); these values are comparable with the observed ones 3,5 . In order to explain the dipole offset, the Hermean-centered quadrupole component, amounting to 40% of the Hermean-centered axial dipole, is required in the case of BU2 and Mercury's field 3 , whereas the quadrupole component is 28% of the dipole in BU1 5 . On the other hand, the TD and SL models result in field morphologies dissimilar to the observed onesnamely, an unacceptably strong and equatorially antisymmetric morphology (Fig. 2c) and an overly weak and dominantly nondipolar morphology (Fig. 2d). These features are also confirmed by the magnetic power spectrum in terms of spherical harmonic degree and order (Fig. 2e, f). Moreover, the mean dipole tilt angle from the spin axis is less than 1°, and slow magnetic field variation like that of Mercury 2,4 is produced in BU1 (Supplementary Fig. 7 and Supplementary Table 1).
Formation of the offset dipolar field. The dipole offset results from the biased dynamo process in the core of the northern hemisphere ( Fig. 3a), which is prompted by the equatorially asymmetric flow structure ( Fig. 3b-d). The asymmetric structures in the convection vortices and magnetic field are temporally stable. The hemispherically averaged relative axial helicity (RAH) of convection is an important quantity with respect to the resultant magnetic field morphology 18 . In BU1, the time-averaged |RAH| is 0.42 for the northern hemisphere and 0.34 for the southern hemisphere (Supplementary Table 1), and these values differ significantly. It has been suggested that the north-south asymmetry in RAH is caused by the antisymmetric mode of convection, the kinetic energy of which is merely~10% of the total 9 .
Analysis of the axial helicity partition indicates that the symmetric and antisymmetric flows interact to yield the hemispherical bias in RAH ( Supplementary Fig. 2a, b). Other dynamo models driven by TD and SL but failing to have dipole offset show an insignificant difference of RAH between hemispheres despite having a fraction of the antisymmetric flow components similar to the BU model, indicating a negligible interaction between different flow modes (Supplementary Fig. 2c-f). More importantly, we found in the present study that an interactive flow structure that yields dipole offset is maintained by the dynamo-generated magnetic field itself. By switching off ri = 0.2 ro = 1 rs = 0.5 Fig. 1 Model geometry in units of core radius. The radius of the solid inner core is r i = 0.2. The interface of the thermally stably stratified layer is set at r s = 0.5. Above r s , the thermally stably stratified region extends to the core mantle boundary, r o = 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-018-08213-7 the feedback effect of the magnetic field on convection in BU1K (Helicity analysis in Methods), a nearly perfectly north-south symmetric flow structure ensues (Fig. 4a), and the quadrupolar magnetic field is generated as a result of kinematic dynamo action (Fig. 4b). The antisymmetric flow component rapidly declines with time ( Supplementary Fig. 8), indicating that the interacting antisymmetric flow mode is driven by the magnetic field. Thus, we conclude that a Mercury-like offset dipolar field is generated by a dynamo, with symmetry-breaking magnetic feedback to yield helicity bias, which we call self-regulation (SR).

Discussion
The present models show that, without taking any speculative external forcing into account, a Mercury-like magnetic field could result from a spontaneous dynamo mechanism driven by thermocompositional convection, in which the compositional buoyancy is caused by the BU compositional process and thermal buoyancy mostly by the internal heating in the liquid core. In a previous study using double-diffusive convection, thermal convection is driven from the bottom (i.e., non-zero ICB heat flux due to latent heat release; Supplementary Fig. 4) with a larger inner core, resulting in stronger, multipolar dynamos 11 . According to runs of BU1-4 and those in a previous study 11 , Mercury-like magnetic morphology is found in cases of zero ICB heat flux with a small inner core and modest flow vigor in terms of a magnetic Reynolds number of~100 (Supplementary Table 1), which is close to the estimate for Mercury's core 9 .
Each model has a velocity field of clearly different structure ( Supplementary Fig. 5): convection confined within a convectively unstable region in the case of BU, fingering-type convection in a thermally stably stratified layer in TD 11 , and faintly layered convection in SL 16 . The radial distributions of the velocity and magnetic fields show the prominence of the axisymmetric toroidal field in BU1 around the stratification boundary (Supplementary Fig. 6). The magnetic field within the convective region is strong enough for SR to work, particularly near the stable-unstable stratification boundary, whereas the weak magnetic field at the surface is due to attenuation associated with the skin effects of the stably stratified layer.
Therefore, in terms of the mechanisms involved in the generation and maintenance of Mercury-like magnetic field morphology, the key findings of the present work are as follows: the axial helicity distribution biased to one hemisphere is due to non-linear interaction of the equatorially symmetric and antisymmetric flow modes; and such interaction is prompted by the Lorentz force powering the antisymmetric flow mode as an SR process. These findings are in contrast to a previous study which explained the helicity asymmetry in terms of linear superposition of the symmetric and antisymmetric modes 9 . The asymmetric flow components also alter convective heat/compositional flux profiles in an asymmetric manner biased to the northern hemisphere ( Supplementary Fig. 9a). It should be mentioned that in spite of the antisymmetric flow mode being much smaller than the symmetric one here, their symmetry-breaking interaction could play a decisive role in generating an asymmetric dynamo 19 . Taking into account the fact that the action of the Lorentz force enhancing the helicity tends to prefer dipolar fields over quadrupolar fields 20,21 , and that the quadrupolar dynamo is obtained in the kinematic run (Fig. 4b), our results suggest that a convection that prefers the quadrupolar field in the absence of SR may require the hemispherical magnetic morphology. If so, feedback of the large-scale strong azimuthal toroidal field may play a role in generating the dipole component 22,23 .
Therefore, we next examined this possibility. The axisymmetric components of the magnetic and velocity fields in Fig. 5 show distinctive differences in structure. In the case of BU1, the asymmetric meridional circulation and zonal flow indicate an SR of the magnetic field acting in the northern hemisphere (Fig. 5a), whereas the almost symmetric flow structures in TD1 and SL1 denote a negligible SR effect on convection (Fig. 5b, c). Compared to BU1K (Fig. 5d), it is clear that the SR alters the flow and magnetic field structures in BU1. Note also that the azimuthal (toroidal) magnetic field generation around the edge of the convection columns and stratification boundary at mid-latitudes in the northern hemisphere is well correlated with the generation of the axial dipole and axial quadrupole components (Supplementary Fig. 10a). In addition, the zonal flow structure is deformed in the northern hemisphere by the SR. Without SR, the symmetric flow structure would appear as in BU1K, and then the symmetric toroidal field would be regenerated only near the equator ( Supplementary Fig. 10d). Eventually, the axial dipole component would also disappear. Although the field generation in the northern polar region is also remarkable, SR effects are unclear there. On the other hand, in TD and SL, magnetic field generation occurs in an almost symmetric way, and therefore there seems to be no SR in operation to feed the magnetically driven antisymmetric flows and break equatorial symmetry. Thus, it is suggested that generation of a biased axisymmetric toroidal field around the stably stratified layer is an important factor if the SR is to work.
To test whether double-diffusive treatment instead of codensity is necessary for the present result, we compared the BU1 to the case of BU1C (Supplementary Note 1) corresponding to the co-density approach, where thermal and compositional diffusivities were set to be equal. The BU1C results in a dipolar magnetic field without significant dipole offset, whose strength is much larger than that of Mercury ( Supplementary Fig. 11 and Supplementary Table 1). As shown in a previous study 11 , dynamos by double-diffusive convection and co-density could be different. Here, it is also demonstrated that double diffusion is an important factor for a Mercury-like field, although the complex physical processes involved must be elucidated in future research.
Asymmetric hemispherical dynamos have also been examined with thermal convection in other bodies such as the Sun 24 , Mars [25][26][27] , and Ganymede 15 . In the case of Sun-like models, the magnetic field contains substantial non-axisymmetric components, and shows wavelike periodic reversals. For the Mars-like models, a north-south asymmetric, heterogeneous outer boundary condition is imposed to create an asymmetric convection structure and magnetic field. In most of these models, the complicated magnetic field morphology contains non-axisymmetric components. Under homogeneous outer thermal boundary conditions, the asymmetric dynamo can be generated by the equatorially antisymmetric, axisymmetric (EAA) flow 28 . Because the EAA mode occupying the kinetic energy comparable with that in the symmetric mode occurs spontaneously as a result of thermal instability, the SR effect is not required to form the asymmetric structure. Dynamo models focusing on Ganymede by TD in the co-density formulation, where no SR seems to manifest, show a regime diagram between the input parameters, stable layer thickness, and magnetic field morphology 15 . The models mentioned should be investigated to see what the most probable mechanism is to spontaneously form a Mercury-like magnetic field. In this regard, whether or not the present findings could be extended towards the parameter regime appropriate to the planetary core remains to be examined. These are issues to be studied by a broader range of parameter survey.
The buoyancy source distribution in our Mercury-like dynamos has implications for the chemical composition and evolution of Mercury's core. Since thermal convection due to internal heating is preferred to that due to bottom heating from the latent heat release upon inner core solidification, slow cooling of the core, retarded inner core growth, and resultantly small inner core are suggested 29,30 . A considerable amount of radioactive heat source, such as potassium and uranium, would then be required to drive thermal convection, maintain a stably stratified layer, and keep the inner core small 31,32 . In this circumstance, compositional convection driven by the BU process due to inner core growth could be modest, consistent with the moderate values of the magnetic Reynolds number in our models. However, in order to form a thermally stably stratified upper layer and a deep convectively unstable layer, there are some issues to be considered regarding the mutual compatibility of buoyancy sources in the BU models. In the case of homogeneously distributed internal heat sources, the heat flux varies in proportion to the radius. Then there is no crossover between the adiabatic and super-adiabatic heat flux, and the entire fluid core is either stable or unstable. In addition, the zero-heat flux condition at the ICB is not compatible with compositional flux. To justify the model assumptions, it could be plausibly argued that the ICB heat flux due to latent heat release is very close to the adiabatic one, or that the radioactive elements are concentrated only in the deep layer, although no obvious justifications are available for the model assumptions. Another possibility could be that the adiabatic temperature gradient varies super-linearly with the radius because of the strong pressure dependence of the thermal expansion coefficient, so that the adiabatic heat flux can exceed the actual heat flux in the upper parts of the core, but falls short of it in the deep parts. These arguments would need verification based on the thermodynamic properties of iron and structural models for Mercury 14,33 .
A thermal and magnetic evolution model with an Fe-Si core allowing formation of a thermally stably stratified layer yields a present-day inner core larger than 800 km, and strong multipolar magnetic fields 33 . These findings may suggest that the core contains some amount of silicon, as well as a low concentration of sulfur (<6%), which prevents TD/SL convection by iron snow 14,34,35 . If so, these facts would pose additional constraints on the thermal and magnetic evolution of Mercury 36 .
The present models could be tested by good-coverage observation in the next BepiColombo mission 37 . This would allow us to conduct a detailed comparison between the improved data of the internal magnetic field and the predictions derived from numerical dynamo models of Mercury.

Methods
Numerical modeling of a planetary dynamo. We consider the convective motion of an electrically conducting, incompressible Boussinesq fluid in a rotating spherical shell of inner and outer radii r i and r o . In most cases, the aspect ratio is χ = r i / r o = 0.2, while an Earth-like value is χ = 0.35. The spherical shell is rotating about the z-axis at an angular rotation rate Ω. The governing equations described in nondimensional form are the Navier-Stokes equation for the velocity field u and nonhydrostatic pressure P, induction equation for the magnetic field B, heat transport equation for the temperature T, and transport equation for the light element concentration C. Length is scaled with the shell thickness D = r o − r i . Time is scaled with the viscous diffusion time D 2 /ν, where ν is the viscosity. The velocity is scaled with ν/D. The magnetic field is scaled with (2ρμηΩ) 1/2 , where ρ is the density, μ is the magnetic permeability in vacuum, and η is the magnetic diffusivity. The temperature and concentration of light elements are scaled with h o T D and h i C D, where h o T is the reference CMB temperature gradient without stable stratification and h i C is the reference ICB compositional gradient. In some models, temperature and light element concentration are scaled using the reference ICB temperature gradient and CMB compositional gradient.
Non-dimensional parameters of the dynamo models are the thermal Prandtl number (Pr T = ν/κ T = 0.1, where κ T is the thermal diffusivity), the compositional Prandtl number (Pr C = ν/κ C = 1, where κ C is the compositional diffusivity), the magnetic Prandtl number (Pm = ν/η = 3), the Ekman number (E = ν/2ΩD 2 = 10 −4 ), the thermal Rayleigh number (Ra T = αgh o T D/2Ων, where α is the rate of thermal expansion, and g is the gravitational acceleration at CMB), and the compositional Rayleigh number (Ra C = βgh i C D/2Ων, where β is the rate of compositional expansion).
The non-dimensional reference temperature profile without a stably stratified layer dT o /dr is described as where r s is the position of the stratification boundary, d s is the thickness of the transition between convecting and stably stratified regions, and Γ 0 is the temperature gradient across the stratified layer. Here, we mostly use r s = 0.5r o , d s = 0.05r o and Γ 0 = 10. Similarly, the radial profile of light element concentration dC o /dr is given by where ε C = −3r i 2 /(r o 3 − r i 3 ) = −0.097 and a C = −r i 2 r o 3 /(r o 3 − r i 3 ) = −0.063. An iron snow layer is represented by superimposing a Gaussian profile onto the basic one as follows: The radial profiles of non-dimensional heat/compositional flux corresponding to the BU, TD, and SL models summarized in Supplementary Table 1 are shown in Supplementary Fig. 1. The present treatment enables us to examine the effects of different modes of compositional convection while preserving the mode of thermal convection, which is an advantage over the co-density models. These profiles are used to implement a stably stratified layer in a mathematically convenient way, because the radial profiles of temperature and composition in Mercury's core are still poorly constrained. Those used in a previous study 11 are shown in Supplementary Fig. 4 for comparison. At both boundaries, the boundary condition for the velocity field is no-slip and insulating for the magnetic field. The influence of treating the inner core as an insulator on the results is insignificant because of its small size 40 . At the outer boundary, we adopted a zero-flux condition for composition and fixed-flux condition for temperature, whereas we assume flux is fixed for composition and either fixed or zero for temperature at the inner boundary. The zero-heat-flux condition is given so as to minimize effects of bottom heating or maximize effects of volumetric internal heating to generate the asymmetric magnetic field 9 .
Initial conditions are given by either random perturbations of the temperature and composition, and an axial dipole field as a seed field or the final result of a run at different parameters. The numerical code used in this study is an extended version of refs. 17,41 . The spatial resolution is 100 or 128 grid points in the radial direction. A spherical harmonics expansion in the angular directions is used up to degree and order 127.
Axial dipole offset. Using the Gauss coefficients ðg M L ; h M L Þup to degree and order two, the eccentricity of the best-fitting dipole from the Hermean center (x 0 , y 0 , z 0 ) is given as 42 where Ignoring non-axisymmetric terms such as Mercury's field, we have x 0 ; y 0 ; z 0 ð Þ$ð 0; 0; Þ. This expression indicates that dipole offset in the axial direction is determined by the ratio of the Hermean-centered axial quadrupole to the centered axial dipole, and also that the offset direction is northward, if these two terms have the same polarity (and southward otherwise).
Diagnostic quantities. Important diagnostic quantities of dynamo simulations used in this study include the magnetic Reynolds number, Rm = Du/η, Elsasser number, Λ = B 2 /(2ρμηΩ), Gauss coefficient of the axial dipole, g 0 1 , axial dipole offset normalized by Hermean radius, D offs =z 0 /R H , dipole tilt angle, Tilt, dipolarity, F dip , fraction of the axisymmetric magnetic field, F axs , fraction of the equatorially antisymmetric flow components in the total kinetic energy, K asym , local Rossby number Ro l , and the hemispherically-averaged relative axial helicity in each hemisphere, RAH N/S , defined as where u z and ω z are the axial components of the velocity and vorticity. Volume integration is taken in terms of either the northern (N) or southern (S) hemisphere.
These expressions suggest that the difference of relative axial helicity in its absolute value is caused by the interaction of the symmetric and antisymmetric flow components when the correlation between the different modes of the velocity (u s z ; u A z ) and the vorticity (ω S z ; ω A z ) is not good. Supplementary Fig. 2 clearly shows that north-south asymmetric helicity distribution in the Mercury-like offset dipolar dynamo can be explained by the evident interaction terms ( Supplementary Fig. 2a, b), while the nearly perfect symmetric helicity distributions in the models without dipole offset are due to the velocity field with negligible interaction ( Supplementary Fig. 2c-f). We therefore conclude from this analysis that asymmetric helicity due to interaction between different flow modes leads to hemispherically biased dynamo action, and consequently a Mercury-like offset dipolar magnetic field is generated.
We then investigate what mechanisms are responsible for an interaction among flow components that maintained the biased helicity distribution. For this purpose, we carry out three additional runs building on the BU1, where simulations are restarted using the final result of BU1 with some modifications. The first run, BU1L, is restarted with the antisymmetric mode of the velocity, temperature, and composition fields reversed with respect to the equator so that RAH is artificially concentrated in the southern hemisphere rather than in the north. The second run, BU1K, is a kinematic dynamo run, where the Lorentz force term is dropped from the momentum equation. The third, BU1M, is a run of the magnetohydrodynamic dynamo resumed with the symmetric components of the magnetic field reversed in terms of the equator. The first run is designed to see effects of the linear and nonlinear hydrodynamic terms and the stably stratified layer, while the second and third runs are intended to show the non-linear effects of the dynamo-generated magnetic field on the core flow. Figure 4 in the main text and Supplementary Fig. 3 clearly show that the helicity bias is sustained by the Lorentz force due to the dynamo action itself; that is, the dipole offset is a natural result of the selfregulation process of the core dynamo prompting symmetry-breaking interaction.
Magnetic field generation. In order to investigate the mechanisms of magnetic field generation, we consider the equation for magnetic energy variation: ∂ ∂t The rightmost term of the equation represents the magnetic energy enhancement achieved by stretching the magnetic field lines due to flow gradient. The stretching terms responsible for generation of the axial dipole, B 0 P1 , axial quadrupole, B 0 P2 , and axisymmetric toroidal components, B 0 T , are examined by calculating Þ , respectively. Time-averaged results are shown in Supplementary Fig. 10.
Code availability. Numerical code for dynamo simulations is available upon request.