© The Authors 2025

1 Introduction

During the past decades, many numerical solutions to the magnetohydrodynamic (MHD) equations in convective spherical shells have been obtained under various assumptions and parameters. Simulations of this type produce an amplification and self-sustaining of magnetic fields that can recover some observed aspects related to a planetary dynamo in gas giants, such as the magnetic field morphology or the latitudinal variations in the atmospheric jets (Jones 2011; Schubert & Soderlund 2011, e.g.,). High-resolution models are getting close to reproducing the internal dynamos of Earth and Jupiter (Schaeffer et al. 2017; Gastine & Wicht 2021, respectively), including stochastic dipole reversals that resemble the geomagnetic field (Glatzmaier & Coe 2007; Meduri et al. 2021). In planetary science, MHD equations are typically expressed using dimensionless dynamo numbers: Rayleigh, Ekman, Prandtl, and magnetic Prandtl numbers, or combinations thereof. These numbers quantify the relative influence of various fluid forces such as the dissipative, buoyant, and Coriolis forces. There is an unavoidable computational caveat: the parameter space accessible through numerical simulations diverges significantly from the physical reality, and some parameters (Rayleigh and Ekman in particular) diverge by many orders of magnitude. The reason is that the range of the relevant spatial scales that are to be followed is too wide, that is, from the microscopical diffusion to the planetary scale for the global rotation or convection patterns.

To partially overcome this intrinsic drawback, some studies have used many numerical models to find scaling laws between the different dimensionless numbers characterizing the dynamo. For example, for rapidly rotating dipole-dominated dynamo solutions under the Boussinesq approximation, Christensen & Aubert (2006) derived scaling laws that connect the dynamo parameters spanning at least two orders of magnitude. These relations should also arguably work in the real planetary regime, as the relative importance of each term in the Navier-Stokes equation (the force balance) is expected to be similar to numerical models (Davidson 2013; Yadav et al. 2016).

For the modeling of gas giants, the anelastic approximation (e.g., Braginsky & Roberts 1995; Gilman & Glatzmaier 1981; Glatzmaier 1984, 1985a,b; Lantz & Fan 1999) is more appropriate than the Boussinesq approximation, as it allows for density variations but still effectively filters out the sound and magnetosonic waves. It relies on using a static, adiabatic, and spherically symmetric background reference state that is specified by density, gravity, temperature, and other thermodynamic variables. In addition, the velocity and magnetic fields are evolved together with the deviations from the background. These equations have been extensively used to model the magnetic field of gas giants and stars. Usually, the no-slip boundary conditions employed by geodynamo models are replaced by stress-free conditions because they reflect the nature of the outer atmosphere of the gas giants better. The simulations are more computationally demanding, and the numerical solutions have several features that are not seen in the Boussinesq approximation. For example, the development of Jupiter-like zonal flows, which promote weaker multipolar fields, and strong dipole fields, which can suppress these zonal flows via Lorentz forces, often compete in these models (e.g., Grote et al. 2000; Simitev & Busse 2003, 2009; Sasaki et al. 2011; Schrinner et al. 2012; Duarte et al. 2018). Another example is the bistability found for not too large Rayleigh numbers: Dipolar-dominated and multipolar solutions coexist with identical parameters (Gastine et al. 2012), where different solutions can be reached by setting different initial conditions (Schrinner et al. 2012). In this context, Yadav et al. (2013) provided scaling laws for dynamo models under the anelastic approximation similar to those of Christensen & Aubert (2006). In this case, dipolar and multipolar solutions as well as a different range of density stratification and radial-dependent diffusivities were used.

For a more realistic modeling of gas giants, the challenge is to incorporate the outer steep gradients of various thermodynamical profiles, since they imply very different timescales farther outward. Moreover, there is a steep outward decrease in the electrical conductivity because hydrogen is not in the metallic state in the outer planetary layers. This behavior was quantified along the Jovian adiabat by French et al. (2012) and was used during the past decade in different dynamo simulations (Gastine & Wicht 2012; Jones 2014; Gastine et al. 2014; Wicht et al. 2019b). These studies provided important results in terms of a comparison with the data provided by the ongoing Juno mission and earlier Jovian missions, and they highlighted the importance of incorporating a realistic background, although with the caveat on the nondimensional dynamo numbers mentioned above. In recent years, some works (Gastine & Wicht 2021; Yadav et al. 2022; Moore et al. 2022) have added a stably stratified layer just below the region in which metallic hydrogen starts to mimic a helium-rain region due to hydrogen-helium demixing (Klepeis et al. 1991; Nettelmann 2015; Nettelmann et al. 2015). This layer helps to naturally obtain alternating east–west zonal winds centered around the equator for the Jovian dynamo, as well as a highly axisymmetric magnetic field for the Saturn model.

By using the aforementioned scaling laws as well as observations, Christensen et al. (2009) determined that for both planets and fast-rotating stars, the energy flux determines the magnetic field strength. To match observational constraints, Reiners et al. (2009) expressed this law in terms of the mass M, luminosity L, and radius R in a more simplified form. Using this scaling law and the analytical evolutionary tracks for substellar objects in Burrows & Liebert (1993); Burrows et al. (2001), Reiners & Christensen (2010) provided a scenario for the evolution of the magnetic field, with which they obtained a steadily weakening (a factor of ≈10 over around 10 Gyr) of the magnetic field at the dynamo surface.

We aim to address the long-term evolution of the dynamo action in Jupiter-like planets through an alternative approach. We perform 3D anelastic dynamo simulations with a background corresponding to different ages of the long-term planetary evolution. By comparing how the solutions change from one age to the next and recalling the intrinsic caveats related to the accessible ranges of nondimensional numbers, we evaluate the trend in morphology and intensity changes that an internal dynamo in cold gas giants can undergo during its evolution over billions of years.

This work is organized as follows: The overall method and the internal thermodynamical profiles from the evolutionary code MESA are described in Sect. 2, where we also summarize our 3D dynamo models that were performed with the MagIC code. In Sect. 3, we show the main results of our parameter exploration, interpret simulations representative of different evolutionary stages, and compare them to results from other works. We finally conclude in Sect. 4.

2 Method

It is currently not possible to simulate the realistic time evolution of the radial dependence of the thermodynamic quantities of a 3D planetary dynamo. The main reason is that the timescale associated with the convection inside gas giants tends to be on the order of years or decades. In contrast, the secular planetary cooling and contraction are appreciable at timescales on the order of gigayears. Due to this timescale separation of at least six orders of magnitude, a set of fixed backgrounds can be considered, and the dynamo models can be evolved for each of them. In other words, we aim at having different snapshots, that is, 3D dynamo solutions, each with fixed radial thermodynamical profiles that correspond to a given age of the 1D long-term evolutionary models. To schematically summarize, we employed a method that followed these steps:

1. Evolve a standard evolutionary 1D model of a contracting, nonirradiated gas giant over 10 Gyr.

2. For a given age of the evolutionary model, implement the radially dependent thermodynamical profiles as the background state of an anelastic spherical shell MHD model, with a given choice of the nondimensional numbers.

3. Evolve the 3D MHD equations in a spherical shell domain under the anelastic approximation and reach a self-sustained dynamo solution. Then let it evolve long enough (typically up to some ten thousand years of physical timescales) to average out the typical fluctuations and ensure statistical significance.

4. Repeat the process from step 2, where some of the nondimensional dynamo numbers rescale (compared to the first simulation) according to the relative variation of the involved thermodynamic quantities. Note that all the nondimensional numbers can be rescaled in this way, so that further assumptions about the rotation rate and viscosity, for instance, are required, as discussed below.

The process can then be repeated for another planetary model or a different choice of the reference values of the nondimensional numbers. In this approach, the trend of the dynamo numbers in the simulation sequence then includes the cooling information.

2.1 Internal structure

2.1.1 Long-term evolution

To model the evolutionary change of radially dependent thermodynamic quantities of gas giants, we used the public code MESA¹ (Paxton et al. 2011, 2013, 2015, 2018, 2019). This is a one-dimensional code that solves the time-dependent stellar structure equations and is capable of evolving low mass bodies, including brown dwarfs and gas giants (see Paxton et al. 2013). The equations we solved are the conservation of mass, hydrostatic equilibrium, energy conservation, and the energy transport equation, respectively, $$\frac{dm}{dr}=4\pi r^2\rho ,$$(1) $$\frac{dP}{dm}=-\frac{Gm}{4\pi r^4},$$(2) $$\frac{dL}{dm}=-T\frac{ds}{dt},$$(3) $$\frac{dT}{dm}=-\frac{GmT}{4\pi r^{4}P}\mathrm{\nabla },$$(4)

where m is the mass enclosed within a radius r, ρ is the density, P is the pressure, G is the gravitational constant, s is the specific entropy, T is the temperature, L is the internal luminosity, and ∇ ≡ d ln T/d ln P is the logarithmic temperature gradient, which was set to the smallest between the adiabatic gradient and the radiative gradient. In the energy Equation (3), the only source term we considered is the gravitational contraction. We neglected additional sources such as stellar irradiation (Guillot et al. 1996) or internal heat deposition (Komacek & Youdin 2017; Thorngren & Fortney 2018) from tidal (Bodenheimer et al. 2001) or Ohmic dissipation (e.g., Batygin & Stevenson 2010; Perna et al. 2010), or chemical processes such as hydrogen dissociation and recombination (Tan & Komacek 2019). These additional terms are fundamental for hot Jupiters (see Fortney et al. 2021 for a review), but negligible for cold, weakly irradiated planets. The set of equations is closed using the MESA equation of state (Paxton et al. 2019), which for the gas giant ranges of interest, is substantially the interpolation of the Saumon-Chabrier-van Horn equation of state for H-He mixtures (Saumon et al. 1995).

These 1D evolutionary models did not incorporate either diluted cores or hydrogen-helium demixing layers, that is, possible stratified layers in the convection interior, of the type mentioned in Sect. 1. For all our models, we assumed an interior inert rocky core of 10 M_⊕ with a homogeneous density ρ_c = 10 g cm⁻³, and a fixed solar composition for the envelope. To illustrate the evolutionary changes, we show in Fig. 1 the different profiles for two different planets with masses 1 and 4 M_J at ages 0.5, 1, and 10 Gyr. For comparison, we also show the profiles from the widely employed results of French et al. (2012) for the interior of Jupiter. As reported by (Paxton et al. 2013), during the evolution of the planet, the radius shrinks slowly, and after a few million years of evolution, this is independent of the chosen initial planetary radius. The planet slowly shrinks, and at early ages, the total radius of the 1 M_J model is therefore larger than that of the current Jupiter. The internal structure is always characterized by a very thin radiative layer with a thick, fully convective isentropic shell that encloses the inert core. The higher planetary mass (4 M_J) shows a larger radius, higher gravity and temperature, and lower thermal expansion coefficients α. However, the trends of the profile with age are similar to the 1 M_J case. Moreover, all models exhibit a nontrivial oscillating behavior of the Grüneisen parameter Γ (bottom panel; see below for the definition).

2.1.2 Background-state implementation

To run an anelastic MHD model, a series of thermodynamical quantities is required to be expressed as a function of radius. We implemented the MESA ρ(r), T(r), g(r), thermal expansion coefficient α(r), and the Grüneisen parameter Γ(r), all shown in Fig. 1, into a 3D model. We obtained α and Γ in terms of the readily available MESA thermodynamic profiles $$\mathrm{\Gamma }={\left(\frac{\mathrm{\partial }\ln T}{\mathrm{\partial }\ln \rho }\right)}_s-1\equiv {\mathrm{\Gamma }}_3-1,\alpha =-\frac{1}{\rho }{\left(\frac{\mathrm{\partial }\rho }{\mathrm{\partial }T}\right)}_P=-\frac{{\chi }_T}{T{\chi }_{\rho }},$$(5)

where $${\chi }_{\rho }\equiv {\left(\frac{\mathrm{\partial }\ln P}{\mathrm{\partial }\ln \rho }\right)}_T,{\chi }_T\equiv {\left(\frac{\mathrm{\partial }\ln P}{\mathrm{\partial }\ln T}\right)}_{\rho }.$$

Transport coefficients, which we discuss and prescribe below, are the other possible thermodynamical quantities that can show a radial dependence.

To define the radial domain for the 3D models, we cut the MESA profiles both at the inner and outer radial domain. The hydrostatic background is mostly isentropic due to the efficient convection. Close to the solid inert core, the profiles show a decrease in entropy, leading to a shallow stratified layer that is due to the boundary conditions. Since a more realistic modeling of the core is beyond the purpose of this study, we considered as the inner boundary of the convective shell for each model the region in which the profile is isentropic and cut out the innermost ≈2–3% (in radius) of the shell case by case. We do not expect these inner cuts to lead to significant differences because the results of Moore et al. (2022) showed that replacing a diluted core with a solid compact core of the same radius has little effect on the strength and morphology of the dynamo.

The density profiles usually spanned more than six orders of magnitude (the typical outermost layer in MESA is at a fraction of a bar), with the largest drop in the 1% outermost part of the planet, where no dynamo is expected. Spherical shell dynamo models cannot handle too large density contrasts. Therefore, as is common in other anelastic dynamo models, we cut the external layers to reduce the density ratio. We ensured that we always kept the pressure at which hydrogen metallization starts (≈1 Mbar) inside the domain. The maximum ratio we considered is ρ_i/ρ_o ≈ 100, and most of our models took values of ρ_i/ρ_o ≈ 20, which corresponds to approximate values for the number of density scale heights, N_ρ := ln (ρ_i/ρ_o), of 4.6 and 3.0, respectively. The thick profiles shown in Fig. 1 have N_ρ ≈ 3, and the thin endings represent an extension up to N_ρ ≈ 4.6. Throughout this work, we use the subindex notation o(i) for the value at the outer (inner) shell. These cuts define r_o, ρ_o, T_o, and P_o as well as their inner values, which are listed in Table 1. The thin radiative outer layer usually ends at about ≈1 bar, well above the external cut. Finally, all profiles except for Γ(r) were normalized to make the outer values equal to unity. This is a common practice for many 3D hydrodynamic codes as the fundamental units are not the physical units. MagIC, for example, works with units of the shell thickness as the length scale and with viscous timescales for units of time (see Sect. 2.2 for the exact details).

After the profiles were cut and normalized to the inner or outer boundary (as required by the code implementation; see below), we used high-degree polynomials to generally fit any shape the profiles can take. For ρ(r), T(r), and g(r), we employed a 20-degree Taylor series, which was enough to smoothly fit all cases tested here. α(r) and Γ(r) have some peaks and valleys that complicated the fitting procedure, however. With polynomials with degrees below ≈100–120, we found that the wiggles were poorly fit. We adopted a 150-degree polynomial for both, $$\begin{array}{cc}(\rho (r),T(r),g(r))& =\sum _{n=0}^{20}({\rho }_n,T_n,g_n)r^n,\\ (\alpha (r),\mathrm{\Gamma }(r))& =\sum _{n=0}^{150}({\alpha }_n,{\mathrm{\Gamma }}_n)r^n,\end{array}$$(6)

where r ranges from r_i to r_o. A change in the reference point of the expansion did not quantitatively improve the fits (we also tried with powers of (r−r_o),(r−r_i) and (r−r_o/2)). Therefore, for simplicity, we opted for powers of r, that is, a MacLaurin series. We also verified that the radial resolution we employed in the MHD code, that is, the number of Chebyshev polynomials (see Sect. 2.2), was more than enough to satisfactorily reproduce even the most complex radial profiles given by the implemented MacLaurin expansion.

Unlike the fit model used by Jones (2014), the background profiles taken from MESA are almost but not exactly isentropic. To quantify this deviation, we considered the quantity |ds/dr| ⋅r/s, which usually takes values of 10⁻⁴−10⁻⁵ with maxima near the outer cut regions of 10⁻³. This might lead to a slight energy imbalance resulting from the radial background profile itself. To ensure that this did not influence the overall dynamics, we analyzed the energy balances, that is, we compared the buoyancy power with viscous and Ohmic dissipation (see Sect. 2.2.3 for more details).

2.1.3 Transport coefficients

The profiles of transport coefficients did not come directly from MESA. Although some important ingredients have evolved, such as the particle density, realistic profiles for diffusivities require proper ab initio calculations. In particular, the electrical and thermal conductivities have to take into account the degenerate state of the electron population. The pressure ionization (rather than thermal) is non-negligible in the dense but relatively cold (compared to stars) convective interior. Moreover, the dynamo region arguably corresponds to the transition to the metallic phase for hydrogen. In this sense, French et al. (2012) calculated the electric conductivity σ for a set of (T, ρ) pairs along the modeled Jupiter adiabat. Further models about transport coefficients in pure H-He mixtures have been developed (usually neglecting the contribution of thermally ionized alkali metals, which becomes relevant where hydrogen is molecular; Kumar et al. 2021), with relative differences in the values of σ by factors of a few (see Bonitz et al. 2024 for a recent review). In any case, the trend is that there is a continuous and steep increase in the conductivity for increasing pressure up to around 1 Mbar, after which the dependence on both temperature and pressure (or density) is much milder.

To capture these fundamental properties, we adopted the electrical conductivity profile first defined in Gómez-Pérez et al. (2010), which consists of an approximately constant conductivity in the innermost hydrogen metallic region, with a polynomial plus exponential decay toward the outer molecular region², $$\frac{1}{\tilde{\lambda }(r)}:=\tilde{\sigma }(r)=\{\begin{array}{ll}1+({\sigma }_m-1){\left(\frac{r-r_i}{r_m-r_i}\right)}^{-a}& r<r_m,\\ {\sigma }_m{e}^{-a\left(\frac{r-r_m}{T_m-r_i}\right)\frac{{\sigma }_m-1}{{\sigma }_m}}& r\ge r_m,\end{array}$$(7)

where $\tilde{\sigma }$ and $\tilde{\lambda }$ are the normalized conductivity and magnetic diffusivity, respectively. The actual physical relation, λ = (μ₀ σ)⁻¹, with μ₀ being the vacuum magnetic permeability, is simplified when the quantities are normalized to their innermost values: $\tilde{\sigma }=\sigma /\sigma \left(r_i\right)$ and $\tilde{\lambda }=\lambda /\lambda \left(r_i\right)$. This expression ensures that both $\tilde{\lambda }$ and $d\tilde{\lambda }/dr$ are continuous at r_m, and that they qualitatively reproduce the main features. Moreover, it allowed us to compare results with several previous works that have employed this profile for gas giant convection and dynamo modeling (Duarte et al. 2013, 2018; Wicht et al. 2019b,a; Gastine & Wicht 2021). They used values of σ_m and a that ranged from approximately 0.9 to 0.01 and from 1 to 25. In our models, we fixed σ_m = 0.1 and a = 7, while r_m was chosen as the radius at which each MESA planetary profile reached 1 Mbar, that is, the approximate pressure above which hydrogen is thought to undergo metallization. The profiles of $\tilde{\sigma }$ are shown in Fig. 2 for the same representative models as in Fig. 1.

On the other hand, for simplicity, we kept the kinematic viscosity v and thermal diffusivity κ constant within the same model. In Sect. 3.6, we return to the impact and caveats of this choice and the related assumption about the Prandtl numbers (Sect. 2.2.4).

2.2 3D numerical dynamo model

The next step is fundamental: We performed 3D MHD spherical shell simulations using the public code MagIC³ with the anelastic approximation (Gastine & Wicht 2012). This is a pseudo-spectral code that uses the spherical harmonic decomposition in the angular directions, that is, θ and θ, and Chebyshev polynomials in the radial direction r. MagIC has been used for both stellar and planetary models, including for Jupiter and Saturn dynamos, for instance, Duarte et al. (2018); Wicht et al. (2019b); Gastine & Wicht (2021); Yadav et al. (2022), and it was tested in convection and dynamo benchmarks Christensen et al. (2001); Jones et al. (2011). We employed the anelastic approximation, which is typically used for modeling the density-stratified low-Mach number convection flows in gas giants and stars Braginsky & Roberts (1995); Lantz & Fan (1999).

The shell is filled with a finitely conducting fluid rotating along the vertical axis $\hat{z}$ with a constant angular velocity Ω, and the background profiles, which we set up as explained in section 2.1.2. The geometry is set by the aspect ratio η := r_i/r_o. We work in dimensionless units using the shell thickness d := (r_o−r_i) as the length unit and the viscous diffusion timescale d²/v as the time unit. The magnetic fields are in units of (ρ_o μ_o λ_i Ω)^1/2, where λ_i = 1/(μ₀ σ(r_i)) is the magnetic diffusivity at the inner boundary. The convection is set by a fixed entropy gradient Δs, and this difference serves as the nondimensional units for s.

The equations we solved were the mass continuity equation, the momentum equation, the entropy equation, the induction equation, and the solenoidal constraint for the magnetic field, $$\mathrm{\nabla }\cdot (\tilde{\rho }\mathbf{u})=0,$$(8) $$\begin{array}{cc}\frac{\mathrm{\partial }\mathbf{u}}{\mathrm{\partial }t}+\mathbf{u}\cdot \mathrm{\nabla }\mathbf{u}=& -\mathrm{\nabla }\left(\frac{p^{\mathrm{\prime }}}{\tilde{\rho }}\right)-\frac{2}{E}{\mathbf{e}}_{\mathbf{z}}\times \mathbf{u}-\frac{Ra}{Pr}\tilde{g}\tilde{\alpha }\tilde{T}{s}^{\mathrm{\prime }}{\mathbf{e}}_{\mathbf{r}}\\ & +\frac{1}{PmE\tilde{\rho }}(\mathrm{\nabla }\times \mathbf{B})\times \mathbf{B}+\frac{1}{\tilde{\rho }}\mathrm{\nabla }\cdot \mathbf{S},\end{array}$$(9) $$\tilde{\rho }\tilde{T}(\frac{\mathrm{\partial }{s}^{\mathrm{\prime }}}{\mathrm{\partial }t}+\mathbf{u}\cdot \mathrm{\nabla }{s}^{\mathrm{\prime }})=\frac{1}{Pr}\mathrm{\nabla }\cdot \left(\tilde{\rho }\tilde{T}\mathrm{\nabla }{s}^{\mathrm{\prime }}\right)+\frac{\mathrm{Pr}Di}{Ra}(Q_v+Q_{\lambda }),$$(10) $$\frac{\mathrm{\partial }\mathbf{B}}{\mathrm{\partial }t}=\mathrm{\nabla }\times (\mathbf{u}\times \mathbf{B})-\frac{1}{P{m}_i}\mathrm{\nabla }\times (\tilde{\lambda }\mathrm{\nabla }\times \mathbf{B}),$$(11) $$\mathrm{\nabla }\cdot \mathbf{B}=0,$$(12)

where the traceless rate-of-strain tension S_ij an the viscous and Ohmic heating terms, Q_ν and Q_λ are defined by $$\begin{array}{cc}{S}_{ij}& \equiv 2\tilde{\rho }(e_{ij}-\frac{1}{3}{\delta }_{ij}\mathbf{\nabla }\cdot \mathbf{u}),e_{ij}\equiv \frac{1}{2}(\frac{\mathrm{\partial }{u}_i}{\mathrm{\partial }{x}_j}+\frac{\mathrm{\partial }{u}_j}{\mathrm{\partial }{x}_i}),\\ Q_{\nu }& \equiv 2\tilde{\rho }(e_{ij}{e}_{ij}-\frac{1}{3}(\mathrm{\nabla }\cdot \mathbf{u}{)}^2),Q_{\lambda }\equiv \frac{\tilde{\lambda }}{P{m}^{2}E}(\mathrm{\nabla }\times \mathbf{B}{)}^2.\end{array}$$

The dimensionless Ekman, Rayleigh, Prandtl, and magnetic Prandtl numbers are respectively defined as $$E\equiv \frac{v}{\mathrm{\Omega }{d}^2},Ra\equiv \frac{{\alpha }_o{g}_o{T}_o{d}^3\mathrm{\Delta }s}{c_{p}v\kappa },Pr\equiv \frac{v}{\kappa },Pm\equiv \frac{v}{{\lambda }_i}.$$

All the quantities marked with a tilde (Sects. 2.1.2 and 2.1.3) are static in time, radially dependent, and normalized to their outer values, except for 1/λ, which because of its fast decay in the outer radial regions, was normalized to its innermost value for practical reasons.

2.2.1 Boundary conditions

We assumed stress-free and impenetrable boundary conditions for the velocity at the inner and outer radii, r = r_i, r_o, $$u_r=\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\frac{u_{\theta }}{r}\right)=\frac{\mathrm{\partial }}{\mathrm{\partial }r}\left(\frac{u_{\varphi }}{r}\right)=0.$$

When stress-free boundary conditions are used at both boundaries, MagIC gives the option to ensure angular momentum conservation, which we have used. We employed constant entropy at both boundary conditions, $$s^{\mathrm{\prime }}(r=r_0)=0,s^{\mathrm{\prime }}(r=r_i)=1.$$

The material outside the outer radius was electrically insulating, that is, the magnetic field matched a potential field. At the inner boundaries, we imposed a perfectly conducting core.

2.2.2 Numerical technique

MagIC solves the set of Equations (8)–(12) with the above boundary conditions by expanding the mass flux and the magnetic fields into poloidal and toroidal potentials, $$\begin{array}{cc}\tilde{\rho }\mathbf{u}& =\mathbf{\nabla }\times (\mathbf{\nabla }\times W{\mathbf{e}}_{\mathbf{r}})+\mathbf{\nabla }\times Z{\mathbf{e}}_{\mathbf{r}},\\ \mathbf{B}& =\mathbf{\nabla }\times (\mathbf{\nabla }\times g{\mathbf{e}}_{\mathbf{r}})+\mathbf{\nabla }\times h{\mathbf{e}}_{\mathbf{r}}.\end{array}$$

The quantities W, Z, g, h, s^′ and p^′ are expanded up to l_max in spherical harmonic degree and N_C in Chebyshev polynomials, where we always set N_C = N_r−1. The equations are time-stepped by advancing nonlinear and Coriolis terms using an explicit second-order Adams–Bashforth scheme, and the remaining terms were time-advanced using the implicit Crank-Nicolson algorithm (for more details, see Glatzmaier 1984; Christensen & Wicht 2007).

2.2.3 Diagnostic parameters

To characterize the numerical dynamo solutions, we made use of several diagnostic quantities. We typically took averages in time or in space, either over the whole volume V or over spherical surfaces, to show the radial dependence, $$\begin{array}{cc}\parallel a\parallel (r,t)& =\int a(r,\theta ,\varphi ,t)sin\theta d\theta d\varphi ,\\ ⟨a⟩(t)& =\frac{1}{V}\int a(r,\theta ,\varphi ,t)dV,\overline{a}=\frac{1}{\mathrm{\Delta }t}{\int }_{t^{\mathrm{\prime }}}^{t^{\mathrm{\prime }}+\mathrm{\Delta }t}a(t)dt.\end{array}$$

We performed the time averages after a stationary state had been reached, and we typically monitored the dimensionless hydrodynamic and magnetic Reynolds numbers, the Rossby number, and the Elsasser number, $$\begin{array}{c}Re=\overline{\sqrt{⟨u^2⟩}},Rm=\overline{\frac{1}{V}{\int }_{r_i}^{r_o}\frac{\sqrt{\parallel {\mathbf{u}}^2\parallel }}{\tilde{\lambda }}{r}^{2}dr},\\ Ro=\frac{RmE}{Pm},\mathrm{\Lambda }=\overline{⟨\frac{{\mathbf{B}}^2}{\tilde{\rho }\tilde{\lambda }}⟩.}\end{array}$$

The total kinetic and magnetic energy, which are in units of ρ₀d⁵E²Ω², were $$E_{kin}=\frac{1}{2}\overline{⟨\tilde{\rho }{\mathbf{u}}^2⟩},E_{mag}=\frac{1}{2}\overline{\frac{1}{EPm}⟨{\mathbf{B}}^2⟩}.$$

We defined the dipole fraction, f_dip = E_mag,l=1/E_mag, as the ratio of the magnetic energy stored in dipolar components (axisymmetric and nonaxisymmetric), divided by the total magnetic energy⁴.

To study the energy dissipation, we used the buoyancy power, $$P_{\nu }(t)\equiv \frac{RaE}{Pr}⟨\tilde{\alpha }\tilde{T}\tilde{g}{s}^{\mathrm{\prime }}{u}_r⟩.$$(13)

We used the subindex v to emphasize that the quantity was calculated in viscous timescales. For comparison with other works, we also used the rotation timescale (see Sect. 3.7). For a well-resolved numerical run, the buoyancy power must be equal to the sum of viscous and Ohmic dissipation rates after a steady-state solution is reached. These are defined as $$D_{visc}(t)\equiv ⟨S^2⟩,D_{ohm}(t)\equiv \frac{1}{EP{m}^2}⟨\tilde{\lambda }(\mathbf{\nabla }\times \mathbf{B}{)}^2⟩.$$(14)

Another quantity to monitor is the fraction of energy dissipated by Joule heating alone, that is, the Ohmic fraction f_ohm = D_ohm/P_γ. After a statistically steady state has been reached, the input buoyant power must balance the viscous and Ohmic diffusion. To evaluate whether the numerical solution has good time and spatial invariance and also if the background state affects the energy balance strongly, we assessed the power imbalance by its proxy $f_P=|\overline{P_v}-\overline{D_{visc}}-\overline{D_{ohm}}|/\overline{P_v}$.

Finally, we studied the time-averaged kinetic and magnetic spectra, that is, the distribution of the energy over different multipoles of order l, which MagIC already implemented as a user-friendly output. We inspected the spectra for each model, in particular, to ensure that the resolution was high enough to resolve the maximum dissipation, that is, l(l+1) E(l), for the volume-integrated spectra and 2D spectra taken at relevant radii, for instance, r_i, r_o, r_m.

2.2.4 Parameter evolution and model descriptions

As explained in 2.1.2, when the MESA profiles were cut close to the desired N_ρ, we extracted ΔT, r_o, r_i, and r_m, from which we deduced η and χ_m := r_m/r_o. The corresponding physical values can be recovered by knowing the units in which each quantity is expressed and the values from the MESA profile, for example, using the thickness of the physical shell thickness d_phys = r_{o, phys}(1 − η). The real quantities of the planet, which reflect the evolutionary changes, were used to evolve the dynamo parameters. Since, as mentioned above, the real physical values E and Ra are computationally inaccessible, we can still use their dependence on the physical values that change during the long-term evolution. In particular, the shell thickness d = r_o−r_i and the temperature difference ΔT enter the definition of the Ekman, E(t) ≈ d(t)⁻², and Rayleigh numbers, Ra(t) ≈ d(t)³ Δ T(t). Therefore, we considered a series of ages, for which, after having found a suitable pair of E₀ and Ra₀ that produces convection and dynamo for the setup with d₀ and ΔT₀ corresponding to a given age, the values E^′ and Ra^′ of the remaining models in that series were set up by scaling with d(t) and Δ T(t), $$E^{\mathrm{\prime }}=E_0\frac{d_0^2}{d^{\mathrm{\prime }2}},R{a}^{\mathrm{\prime }}=R{a}_0\frac{d^{\mathrm{\prime }3}\mathrm{\Delta }{T}^{\mathrm{\prime }}}{d_0^3\mathrm{\Delta }{T}_0}.$$

We made two assumptions: (i) The diffusivities at a given radius remain constant in time, which implies that we considered the same values of Pr and Pm along a sequence and that the change in E only comes from the contraction of the planet; and (ii) planetary rotation is constant in time. The latter assumption is well justified when we consider the possible relevant torque acting on a gas giant. Batygin (2018) studied the evolution of rotation, considering the magnetic coupling between the planetary interior and the quasi-Keplerian motion of the disk in the planetary formation stages. This resulted in efficient braking of the planetary spin that ceased to evolve after ≈1 Myr, when it reached a terminal rotation rate (probably similar that of Jupiter), which can hardly change later. Our earliest model is at 100 Myr, for which the rotation can be safely considered constant. In this sense, cold giants are expected to be fast rotators (Ro < 0.12), and they might host planetary dynamos similar to those found in dipole-dominated numerical solutions with very low E (Davidson 2013; Yadav et al. 2016; Schwaiger et al. 2019). In this scenario, there is a quasi-geostrophic balance (Coriolis and pressure forces) at the largest scales, followed by an ageostrophic magneto-Archimedean-Coriolis balance (Coriolis, buoyancy, and Lorentz forces).

With these assumptions, we considered five sets of dynamo models: (i) A long series with a total of 12 evolutionary stages, ranging from 0.1 to 10 Gyr for a 1 M_J planet; (ii) different density ratios, that is, different cutoff radii, for the same 1 M_J model at 1 Gy; (iii) different planetary masses with N_ρ ≈ 3.0; (iv) several models with Pm and Pr different from 1; and (v) a 4 M_J mass series with N_ρ ≈ 4.6. When a different mass was chosen, the radial profiles changed (see Fig. 2), so that using Eq. (7) with the above-mentioned values of the free parameters a and ${\tilde{\sigma }}_m$, we had to adapt the density contrast N_ρ to include the drop in the conductivity in the outer layers of our shell, without at the same time, having too low values of σ. For this reason, the series of 4 M_J has a higher N_ρ (the exponential drop of σ would have been cut out with N_ρ ≈ 3). For the same reason, the 0.3 M_J model has a lower contrast, N_ρ ≈ 1.1. The alternative would have been to consider an equally arbitrary change of the free parameters in Eq. (7). In Table 1 we show the input values for the 3D simulations together with the parameters of the background profiles from the MESA 1D long-term evolution.

3 Results

3.1 Preliminary exploration of parameters

Our first goal was to determine the dependence of the dynamo solutions on mass and age for a given evolutionary sequence of background thermodynamic setups (Table 1). Therefore, the first requirement was to determine a range of parameters for which both convection and dynamo operated for all models. This practically meant that we needed to find a feasible range of Ra and E (which, within a sequence of models, have relative variations set by Δ T, T_o and d, Table 1), for which convection and dynamo action are present in the entire sequence, keeping in mind that the chosen values were orders of magnitude away from the realistic values, as mentioned above. To determine these feasible ranges, we needed to perform a preliminary parameter exploration. At the same time, we assessed the sensitivity of results on other parameters. We summarize in this subsection the four main steps we took for this overall assessment.

First, to locate the region with viable dynamo solutions, we performed low- and medium-resolution runs for one specific model, the model with 1 M_J 10 Gyr. This is the oldest and coldest model of the 1 M_J sequence, that is, with the lowest value of ΔT, which also implies the lowest value Ra of the series, that is, it is least favorable to convection. We spanned the ranges 10⁻⁵ < E < 10⁻³, 10⁶ < Ra < 10¹⁰ with Pm = Pr = 1 and a relatively low resolution, (N_r, N_θ, N_θ) = (193, 192, 384). For E ≤ 10⁻⁴, we obtained convection for Ra ≳ 10⁷, and, additionally, magnetic field growth for Ra ≳ 5 ⋅ 10⁷.

Second, for the 1 M_J 10 Gyr case with E = 10⁻⁵, Ra = 5 ⋅ 10⁸ (the reference model in the rest of this subsection), we explored the Prandtl numbers in a relatively easily accessible range 0.25 < Pm, Pr < 4. This was done to discuss the impact of our assumptions of constant-in-time diffusivities Sect. 2.1.3). The results of this exploration are shown in Sect. 3.6 for high-resolution models and different evolutionary ages.

Third, for the same reference model (E = 10⁻⁵, Ra = 5 ⋅ 10⁸, Pm = Pr = 1), we explored the sensitivity to the parameters a and σ_m, which define the slope of $\tilde{\lambda }(r)$ in the outermost layers. For a wide range of values (i.e., 0.07 < σ_m < 0.9, 5 < a < 15), we indeed recovered the results of Duarte et al. (2013), that is, the strong equatorial jet remains confined to the weaker conducting outer region and does not interfere with the deeper dynamo action. As previously mentioned, we finally opted for rather low values for both σ_m of 0.1 and a of 7. These values are similar to those used by Gastine & Wicht (2021), whose model approximately reproduces the French et al. (2012) profiles. For numerical stability reasons, we did not use higher values of a, that is, steeper exponential drops, because the hydrogen metallic region of some of our models is already quite deep (χ_m > 0.85), which leads to too high values of $\tilde{\lambda }(r)$ near the outer surface. Similarly, we studied the effect of N_ρ on several diagnostic quantities (see Appendix A for more details). Given the nonnegligible effects of choosing different values of N_ρ, we compare below models with the same N_ρ, to avoid additional biases.

Fourth, we tested different boundary conditions again for the same reference model. Integrated quantities such as Rm or E_kin and convection patterns did not change appreciably when rigid boundary conditions were applied at the inner core. The radial distributions showed a drop in velocity in a very thin region (<1%) of the radius. For the magnetic field, we tested for a perfect conductor and insulator for the inner core and an insulating, perfect conductor and pseudo-vacuum at the outer radii. We found no relevant differences in the internal dynamo.

Considering the explored values of Ra and E for these low-resolution test, we set Ra = 1.3 ⋅ 10⁻⁵ and E ≈ 5 ⋅ 10⁸ for the 1 M_J 10 Gyr model with a N_ρ ≈ 3. This choice set the rescaled values for the remaining sequence in a range that allowed dynamo action. Using the definition of E, Ra, we scaled their values in each run with the corresponding values of T, ΔT, d, as described in Sect. 2.2.4. Most models use a resolution of (N_r, N_θ, N_ϕ) = (289, 256, 512). As an exception, the 4 M_J N_ρ ≈ 4.6 runs use a grid of (385, 320, 640).

Finally, we employed a strategy to save computational time that was inspired by other spherical shell dynamo works (Christensen & Aubert 2006). The idea is that the final solution does not depend on whether the initial conditions are taken from another saturated dynamo model or if they are the usual u = 0 with a weak perturbation of both B and T^′/s^′. This suits our case in particular because the relative changes in the dynamo parameters from one setup to the next are small, and the solution of the new setup is reached much faster than starting from a u = 0 state. In Appendix B we show the results of some specific numerical experiments that support this strategy.

We now proceed to the main results. We refer to Appendix D for a table with detailed values of the time-averaged output nondimensional numbers and the other quantities we used as diagnostics.

3.2 Dynamo solutions: General behavior

In Fig. 3, we show snapshots (maps and slices of some velocity and magnetic field components) of the 1 M_J 1 Gyr saturated dynamo solution, which is a representative case. The other models we obtained are qualitatively similar among themselves in terms of the morphology of the velocity and magnetic fields. The largest differences are the relative average strengths of u and B and the magnetic field dipolarity, which we discuss below. All the models have a strong equatorial flow that reaches deep, down to the dynamo region. The velocity and magnetic fields show a westward drift in the inner parts. The magnetic field is mostly constrained under r < r_m, where convection is also stronger, as seen in u_r. As expected from rotation-dominated convection, we found columnar structures in the direction of the rotation axis (Zhang & Busse 1987; Ardes et al. 1997; Simitev & Busse 2003), as shown in the meridional slices. Generally speaking, our numerical solutions are similar to those reported by other works for Jovian-like dynamos with a nonconstant electrical conductivity (Jones 2014; Duarte et al. 2018)⁵.

The dynamo solutions we found are generally less dipole-dominated than their incompressible (Boussinesq) counterparts with similar dynamo parameters (i.e., Ra, E, Pr, and Pm). However, we note that we did not explore values for Pr lower than 0.1 with Pm > 1, where dipole-dominated solutions have been found (Jones 2014; Tsang & Jones 2020). We restricted ourselves to a less demanding parameter space with a wider liberty of parameter exploration, but with the caveat that we might obtain less dipole-dominated models. Similarly to Yadav et al. (2013), for all runs shown here, we also obtained Nu > 2 at the top and bottom surfaces, which ensures a fully developed convection. The Nusselt number Nu is the ratio of the total transported heat flux to the conducted heat flux.

Moreover, for gas giants, the definition of the dynamo surface is not absolute. As hydrogen gradually transitions outward from metallic to molecular, the electrical conductivity and electrical currents are quickly (but not abruptly) damped over a finite region. For our models, the most obvious choice for the dynamo surface is the radius at which the exponential decay for σ starts, that is, r_m. To obtain a more physically justified definition for the dynamo surface, we followed Tsang & Jones (2020) by computing the magnetic energy spectra at different radii, F_l(r). They defined the dynamo surface, r_dyn, as the radius within which the slope of the Lowes spectrum (i.e., a potential solution extrapolated back from the outermost layer to the interior) diverges from the slope of the simulated F_l(r). A similar analysis for some of our models is shown in Appendix C, where we find that this definition of the dynamo surface always gives values very close to r_m.

3.3 Evolutionary changes

We now focus on the variation in the solutions along the longest sequence of models, the 1 M_J planet with N_ρ ≈ 3 and Pm = Pr = 1. The shell-averaged spectral distribution and radial distribution of the magnetic (solid lines) and kinetic (dashed lines) energy are shown in Fig. 4. We focus on three representative ages: 0.4, 2.1, and 6.5 Gyr.

In the left panel, the kinetic spectra show a drop of about 1.5 orders in magnitude or more from the integral to viscous scales, while the magnetic spectra decrease by at least 3 orders. The sawtooth shape on the lowest multipole side is associated with the external jet that we see in all of our runs. This behavior is not seen for the m spectrum, which is dominated by the zonal flows, m = 0. For higher multipoles, the spectrum plateaus before it reaches the viscous scale and drops off. A comparison of the three different models shows that the overall shape of the kinetic spectra does not significantly change, other than a constant decrease in time throughout all harmonic degrees. The magnetic spectra show a similar diffusive scale, which is approximately located at the same l as the viscous diffusive scale (compatible with Pr = Pm = 1), although the knee is less pronounced. With the usual measure of dipolarity, it ranges from $0.3>f_{dip,l<12}^{axi,surf}=E_{mag}(r_0{)}_{l=1,m=0}/E_{mag}(r_0{)}_{l\le 12}>0.8$. Depending on the work, this could be considered multipolar or dipolar (Christensen & Aubert 2006; Yadav et al. 2013; Zaire et al. 2022). The magnetic spectra show a clear evolution with age: The relative weights of high multipoles tend to decrease, while the strength of the large scales slightly increases. This inversion leads to an increase in the total dipolarity (see the discussion below).

The radial energy distributions, shown in the right panel of Fig. 4, shows that E_kin(r) increases almost monotonically outward, but the steepest changes occur in the outermost layers, r ≳ r_m, due to the appearance of the equatorial zonal wind, where the density is lower and the magnetic drag is weaker than in the interior. Similarly to the spectra, the kinetic radial distribution does not show a clear variation with age. The radial profile E_mag(r) and its change with age shown also in Fig. 4 is instead more complex. The radial profile peaks at radii slightly smaller than r_m, after which it significantly drops, following the σ(r) profiles (Fig. 2). A comparison of different evolutionary ages shows that the innermost region of the radial distribution does not show a clear trend, with a slight increase in the deepest regions for late-age models. On the other hand, the layers ≈10–20% below r_m show a steady decrease with age (see Sect. 3.5).

The overall changes between the runs are gradual, and we found that a few models can already predict the general behavior. Several time and volume-averaged diagnostic quantities are shown as a function of evolutionary models in Fig. 5. As expected during the planetary cool-down of a gas giant, Rm decays approximately like a power law in time. Ro and P_ν also behave similarly, and they are therefore not shown to avoid repetition. All of these quantities are dependent on u_rms or at least on one of its components. This reflects the fact that the mean velocity is dictated by the buoyancy input parameter Ra, which is proportional to the temperature difference in the convective shell. Fig. 5 also shows that the dipolarity f_dip does not seem to change significantly in time, except for a mild increase between 1.5 and 3.8 Gyr. This appears to indicate a transition between a multipolar or weakly dipolar-dominated regime to a strong dipolar-dominated regime. This increasing trend is obscure, with $f_{dip,l<12}^{axi,surf}$ (for the 1 M_J series, they take values from 0.3 to 0.8). As an alternative, we obtained the average f_dip over the volume 5% near the dynamo surface, f_dip,dyn. Fig. 5 shows that f_dip,dyn grows gradually, with a similar jump seen in f_dip.

This transition from multipolar to dipolar dynamo was also observed by Zaire et al. (2022) for dynamos in stratified stellar interiors, which they modeled with shallower shells (r_i/r_o = 0.6) and different N_ρ and Ra. They obtained a threshold F_I/F_L (i.e., the relative importance of the inertial over Lorentz forces) below which multipolar dynamos collapse into dipolar dynamos. They also reported that E_kin/E_mag can equally well capture this magnetic morphology transition and found this transition at about E_kin/E_mag = 0.7. In Fig. 6 we show both $f_{dip,l<12}^{axi,surf}$ and f_dip,dyn as a function of E_kin/E_mag for the 1 M_J, Pr = Pm = 1 series. We also report two distinctive populated areas, low dipolarity with high E_kin/E_mag and high dipolarity with lower E_kin/E_mag. This abrupt change in the magnetic field morphology seems to be better reflected with f_dip,dyn. A good definition for a dipole-dominated dynamo could be f_dip,dyn > 0.1, which might be compatible with the definitions of Yadav et al. (2013) or Zaire et al. (2022) of $f_{dip,l<12}^{axi,surf}>0.3$ and 0.5, respectively. These dipole-related quantities are usually the most fluctuating integrated quantities in saturated dynamo solutions because they are susceptible to the specific magnetic field configuration.

In any case, our strongest evolutionary trend shows a transition from a multipolar to a dipolar regime in the middle of our series. Within the more multipolar part of the series, Λ and f_ohm show a minor decay in time, and E_mag/E_kin also seem to decrease, but more subtly. In the dipolar regime, these trends reverse: Λ and f_ohm plateau and show a slight increase; the ratio E_mag/E_kin grows noticeably. Following the magnetic fields of Jupiter and Saturn (Connerney et al. 2022; Cao et al. 2023), gas-giant dynamos are expected to live in a parameter space region with dipole-dominated solutions, and thus, the expected evolution would be of the latter part of our series. Moreover, the Jovian value for Pm is expected to be ≈10⁻⁶, meaning that D_ohm ≫ D_visc, and thus, f_ohm is expected to be close to one. This would mean that our predicted f_ohm trend is not physically noticeable. An increase in E_mag/E_kin and Λ, the integrated nondimensional magnetic energy, might contradict the aforementioned scaling laws, but we show their compatibility in Sect. 3.5.

The values of E_mag/E_kin that we show are below equipartition (i.e., 0.1 < E_mag/E_kin < 0.6 for the long 1 M_J sequence) because we analyzed a volume including the nonconducting outer layer. When we restrict the energy integration within the metallic region r < r_m, then 0.25 < E_mag/E_kin|_{r < rm} < 1.4. This tendency can be sensed from the radial distributions in Figs. 4 and 7.

3.4 Dependence on planetary mass

We have obtained saturated models for five different masses at 5 Gyr. We directly compared only four of them because the conductivity of the 0.3 M_J model was not numerically feasible with N_ρ ≈ 3 (see above). The model of 4 M_J is slightly under-resolved, possibly because of its higher Ra and a very little conductivity drop, which does not help to stabilize the stress-free boundary conditions. We did not use higher-mass 1D models (specifically, the 8 and 12 M_J) because of the resolution constraints required by the E, Ra combination.

The main dimensionless diagnostics for these runs were already shown in Fig. 5. The overall trends are dictated by the increase of Ra with a decrease in E that comes with the 1D profiles themselves. Therefore, as the mass increases, both f_dip and f_dip,dyn decrease, while f_ohm, Λ, and Rm increase. The energy ratio is approximately maintained.

In Fig. 7 we plot the energy spectra and radial distribution for these models. The key features are very similar to the feature described above (Fig. 4). A noticeable difference is that the magnetic spectra seem to become flatter for higher masses. The sawtooth shape in the kinetic spectra diminishes with mass because depth of the outer nonconductive layer decreases as the hydrogen metallization pressures are reached faster. To overcome this difference and to observe whether the evolutionary trends shown in Fig. 5 changed for other masses, we obtained saturated dynamos for the 4 M_J 0.5, 1, and 10 Gyr N_ρ ≈ 4.6 models. For these three, we obtained similar dynamos with larger equatorial jets, in other words, we recovered the sawtooth shape seen in the kinetic spectrum. The evolutionary trends match the multipolar side of the 1 M_J long series.

3.5 Evolution of the magnetic field strength at the dynamo surface

The scaling law provided by Reiners et al. (2009) gives the mean strength of the magnetic field at the dynamo surface, B_dyn, in terms of the mass M, luminosity L, and radius R of the substellar object, $$B_{dyn}={4.8}_{-2.8}^{+3.2}{\left(\frac{M}{M_{\odot }}\right)}^{1/6}{\left(\frac{L}{L_{\odot }}\right)}^{1/3}{\left(\frac{R_{\odot }}{R}\right)}^{7/6}\mathrm{k}\mathrm{G}.$$(15)

We evaluated B_dyn by using the L(t), R(t) output from our MESA simulations (pink lines in Fig. 8).

Another slightly different estimate comes from inserting the analytical expressions for L(t) and R(t) of Burrows & Liebert (1993); Burrows et al. (2001) in Eq. (15) as given for a substellar-mass solar-metallicity object, $$L\approx 4\cdot {10}^5{L}_{\odot }{\left(\frac{1Gy}{t}\right)}^{1.3}{\left(\frac{M}{0.05M_{\odot }}\right)}^{2.64},$$(16) $$R\approx 6.7\cdot {10}^4\mathrm{k}\mathrm{m}{\left(\frac{{10}^5\mathrm{c}\mathrm{m}{\mathrm{s}}^{-2}}{g}\right)}^{0.18}{\left(\frac{T_{eff}}{1000\mathrm{K}}\right)}^{0.11}.$$(17)

Using these estimates, Reiners & Christensen (2010) obtained a slow decay of the dynamo magnetic field of about one order of magnitude in about 10 Gyr. Fig. 1 in their paper shows that the approximate power-law relation is B_dyn ≈ t^−0.3 (marked with a gray line in Fig. 8).

We then compared these methods with ours. To do this, we evaluated the values of B_dyn for our models. As mentioned above, we made use of the fact that the effective dynamo surface is located at r_m. We decided to obtain the volume average of E_mag over a spherical shell from r_m to some not-too-deep layer, $$E_{mag,dyn}(q)=\frac{1}{r_m-r_m^{\mathrm{\prime }}}{\int }_{r_m^{\mathrm{\prime }}}^{r_m}{E}_{mag}(r)dr,$$(18)

where $r_m^{\mathrm{\prime }}(q)=({\chi }_m-q)r_o$, and E_mag(r) was obtained from the previously shown radial distributions. The shell thickness is therefore controlled by the parameter q, which in turn allowed us to evaluate the surface dynamo field as $B_{dyn}(q)=\sqrt{2{\mu }_0{E}_{mag,dyn}(q)/MV}$. In Fig. 8 we show the estimated B_dyn for q = 0.05, 0.1, 0.15, and 0.2 (colored points in Fig. 8, with the related statistical error). The best-fitting slope (dotted lines) decreases with thicker integrating regions, that is, larger q. Within the standard deviations, we recover the previously mentioned slope of ≈t^−0.3 for the most uncertain slope (thinnest averaging region). The others are slightly shallower than the trends obtained by Reiners & Christensen (2010) and Reiners et al. (2009).

3.6 Dependence on the Prandtl numbers

To assess the impact of the assumption of constant Pr and Pm, we obtained several saturated dynamo states with Pr, Pm ≠ 1 for some 1 M_J models. We performed runs with 0.5, 1, and 10 Gyr, which we deemed enough for assessing general properties of the trends. We investigated within the following range: 0.5 < Pm < 4 and 0.5 < Pr < 2. The results are shown in Fig. 9.

The increase in Pm can be understood as lowering λ while keeping v constant, or similarly, increasing v with constant λ. Both effects lead to an increase in the magnetic energy in relation to the available buoyant power. Therefore, Fig. 9 shows that Rm and Λ tend to increase with Pm, which is a more efficient dynamo mechanism (Elias-López et al. 2024). The same applies for E_mag/E_kin and f_ohm, as the decreasing λ increases the magnetic energy percentage and the Ohmic dissipation contribution. The surface and volumetric dipolarities both decrease with increasing Pm, which is compatible with what was found by Tsang & Jones (2020): a higher Pm means a less steep magnetic spectrum, or in other words, a more weakly dipole-dominated magnetic spectrum.

In contrast, Rm, Λ, and f_ohm decrease for increasing Pr. This can easily be understood by considering a higher Pr as increasing values of v in comparison to κ. A higher viscosity will lead to lower kinetic as well as magnetic energy, and therefore, to lower Rm and Λ. The decrease in f_ohm means that Ohmic dissipation becomes less important than viscous dissipation. The magnetic energy ratio, E_mag/E_kin, and the two dipolarities, f_dip and f_dip,dyn, increase with Pr because it leads to a more efficient dynamo mechanism.

With these trends in mind, we describe the effect of the constant Pr and Pm assumptions on the obtained trends in Figs. 5 and 8. A slight increase in Pr is expected to occur during the long-term evolution of the planets because while it cools down the ratio of the thermal and electrical conductivities (inversely proportional to Pr) decreases according to the Widemann-Franz law, which is valid in the metallic region (French et al. 2012). In contrast, the viscosity and conductivity themselves are not thought to vary appreciably with temperature (i.e., in time) (French et al. 2012; Bonitz et al. 2024). Therefore, by using the trend of Λ with Pr, we expect that for an evolutionary change in Pr, the trend B_dyn(t) would be slightly steeper, which might agree even better with the Reiners & Christensen (2010) scaling law trend. However, a firm conclusion about this based on a modified setup for diffusivities and Pr evolution is left for future work.

In general, the evolutionary trends shown in Fig. 5 are maintained regardless of Pr and Pm. Rm decreases similarly for all set runs. Finally, E_mag/E_kin, Λ, and f_ohm show a different behavior depending on their dipolarity. The two most strongly dipolar solutions (Pm = 1, Pr = 2 and Pm = 2, Pr = 2) consistently show the same behavior as noted above, that is, an increase in E_mag/E_kin and f_ohm in time, and a plateau or mild decrease in Λ. In contrast, the most multipolar set of Prandtl numbers (Pm = 4, Pr = 1) decreases in f_ohm and Λ and increases slightly less in E_mag/E_kin. The other sets of runs are consistent with a transition from multipolar to dipolar similar to Fig. 5.

3.7 Scaling laws

Yadav et al. (2013) used a large set of anelastic dynamo numerical solutions to derive scaling laws by relating several representative dimensionless diagnostic parameters. Their dynamos covered a large space of parameters: 0 ≤ N_ρ ≲ 5.5, 0.1 ≤ η ≤ 0.75, 0.3 ≤ Pr ≤ 10, 0.2 ≤ Pm ≤ 20, 10⁻ 6 ≤ E ≤ 10⁻ 3, and 2.5. 10⁵ ≤ Ra ≤ 2.5 ⋅ 10⁹. The scattering plots shown in Figs. 10 and 11 superimpose the results of Yadav et al. (2013) with the data representing our models.

To compare our results, we used the inverse rotation frequency Ω⁻¹ as the time unit and the magnetic field units of $\mathrm{\Omega }D\sqrt{{\mu }_o{\rho }_o}$. The buoyancy power P is $$P=\frac{Ra{E}^3}{Pr}\frac{⟨\tilde{\alpha }\tilde{T}\tilde{g}{s}^{\mathrm{\prime }}{u}_r⟩}{M},$$(19)

where M is the dimensionless mass of the shell and is listed in Table 1 for our models. In Fig. 10 we show Ro as a function of P/Pm^13/45. For the runs of Yadav et al. (2013), there is a clear separation between the models with constant σ, which lie close to the best-fitting power law (gray line), and the runs with a decaying σ profile near the surface, which lie slightly above but parallel to the trend. The reason is that for a decaying σ, the strong jets that appear in the external nonconductive layer tend to increase the total kinetic energy, and thus, Ro for the same amount of available P. Our definition of Ro differs by a factor of $1/\tilde{\lambda }$, which erases the outer jet contribution. Our runs therefore mostly lie above the scaling law itself. Our models represent different evolutionary stages of planetary dynamos that move through this dimensionless space. This progression is parallel to the power law Ro = 2.47 P^0.45 Pm^−0.13 and reaches from higher to lower values: For a single model, Ro decreases by about half an order of magnitude, while P/Pm^13/45 decreases by one order of magnitude (which corresponds to the 0.45 exponent). Physically, this evolution should be positioned orders of magnitude away, but if the scaling law holds, so should this trend.

We also defined the Lorentz number, Lo, as the nondimensional magnetic field strength in Ω⁻¹ time units per unit of mass. In this case, B is in units of $\mathrm{\Omega }d\sqrt{{\rho }_o{\mu }_o}$, and we relate the definition of Lo with Λ, E and Pm, $$Lo=\frac{B_{rms}}{\mathrm{\Omega }d\sqrt{{\rho }_o{\mu }_o}}\frac{1}{\sqrt{M}}=\sqrt{\frac{2E_{mag}{E}^2}{M}}.$$(20)

The other magnetically related scaling law involves the characteristic timescale of magnetic energy dissipation τ_mag, which is defined as the nondimensional magnetic energy divided by the joule heat dissipation (all in units of Ω⁻¹), $${\tau }_{mag}=\frac{E_{mag}{E}^2}{P{f}_{ohm}M}.$$(21)

In Fig. 11, we show the scaling laws for Lo and τ_mag with the same legend as in Fig. 10. Yadav et al. (2013) reported that dipolar- and multipolar-dominated solutions take similar but parallel trends. We overplot our runs with the joined dipolar and multipolar branches for the two scaling laws (they show them separately). As suspected from Sect. 3.3, our series evolves from the multipolar to the dipolar branch in both diagrams, but this is more clearly visible in the τ_mag plot. This is also an argument in favor for this possible multipolar to dipolar transition.

The power-law relations shown above are purely fits obtained from Yadav et al. (2013). The velocity scaling of Ro ∝ P^α with α somewhat larger than 0.4 was theoretically justified by force balances by some authors (Aubert et al. 2001; Davidson 2013; Starchenko & Jones 2002), but none derived a Pm dependence. Similarly, τ_mag ∝ Ro^α was also discussed, where α ≲−1 (Christensen & Tilgner 2004; Stelzer & Jackson 2013) or α ≈ − 0.75 (Davidson 2013). Finally, if the magnetic field is only a function of power, dimensional arguments dictate that it must depend on the cubic root of the power, that is, Lo ∝ P^1/3 (Kunnen et al. 2010; Christensen et al. 2009; Davidson 2013).

4 Conclusions

We used radial profiles taken from gas giant evolutionary models to obtain sequences of 3D MHD spherical shell dynamo models. From the public code MESA, we obtained the radial hydrostatic profiles for planets with different masses at different stages of evolution: 0.3 M_J ≤ M_P ≤ 4 M_J and 0.2 Gyr ≤ t ≤ 10 Gyr, respectively. From the evolutionary tracks, we derived the trends for the dynamo parameters. Using the radial profiles as the background state, we solved the resistive MHD equations under the anelastic approximation with the pseudospectral spherical shell code MagIC. We obtained saturated dynamo solutions and interpreted them as different snapshots of planetary dynamos during their long-term evolution.

For our longest set of runs that represents different evolutionary times of a 1 M_J planet, we find a transition from a multipolar- to a dipolar-dominated dynamo regime. Within the dipolar or multipolar branch, a few snapshots are enough to generally assess the behavior that cannot be straightforwardly derived from scaling laws. As the planet evolves and cools, we obtained a steady decrease for Rm, P, and Ro as well as an increase in the volumetric and surface dipolarities. For multipolar dynamo solutions Λ, f_ohm and E_mag/E_kin decrease with time, whereas they increase for dipolar dynamos. These quantities are a proxy for the magnetic field energy, the dynamics of the power dissipation, and the energy ratio. These trends hold for different Prandtl numbers and for the 4 M_J models. Future studies are needed to confirm that these trends are also observed at lower Ekman numbers because the physical nondimensional parameter space is currently computationally inaccessible.

The decay of the magnetic field strength at the dynamo surface (Fig. 8) is roughly compatible with existing scaling-law estimates. Our models are based on a sequence of realistic backgrounds, and thus, they can be representative of the real long-term dynamo evolution. The trends in mass and age are also expected to hold for the realistic but computationally unfeasible range of nondimensional numbers (which is the usual intrinsic caveat of any dynamo study). We compared our results with the anelastic scaling laws of Yadav et al. (2013) between nonlinear combinations of dimensionless diagnostics, and we showed that the long-term evolution of a cold Jupiter dynamo evolves in this parameter space. Consistently with the diagnostics on the morphology described above, some age sequences transition from the multipolar-dominated to the dipolar-dominated family of solutions, which are separated in this parameter space, as shown by Yadav et al. (2013).

Additionally, in the sequences, we considered fixed values of Pr and Pm, and the values of Ra and E changed only due to the evolving values of ΔT, T_o, and the shell thickness. Taking into account the likely slight decrease in Pr (due to the long-term cooling, which decreases the thermal conductivity), we would expect a slightly steeper decrease in the magnetic field at the dynamo surface. In the presence of non-negligible external torques that would spin down the planet (e.g., tidal frictions with large satellites), E would (slightly) increase, which for the trends seen in our sets would lead to a (slight) enhancement of the slow magnetic decay. We leave this fine-tuned exploration for a follow-up work.

Taking Jupiter and Saturn as prototypes, we expect gasgiant dynamos to be dipolar-dominated at an advanced evolutionary age. However, according to our results, multipolar gas-giant dynamos could exist in the early stages of planetary dynamos, which would then evolve into a dipolar regime. Another possibility is that gas-giant dynamos are already born in a dipole-dominated parameter space region and remain so. These arguments can also be applied to mildly irradiated gas giants and brown dwarfs because they have similar low-Ro dynamos that follow the scaling laws of Christensen et al. (2009) because the orbital distance is large enough to prevent tidal interactions or inflation from dominating their energy budget. In contrast, these trends cannot be applied to rapidly rotating main-sequence stars as their evolution is dictated by hydrogen burning and not by a slow cooling.

Finally, the NASA exoplanet archive currently (March 2025) lists more than 650 confirmed gas-giant candidates with M sin (i) > 0.2 M_J and P_orb > 200 days, for which irradiation and tidal synchronization are negligible. Only four exoplanets under 10 pc have an estimate for the Solar System age. The closest and youngest exoplanet (0.6 ± 0.2 Gyr) is ϵ Eridani b (Hatzes et al. 2000), with a mass of ${0.66}_{-0.09}^{0.12}{M}_J$, which is a good candidate for an intense multipolar-dominated dynamo. The other three candidates, Gliese 832 b (Bailey et al. 2009), HD 219134 h (Motalebi et al. 2015; Vogt et al. 2015), and GJ 3512 b (Morales et al. 2019), with predicted ages >5 Gyr, are likely to instead host Jupiter-like dipolar dynamos. If their magnetic field is intense enough to produce electron-cyclotron maser emission that is detectable from ground (with an associated gyro-frequency, v ≃ 2.8 B[G] MHz, higher than the ionospheric ≈10 MHz cutoff; Zarka 1998), current (LOFAR) and next-generation (SKA-low) low-frequency radio interferometers might eventually detect their exoplanetary radio emission. This might provide indications of the intensity of the magnetic fields, and, possibly, their morphology.

Acknowledgements

AEL and CSG’s work has been carried out within the framework of the doctoral program in Physics of the Universitat Autònoma de Barcelona. AEL, DV, CSG, and TA are supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (ERC Starting Grant "IMAGINE" No. 948582, PI: DV). FDS acknowledges support from a Marie Curie Action of the European Union (Grant agreement 101030103). AEL, FDS, DV, CSG, and TA acknowledge the support from the "María de Maeztu" award to the Institut de Ciències de l’Espai (CEX2020-001058-M). We acknowledge the use of the MareNostrum BSC supercomputer of the Spanish Supercomputing Network, via projects RES/BSC Call AECT-2024-2-0011 (PI AEL), AECT-2023-2-0034 (PI FDS), AECT-2024-2-0003 (PI FDS). AEL acknowledges support and hospitality from the Simons Foundation through the predoctoral program at the Center for Computational Astrophysics, Flatiron Institute. We are grateful to Thomas Gastine, for guidance on the use of MagIC. The authors acknowledge insightful comments by the anonymous referee which helped to significantly improve the manuscript.

Appendix A Sensitivity on the density ratio

To evaluate how the external cut applied to the 1D MESA profile influences the overall dynamo behavior, we compare the dynamo corresponding to N_ρ ≈ 1.1, 3.0, 3.7, 4.6 (i.e., ρ_o/ρ_i ≈ 10, 20, 40, 100) for the 1 M_J 1 Gyr model keeping the same 3D resolution of (N_r, N_θ, N_ϕ) = (289, 256, 512). By looking at the spectra, the model with N_ρ ≈ 4.6 seems slightly under-resolved (there is an overall drop of only 1 order of magnitude, less than what indicates a large enough grid), but the overall quantities seem to follow the same trend as with the other distinct N_ρ models. We could not explore properly higher values of N_ρ ≳ 5.3, due to the excessive resolution required. In Table A.1 we show the different diagnostics.

Generally, many of the dimensionless quantities are affected by the value of N_ρ. Except for the N_ρ ≈ 100 model, the overall kinetic energy seems to plateau (in code units). However, since a dominant fraction of the kinetic energy is located in the nonconductive outer layers (see Fig. 4), all magnitudes containing u_rms may differ substantially. As we capture more density ratio, the aforementioned zonal flows gain importance and compete against the magnetic field in the interior affecting the overall dynamics. For example, even though Ra increases, the total buoyant power, P, decreases because it depends only on u_r. Similarly, Rm and Ro also decrease even though they depend on u_rms. This is due to the decreasing dimensionless conductivity $1/\tilde{\lambda }$, which erases the zonal flow contribution and only captures the suppression of the internal convection with higher N_ρ. Consequently, the Elsasser number Λ and E_mag, also N_ρ also decrease for the same suppression reasons. Correcting factors for the dimensionless mass M and volume V do not mitigate the differences among the N_ρ series.

To overcome this systematic effect, we tend to compare models with a similar N_ρ, independently of the mass and age of the model. For most models, we choose N_ρ ≈ 3, as it is more computationally feasible but still has a relevant nonconductive outer layer where the zonal jet develops. This restriction allows us to analyze the different saturated models for an evolutionary sequence.

Appendix B Initial conditions from previous models

To reduce computational resources and avoid starting each model with random initial conditions, we have usually employed already saturated solutions as initial conditions for other models. As we are interpreting the different models as stages in planetary evolution an obvious choice would be to use, for example, the saturated state of the 1 M_J 0.5 Gyr model as initial conditions for the M_J 0.7 Gy, and so on.

As a proof of concept, we made two tests: we used the final saturated state of a 0.5 Gyr model as initial conditions for a 1 Gyr model, and the same for a 1 Gyr and 10 Gyr models. In Fig. B.1 we show the kinetic and magnetic energy time series for this transitions. The steady states reached are indistinguishable, that is one cannot discern from the spectra, radial distribution, or final diagnostics in which initial conditions were used. The only quantities that showed a noticeable difference are the dipolarity indicators (f_dip), but they are within one standard deviation from each other. To obtain satisfactorily similar means much longer computing times would be probably required. Overall, there is convergence, starting from different initial conditions.

The activation of convection, followed by the dynamo kinematic phase and finally to the saturated phase where Lorentz forces become relevant is a lengthy computation. Starting from an already saturated solution not too far from the expected one highly reduces by more than a factor of 5 the computing time needed to reach the new steady state. Thus in most of all models, we have used a high Ra model as initial conditions, specifically the saturated 1 M_J 0.5 Gyr model with Pm = Pr = 1.

Appendix C Spectral-radial distributions and dynamo surface definition

To compare our results with (Tsang & Jones 2020), we obtained the spectral-radial energy distribution. In the left and middle panels of Fig. C.1 we show these spectra for the saturated 1 M_J 2.1 Gyr model. The radial integration of both spectra is already shown in Fig. 4. It can be observed that the radially dependent magnetic energy spectra, F_l(r), decay rapidly for r > r_m = 0.88 R_J. At the same region, the kinetic spectra start showing the equatorial jet pattern.

In the region where J = 0, we define the scalar potential V, for which B = −∇ V, with its usual spherical expansion: $$V\equiv a\sum _{l=1}^{l_{max}}{\left(\frac{a}{r}\right)}^{l+1}\sum _{m=0}^l{P}_l^m(cos\theta )[g_l^{m}cos(m\varphi )+h_l^{m}sin(m\varphi )],$$(C.1)

where a is taken as the planetary radius, $P_l^m$ are the Schmidt semi-normalized associated Legendre polynomials, and the Gauss coefficients $g_l^m$ and $h_l^m$ are obtained from measurements. Then the Lowes spectrum (Mauersberger 1956; Lowes 1974) on the planetary surface is defined as: $$R_l(a)\equiv (l+1)\sum _{m=0}^l[(g_l^m{)}^2+(h_l^m{)}^2].$$(C.2)

For any other radii, each degree l has a different contribution as they decay differently with distance r: $$R_l(r)={\left(\frac{a}{r}\right)}^{2l+4}{R}_l(a).$$(C.3)

This expression gives the Lowes spectrum in the interior of the planet, that is, the downward extrapolation which would be valid in the absence of electrical current (potential field). Therefore, we expect to have Lowes spectra only outside the dynamo surface, since inside it B stops being potential. Moreover, from equipartition reasons, it is usually assumed that the magnetic field in the dynamo region is equally distributing among different scales until the diffusive scale (‘white source hypothesis’, Backus et al. 1996). A flat spectrum at the dynamo surface would imply that at the planetary surface, there will be a linear relation log₁₀ R_l(r) ≈−β(r). For our saturated dynamo solutions, we define the Lowes spectrum as the magnetic spectra at the outermost radius, where we impose potential boundary conditions: $$R_l(r_o)=F_l(r_o).$$

On the right panel of Fig. C.1 we show both R_l(r) and F_l(r) at different radii for one specific saturated solution. The black line corresponds to r_o and the others are at different depths. At some specific depth R_l(r) stops being similar to F_l(r) and even has an unphysical negative slope. The Lowes radius is defined where β(r) = 0 and it is usually taken as the depth where the dynamo starts. As (Tsang & Jones 2020) noted, R_l(r) is already quite different from F_l(r) at such radius. To surpass this discrepancy, they similarly defined the spectra slope for F_l(r), that is, log₁₀F_l(r) ≈−α(r). The slope α(r) is almost the same as β(r) outside the dynamo region but becomes more or less flat at positive values inside. They argue that radius where α(r) and β(r) show a discrepancy is where the effective dynamo surface is.

We similarly obtained the spectra slopes α(r) and β(r) by a least square minimization between multipoles 10 and 50, which is the region where spectra are exponential and have not reached the dissipation scales set by our resolution. In Fig. C.2 we show the resulting spectra slopes, for the same models shown in Fig. 4. For some of our models, we do not obtain a flat spectrum inside, but if we take the average value of α(r) in the dynamo region and interpolate it for the decaying outer part, we obtain values very similar to the ones of r_m, shown as vertical lines. Therefore, we can safely assume that the radial position where the conductivity starts the exponential decay, r_m, is a good definition as the dynamo surface for our models.

Appendix D Diagnostics for all the models

In the Table D.1 we show the output parameters graphically shown in Figs. 5, 9.

References

---

All Tables

All Figures