Influence of reference states on Jupiter’s dynamo simulations

1.   Introduction

It is widely accepted that planetary magnetic fields are generated in an electrically conducting fluid layer in their deep interiors through the dynamo process (Jones, 2011). Therefore, planetary magnetic fields provide a window to infer the internal structure and dynamics of planets (Stevenson, 1982a, 2003; Guillot, 2005; Schubert and Soderlund, 2011; Liu SH et al., 2020). Among our solar system’s planets, Jupiter exhibits the strongest magnetic field, which is thought to be generated in its metallic hydrogen region (Stevenson, 1982a; Guillot, 2005; Jones et al., 2011). Since 2016, NASA’s Juno spacecraft has been measuring Jupiter’s magnetic field and gravity field with unprecedented precision, creating data that can be used to infer the internal structure and dynamics of the planet (Bolton et al., 2017a, b).

Pre-Juno observations revealed that Jupiter bears an Earth-like dipole-dominated magnetic field, though the dipole moment of Jupiter is much larger than that of Earth (Connerney et al., 1998). Recent studies based on Juno’s observations have found that Jupiter’s magnetic field also has a significant non-dipolar part in the northern hemisphere (Moore et al., 2018). Numerical simulations of Jupiter’s dynamo action can produce Jupiter-like magnetic fields (Gastine et al., 2014; Jones, 2014; Duarte et al., 2018; Glatzmaier, 2018; Wicht et al., 2019; Tsang and Jones, 2020), though numerical models are intentionally simplified and use unrealistic control parameters due to the restrictions of computational power. Nevertheless, numerical modeling remains a useful approach to understand the physical processes and mechanisms of Jupiter’s magnetic field generation.

Jupiter’s dynamo is modeled as convective perturbations around a nearly adiabatic and hydrostatic reference state — conditions prescribed to make numerical simulation feasible (Jones, 2014). The reference state depends on physical properties and the equation of state of the planetary interior, which are not well constrained. Two types of fully convective reference state models are usually used in Jupiter’s dynamo simulations. One is based on ab initio calculations of Jupiter’s interior (Jones, 2014); the other type is based on the polytropic and ideal gas assumptions (Gastine et al., 2014). Different reference models result in different densities and temperature gradients within the planet, and thus lead to different buoyancy forces when other control parameters and conditions are the same. Therefore, the reference state prescribed can influence the operation of dynamo action and the generated magnetic fields. Given the uncertainty of the background state of Jupiter, it is important to characterize how sensitive the dynamo process is to the reference state model employed in each simulation.

Previous studies have used a variety of reference states in Jupiter’s dynamo simulations including both the ab initio model (Jones, 2014; Dietrich and Jones, 2018; Glatzmaier, 2018; Tsang and Jones, 2020) and polytropic models (Gastine and Wicht, 2012; Gastine et al., 2014; Duarte et al., 2018; Wicht et al., 2019). Gastine et al. (2012) studied the influence of the density stratification on dynamos and found that stronger density stratification leads to non-dipolar magnetic fields, but they used polytropic reference states in all simulations. Duarte et al. (2018) performed more than 60 simulations with varying reference states (including both ab initio models and polytropic models), electrical conductivity profiles, and control parameters to explore the sensitivity of the dynamo behavior to those factors. They found that the reference model tends to be not so critical to the dynamo action compared to the other internal properties, though detailed analysis and comparisons are required.

In this study, we investigate how choice of the reference density and temperature profiles influence the inferred dynamos and generated magnetic fields. We carry out two groups of convection-driven dynamo based on the anelastic approximation. Within each group, we test different reference models of both ab initio and polytropic types but keep all other control parameters and conditions identical. It turns out that the influence of the reference state on the dynamo action depends on the regime of force balance. When the buoyancy is dominant over the Lorentz force, different reference states can lead to very different dynamical outcomes. However, when the Lorentz force is dominant over the buoyancy, the dynamo action is insensitive to the reference state.

The remaining part of the paper is organized as follows. Section 2 introduces the numerical model of dynamo simulations and the reference state. Section 3 presents the numerical results of the two groups. The paper is closed with a discussion and conclusion in Section 4.

2.   Numerical Model

2.1   Governing Equations

We model Jupiter’s dynamo region as a rotating spherical shell with inner radius ${r}_{i}$ and outer radius ${r}_{o}$ (i.e., the radius ratio $\eta ={r}_{i}/{r}_{o}$), which rotates at $\boldsymbol{\varOmega} =\varOmega {\boldsymbol e}_{z}$. The spherical shell is filled with fluid of homogenous viscosity $\nu$ and thermal diffusivity $\kappa$, but radially varied electrical conductivity $\sigma \left(r\right),$ which mimics the metallic hydrogen region of Jupiter. The magnetic diffusivity $\lambda \left(r\right)= 1/\left(\mu_0\sigma \right)$ is thus also a function of the radius, where ${\mu }_{0}$ is the vacuum magnetic permeability. The dynamo action is driven by thermal convection under the anelastic approximation, which considers the density $\overline{\rho }\left(r\right)$ and temperature $\overline{T}\left(r\right)$ variation along the radius but filters out sound waves to avoid the smaller time step (Braginsky and Roberts, 1995; Lantz and Fan Y, 1999). The convection-driven dynamo process is described as small perturbations around an adiabatic and hydrostatic reference state that is spherically symmetric, i.e., depending on the radius only.

The dimensional equation of motion is given as (Jones et al., 2011),

where $\boldsymbol{u}$ is the velocity, ω is the vorticity, $\boldsymbol{B}$ is the magnetic field, $p{'}$ is the pressure perturbation, $\boldsymbol{j}=\dfrac{1}{{\mu }_{0}}\nabla \times\boldsymbol{B}$ is the current density and $S$ is the rate-of-strain tensor given by:

The mass conservation under the anelastic approximation is given as

Conservation of energy gives

where $s$ is the specific entropy, and κ is the entropy diffusivity, ${Q}_{v}$ and ${Q}_{j}$ are viscous heating and Joule heating, respectively.

The magnetic induction equation reads

and the magnetic field is divergence-free:

The above equations are solved in a dimensionless form. Using the shell thickness $d={r}_{o}-{r}_{i}$ as the length scale, the viscous diffusion time ${\tau }_{\nu }={d}^{2}/\nu$ as the time scale, $v/d$ as the velocity scale and the magnetic unit (ρ₀ μ₀ λ_i Ω)^1/2, the dimensionless equations governing the anelastic convection dynamo are given as (Jones et al., 2011; Duarte et al., 2018).

The dimensionless density $\widetilde {\rho }=\overline{\rho }\left(r\right)/{\rho }_{0}$, temperature $\widetilde {T}=\overline{T}\left(r\right)/{T}_{0}$, and magnetic diffusivity $\widetilde {\lambda }=\lambda \left(r\right)/{\lambda }_{i}=\sigma \left(r\right)/{\sigma }_{i}=1/\widetilde {\sigma }$ represent the dimensionless reference state profiles. Note that the density and temperature are normalized by the dimensional values at the outer boundary while the magnetic diffusivity is normalized by the dimensional magnetic diffusivity at the inner boundary ${\lambda }_{i}$. The detailed reference states will be presented in Section 2.3.

From now on, all physical quantities are given in dimensionless form unless otherwise stated.

Based on the aforementioned normalization, the dynamo model is controlled by the following dimensionless parameters, i.e., Rayleigh number $Ra$, Ekman number $E$, Prandtl number $Pr$, magnetic Prandtl number $Pm$, the dissipation number:

where ${\alpha }_{o}$ is the thermal expansion coefficient at the outer boundary, $\Delta s$ is the entropy difference between the inner and outer boundaries and ${g}_{o}$ is the gravity at the outer boundary.

The outer boundary is set as stress-free and insulating, while the inner core is regarded as rigid and finitely conductive; constant entropy is applied on both boundaries.

2.2   Reference States

As mentioned before, the dynamo process is modeled as perturbations around a hydrostatic reference state, which remains poorly constrained. The main objective of this paper is to investigate the influence of different reference states on the dynamo process. Previous studies usually used two different types of reference state models (Gastine et al., 2014; Jones, 2014; Duarte et al., 2018; Wicht et al., 2019; Tsang and Jones, 2020). The first type of reference state is the so-called ab initio model, which is based on ab initio simulation of Jupiter’s interior. The reference density $\widetilde {\rho }$ and temperature $\widetilde {T}$ profiles are given by polynomials fitted to the ab initio simulations results of French et al. (2012). The gravity is self-consistently determined by the density profile, and the buoyancy force is directly related to the temperature gradient

where $\widetilde {\alpha }=\alpha /{\alpha }_{o}$ is the normalized thermal expansivity, $\widetilde {g}=g/{g}_{o}$ is the normalized gravity.

The second type of reference state is based on the polytropic model, i.e., $\widetilde {p}={\widetilde {\rho }}^{1+1/m},$ and the ideal gas assumption $\widetilde {p}=\widetilde {\rho }\widetilde {T}$, leading to $\widetilde {\rho }={\widetilde {T}}^{\,m}$, where $m=1/\left({c}_{p}/{c}_{\nu }-1\right)$ is the polytropic index. In order to determine the density and temperature profiles, we need to make use of the relation between the temperature gradient and the gravity given in Equation (12) and prescribe a gravity profile. Following Duarte et al. (2018), we assume the gravity is proportional to the radius $r$, which suggests that the temperature gradient also is proportional to the radius, given $\widetilde {\alpha }=1/\widetilde {T}$ for the ideal gas. As noted by Duarte et al. (2018), such a gravity profile is unrealistic for Jupiter. Nevertheless, it serves as a simplified model to investigate the influence of the reference state model on the dynamo process. Based on these assumptions, the reference density, and temperature profiles are then uniquely determined by the polytropic index m and the number of density scale heights ${N}_{\rho }=\rm{ln}\left[\dfrac{\widetilde {\rho }\left({r}_{i}\right)}{\widetilde {\rho }\left({r}_{o}\right)}\right]$ in the layer.

In this study, we perform two groups of dynamo simulations. In each group, we consider different reference states while keeping other control parameters and conditions the same. The reference state profiles of $\widetilde {T}\left(r\right),\widetilde {\rho }\left(r\right)$ and $\rm{d}\widetilde {T}\left(r\right)/\rm{d}r$ used in Group A and Group B are shown in Figure 1a and 1b, respectively.

Figure  1.  Radial profiles of normalized temperature, density, and buoyancy term. Cyan line represents polynomial fitting to ab initio simulation results of French et al. (2012). (a) The simulations in Group A are cut at $0.99{R}_{J}$; $\eta =0.2$ is the ratio of inner core radius to outer boundary radius. (b) In Group B, the cut radius is ${0.9586R}_{J}$ and the radius ratio $\eta =0.0963$.

DownLoad: Full-Size Img PowerPoint

For the electrical conductivity, we use profiles given by continuous equations (Gómez-Pérez et al., 2010):

where ${r}_{m}$ is the transition radius, $a$ is the decay rate, ${\sigma }_{m}$ is the conductivity at the transition point. Figure 2 shows the electrical conductivity profiles we used in Group A (blue curve) and Group B (red curve). Detailed parameters of the reference state models are also given in Table 1.

Figure  2.  Normalized electrical conductivity profile ${\tilde \sigma} (r)$. The dotted line represents the conductivity results from the ab initio simulation of French et al. (2012). The blue (red) line represents the conductivity profile used in Group A (Group B), which is determined by polynomial Equation (13) with controlling parameters $a=13.0,{\sigma }_{m}=0.2,{r}_{m}=90$ (Duarte et al., 2018).

DownLoad: Full-Size Img PowerPoint

Table  1.  Input dimensionless parameters and time-averaged results. The polytropic reference models are controlled by the number of density scale heights ${N}_{\rho }$ and polytropic index m. ${R}_{m}$ and $\varLambda$ are time-averaged magnetic Reynolds number and Elsasser number.

	Model	$N_\rho$	m (polind)	η	Initial B field	Rm	Λ	Dynamo regime
Group ARa = 4.78e+07; E = 1.0e−04; Pm = 2.0; Pr = 1.0	A1	–	–	0.2	Weak dipole	4.30E+01	–	No dynamo
	A2	4	2	0.2	Weak dipole	2.77E+02	1.54E+00	Multipolar dynamo
	A3	5	2	0.2	Weak dipole	1.31E+02	6.21E-01	Multipolar dynamo
Group B Ra = 3.60e+07; E = 2.5e−05; Pm = 6.2; Pr = 0.1	B1	–	–	0.0963	Strong dipole	7.55E+02	2.16E+01	Strong dipole
	B2	2.95	1.879	0.0963	Strong dipole	1.27E+03	1.47E+01	Strong dipole

 | Show Table

DownLoad: CSV

2.3   Numerical Method

The fully nonlinear MHD Equations (1)–(3) are solved using an open-source MHD code, MagIC (https://magic-sph.github.io/), which employs the pseudo-spectral method based on an expansion of spherical harmonics in angular direction, and Chebyshev polynomials in the radial direction (Gastine et al., 2012). The code makes use of the fast spherical harmonic library SHTns (https://nschaeff.bitbucket.io/shtns/) (Schaeffer, 2013).

To avoid the time step getting too small, the simulations of Group A are cut at $0.99{R}_{J}$. The three simulations of Group A begin with a weak magnetic field. We use 73 radial grid points and 320 azimuthal grid points in all Group A simulations. The two models in Group B use similar parameters as in the case A of Jones (2014), in which the cut-off radius is $0.958{R}_{J}$. We use the 161 radial grid points and 384 azimuthal grid points in both Group B simulations. The convergence of the numerical simulations was carefully checked by examining the energy spectra of the simulations and by increasing the resolution for some simulations.

2.4   Diagnostic Parameters

The kinetic energy and magnetic energy are defined as the volume integral over spherical shell (Jones, 2011).

Three dimensionless diagnostic parameters are used here to qualitatively characterize the properties of the solutions: the Reynolds number ${R}{e}$, the magnetic Reynolds number ${R}{m}$, the Elsasser number $\varLambda$:

where the ${u}_{\rm{rms}}$ and ${B}_{\rm{rms}}$ are the RMS values of velocity and magnetic field ${\lambda }_{i}$.

3.   Numerical Results

We perform two groups of simulations to explore the influence of reference states. In each group, we use the same control parameters, electrical conductivity profiles, and initial conditions, but different reference states of the density and temperature are considered. Detailed information of all simulations is summarized in Table 1.

3.1   Energy

Figure 3 shows the time evolution of the magnetic energy and kinetic energy of all simulations, which gives an overview of our numerical results. In Group A, dynamo action is sustained for the polytropic reference states (A2 and A3) while the case using the ab initio reference state (A1) fails to sustain a dynamo (Figure 3a). The saturated kinetic energy of the ab initio case is also much smaller than that of the polytropic cases, as we can see from Figure 3b. These differences are essentially due to the effective buoyancy force being weaker in A1 compared to A2 and A3, as one can see from Figure 1a, despite the nominal Ra being the same in all cases.

In Group B, both the polytropic and ab initio reference states (B1 and B2) lead to sustained dynamo action, with comparable magnetic energy (Figure 3c), although the kinetic energy in B1 (ab initio) is slightly smaller than that in B2 (polytropic). Note that the effective buoyancy forces are quite different between B1 and B2, as one can see from Figure 1b.

Figure  3.  Time series of energy of the Group A (upper panel) and Group B (lower panel). Time is in units of the viscous time scale.

DownLoad: Full-Size Img PowerPoint

As we can see from the above results, the influence of different reference states on dynamo action can be either significant (Group A) or negligible (Group B). It turns out that the different behavior in Groups 1 and 2 may be explained by the different regimes of the force balance, as we will discuss in Section 3.3.

3.2   Morphology of Flows and Magnetic Fields

Figure 4 shows the radial velocity in the equatorial plane and the zonal velocity (i.e., axisymmetric ${u}_{\varphi }$) in the meridional plane of the three models in Group A. Convective motions in A1 take place mainly in the outer envelope and are very weak in the deep region (Figure 4a) due to the effective buoyancy force associated with the ab initio reference state. Indeed, such convective flows ($R_m=43$) are not sufficient to sustain a dynamo, as we have shown in Figure 4a. One can see more vigorous convection taking place in A2 and A3, which use the polytropic reference states, though the two models exhibit slightly different flow structures. In A3, with stronger density stratification, convection flow with small scale and high amplitude occurs primarily in the outer regions, because of the increasing buoyancy along the radius (Gastine and Wicht, 2012).

Figure  4.  Snapshots of radial velocity u_r (in dimensionless units of the Reynolds number) in the equatorial plane (left panel) and zonal flow u_φ in the meridional plane (right panel) of three models in Group A. (a−b): A1; (c−d): A2; (e−f): A3.

DownLoad: Full-Size Img PowerPoint

The polytropic models can produce very similar temperature and density profiles as the ab initio model at some parameters, but their buoyancy force term is generally larger because of their larger ‘artificial’ temperature gradient. Therefore, stronger convection and fluid velocities tend to be generated in the polytropic models.

In all cases, mean zonal flows are developed due to the nonlinear interactions of convective flows (Zhang K, 1992; Christensen, 2002; Heimpel et al., 2005; Jones and Kuzanyan, 2009). The zonal flow in A1 is nearly geostrophic (Figure 4b) without the influence of the Lorentz force, whereas zonal flows in A2 and A3 are influenced by the Lorentz force, leading to less geostrophic zonal flows (Figure 4d and 4f).

Figure 5 shows the radial velocity in the equatorial plane and the zonal velocity in the meridional plane of two models (B1 and B2) in Group B. In this group, the two models exhibit similar flow structures with vigorous convection in the whole domain and strong prograde zonal flows near the equator, although different reference states are used. The effective buoyancy forces resulting from the ab initio reference state (B1) and the polytropic reference state (B2) are rather different (see Figure 1b), but both models produce sufficiently vigorous fluid motions ($R_m=755$ and $1270$ respectively) to drive dynamos and sustain strong magnetic fields ($\varLambda =21.6$ and 14.7 respectively). As we will show in section 3.3, the Lorentz force plays a dominant role in both models in Group B. Therefore, fluid motions are strongly influenced by the Lorentz force, but are less sensitive to the reference state that determines the effective buoyancy force.

Figure  5.  Snapshots of radial velocity u_r (in dimensionless units of the Reynolds number) in the equatorial plane (left panel) and zonal flow u_φ in the meridional plane (right panel) of three models in Group B. (a–b): B1; (c–d): B2.

DownLoad: Full-Size Img PowerPoint

Figure 6 shows the radial magnetic field at the outer surface and in a meridional plane for the two cases with successful dynamos in Group A. We can see that A2 and A3 produce dynamos with a multipolar feature. Figure 7 shows the radial magnetic field for the two models in Group B, both of which sustain similar dipole-dominant magnetic fields at the outer surface.

Figure  6.  Snapshots of the radial magnetic field B_r (in dimensionless units of the Elsasser number) at the outer surface (left panel) and the azimuthal plane (right panel) for A2 (a–b) and A3 (c–d) in Group A.

DownLoad: Full-Size Img PowerPoint

Figure  7.  Snapshots of the radial magnetic field B_r (in dimensionless units of the Elsasser number) at the outer surface (left panel) and the azimuthal plane (right panel) for B1 (a–b) and B2 (c–d) in Group B.

DownLoad: Full-Size Img PowerPoint

3.3   Force Balance

As we have shown, the influence of the reference states can be very different for simulations in Group A and Group B. In this section, we try to understand the different consequences of the reference state from the point of view of force balance. In doing so, the RMS of each force term F in the momentum equation is calculated and decomposed into a sum of spherical harmonic contributions ${F}_{{l}}$ (Aubert et al., 2017; Schwaiger et al., 2019):

Figure 8 shows the time-averaged RMS force terms as a function of the spherical harmonic degree ${{l}}$ for inertial force (${F}_{I}$), Coriolis force (${F}_{C}$), Lorentz force (${F}_{L}$), buoyancy force (${F}_{B}$) pressure gradient (${F}_{P}$). The Inertial force and viscous force are not shown here due to their subdominant roles at smaller ${{l}}$ (larger scales).

Figure  8.  Force balance spectrum plotted as a function of spherical harmonic degree ${{l}}$. The mean values of different forces are represented by solid lines; corresponding one-standard deviations across the mean are denoted by the shaded area. (a): A2; (b): A3; (c): B1 (d): B2.

DownLoad: Full-Size Img PowerPoint

In all simulations, the leading order force balance at a large scale is the so-called quasi-geostrophic balance, i.e., a balance between the Coriolis force and the pressure gradient, suggesting that rotation plays a dominant role. The major difference between Group A and Group B is the relative contribution of the Lorentz force and the buoyancy force. In Group A, the buoyancy force is dominant over the Lorentz force for ${{l}}<20$ (A2) and ${{l}}<30$ (A3). As A1 fails to sustain a dynamo, the force balance of this case is not shown here. In contrast to Group A, the Lorentz force dominates over the buoyancy force at almost at all scales in both B1 and B2 in Group B.

These results suggest that the influence of the reference state depends on the regime of force balance. In Group A, the buoyancy force plays a dominant role apart from the leading order quasi-geostrophic balance. Therefore, the contrast of the effective buoyancy caused by the ab initio reference state and polytropic reference states lead to very different dynamical outcomes in Group A. When the Lorentz force is more dominant than the buoyancy force, different reference states lead to rather different buoyancy forces, but the dynamics are mainly controlled by Lorentz force apart from the leading order quasi-geostrophic balance. Therefore, the reference states tend to be much less influential in the dynamical outcomes. This is broadly in line with theoretical expectations that the Lorentz force determines the flow structures and convective length scales (Zhang KK and Schubert, 2000), akin to magnetoconvection. However, we should note that, in order to reach this regime, the buoyancy force needs to be strong enough to sustain a strong field dynamo. Our results in Group B are in line with previous studies (Gastine and Wicht, 2012; Jones, 2014; Duarte et al., 2018), which found that dynamo action is not sensitive to the reference density and temperature profiles used in the simulations.

4.   Discussion and Conclusion

This paper has investigated the influence of different reference state models on Jupiter-like dynamo simulations. We performed two groups of dynamo simulations using different reference state models. In each group, we considered both the ab initio and the polytropic types of reference model for the background density and temperature while we kept other control parameters and conditions the same.

We found that the influence of the reference state can be either substantial (models in Group A) or minor (models in Group B), depending on the regime of the force balance. Different reference profiles of the density and temperature lead to a rather different buoyancy force despite the same Ra being used. When the buoyancy is dominant over the Lorentz force, varying the reference state leads to very different convective motions and dynamo processes. However, when a strong field dynamo is sustained and the Lorentz force becomes dominant, convective motions and dynamo actions become less sensitive to the buoyancy force. Therefore, both convective motions and generated magnetic fields are similar despite rather different reference states being used in a strong field regime.

Our numerical results may have significant implications for dynamo simulations of Jupiter and other giant planets. The internal states and structures of giant planets remain poorly constrained (Stevenson, 2020). Jupiter’s dynamo is expected to be in a strong field regime (Aurnou and King, 2017), i.e., the Lorentz force plays a dominant role, so our numerical results suggest that the uncertainty of the internal structure has little influence on the planet’s dynamo process. However, we note that there are several caveats in this study that should be considered in the future. Firstly, we have not considered the influence of the electrical conductivity profiles when changing the reference density and temperature profiles. Indeed, previous studies (Dietrich and Jones, 2018; Duarte et al., 2018; Wicht et al., 2019) have found that the conductivity profiles can be very influential on the dynamo action. Secondly, although different density and temperature profiles are employed in this study, we assume that the whole simulation domain is convectively unstable. Stably stratified layers may exist within giant planets; if so, they could change the magnetohydrodynamics significantly (Stevenson, 1982b; Cao and Stevenson, 2017; Christensen et al., 2020). Finally, our numerical simulations are limited to a small range of control parameters with moderate Ekman numbers. Further simulations over a wide range of control parameters at lower Ekman number are required, though such simulations will be computationally demanding.

Acknowledgments

We would like to thank Dali Kong and Hao Cao for very useful discussions and suggestions. This study was supported by the B-type Strategic Priority Program of the CAS (XDB41000000) and the pre-research project on Civil Aerospace Technologies of CNSA (D020308) and the Macau Foundation. CAJ acknowledges support from STFC, grant number ST/S00047X/1 held at the University of Leeds. Simulations were performed on the Taiyi cluster supported by the Center for Computational Science and Engineering of Southern University of Science and Technology. This study made use of the open-source code MagIC, which is freely available at https://magic-sph.github.io.