Dynamo Simulations of Jupiter's Magnetic Field: The Role of Stable Stratification and a Dilute Core

Understanding Jupiter's present‐day interior structure and dynamics is key to constraining planetary accretion models. In particular, the extent of stable stratification (i.e., non‐convective regions) in the planet strongly influences long‐term cooling processes, and may record primordial heavy element gradients from early in a planet's formation. Because the Galileo entry probe measured a subsolar helium abundance, Jupiter interior models often invoke an outer stably stratified region due to helium rain. Additionally, Juno gravity data suggest a deeper, potentially stratified dilute core extending halfway through the planet. However, fits to Jupiter's gravitational data are non‐unique, and outstanding uncertainty over the equations of state for hydrogen and helium remain. Here, we use high‐resolution numerical magnetohydrodynamic simulations of Jupiter's magnetic field to place constraints on the extent of stable stratification within the planet. We find that compared to traditional interior models, an upper stably stratified layer between 0.9 and 0.95 Jupiter radii (RJ) helps to explain both Jupiter's dipolar magnetic field and zonal winds. In contrast, an extended dilute core that is entirely stably stratified (no convective layers) yields significantly worse fits to both. However, our models with extended deep stratification still generate dipolar magnetic fields if an upper stratified region is also present. Overall, we find that a planet with a dilute core i.e., strongly stably stratified is increasingly challenging to reconcile with Jupiter's magnetic field and winds. Thus if a dilute core is present, alternative modalities such as a fully convective dilute core, a complex multilayered interior structure, or double diffusive convection may be required.

The possible existence and dynamics of an inhomogenous, stratified dilute core are challenging to reconcile with traditional planetary accretion theories . Conceivably, a dilute core could either be a relic of formation (from initial mixing during formation or late-falling heavy elements), or due to erosion of the central core during the planet's subsequent 4 billion years of evolution. However, formation and evolution models show that either there is not enough initial energy to start dissolving the core at all, or there is too much energy and the heavy element gradients are fully erased within a few million years of formation (Müller et al., 2020). Explaining the contrast between in situ planetary formation models and present-day gravity analysis may require turning to alternative scenarios such as a giant impact (Liu et al., 2019) after Jupiter has cooled substantially, or migration from a colder, heavy-element rich location at 10+ AU (Shibata & Helled, 2022). If a dilute core is indeed present, the exact density gradients within it would provide key constraints on Jupiter's primordial accretion and subsequent evolution. Thus, finding additional methods to constrain the density gradients within Jupiter's interior is essential to understanding the formation of our solar system.
In addition to the implications for Jupiter's core, the physics of hydrogen and helium at high pressures may lead to additional regions of stable stratification within Jupiter. Specifically, as hydrogen transitions to a liquid metal near the surface of the planet (∼0.85-0.9 R J , where 1 R J is Jupiter's radius), helium may become briefly insoluble and rain out of solution (Stevenson & Salpeter, 1977). Physically, the heavier helium droplets sink inwards (releasing gravitational energy), leaving behind a lighter layer of helium-depleted hydrogen. Depending upon the amount of depletion, stably stratified region could develop, which would inhibit heat transfer from the deeper interior. However, there is significant disagreement regarding the phase diagram of H/He immiscibility (Brygoo et al., 2021;Fortney & Hubbard, 2003;Lorenzen et al., 2011;Morales et al., 2013;Schöttler & Redmer, 2018;Stevenson & Salpeter, 1977), though these predictions broadly suggest it would occur at 1 Mbar or deeper, that is, <0.86 RJ (Hubbard & Militzer, 2016;Nettelmann et al., 2015) (see Figure 1). Still, the uncertainties in the H/He phase diagram map onto a wide range of predictions for the onset timing and extent of helium rainout . Searching for and constraining the extent of stratification inside of Jupiter is thus critical to modeling the planet's evolution-and ultimately, the processes driving its original formation.
Finally, Juno gravity data has highlighted one last major puzzle for stratification within Jupiter. While it has long been known that interior models require atmospheric heavy element abundances that are lower than protosolar values unless ad hoc changes to the H/He equation of state are introduced (Hubbard & Militzer, 2016), this issue has become particularly apparent with the increasing precision of the Juno gravity data and their coupling to wind models. Perturbations of the H/He-adiabat (Nettelmann et al., 2021;Wahl et al., 2017), superadiabatic or double diffusive regions (Debras & Chabrier, 2019), or artificially increased surface temperatures (Miguel et al., 2022;Nettelmann et al., 2021) have been invoked. Those models commonly suggest that the outermost 10.1029/2022JE007479 3 of 20 region of Jupiter is not simply adiabatic and convective. Independently Christensen et al. (2020), suggest the presence of a convectively stable shell around 0.95-0.97 R J in order to explain the attenuation of the zonal wind speeds there. However, this stable region is located shallower (closer to the planetary surface) than the region commonly invoked for helium rain dissolution. Thus the source for a stably stratified region out at ∼0.95 R J remains an outstanding question.

Constraints From Jupiter's Magnetic Field
Above, we highlighted ongoing challenges between formation models and in situ gravity analysis. In this paper, we attempt to resolve these open questions by drawing from a complementary discipline: the study of planetary magnetic fields (dynamos). These fields are highly sensitive to many unique aspects of a planet's interior structure; in particular, the morphology and time-variation of planetary dynamos are very sensitive to the flow dynamics within a planet-including which regions are convective versus which are stably stratified. Thus, planetary magnetic field studies are particularly well-suited to investigating the presence and extent of stably stratified layers (SSLs) within Jupiter.
While Jupiter's magnetic field is primarily dipolar, the NASA Juno Mission has revealed additional complexity that is challenging to explain. Specifically, the field has strong north-south asymmetry (Moore et al., 2018). It is also dominated by a small number of regions with highly intense, localized field features such as the Great Blue Spot Moore et al., 2017) which may be sheared by the planet's zonal winds over time (Connerney et al., 2021;Moore et al., 2019). Previous dynamo models featuring an interior structure consisting of a solid inner core and a fully convective hydrogen shell tend to result in magnetic fields that are either too dipolar or too multipolar (Dietrich & Jones, 2018;Duarte et al., 2018;Jones, 2014). This is inconsistent with the intermediate-level complexity in Jupiter's magnetic field as observed by Juno. Overall, this suggests the need for alternative interior models, such as those with SSLs due to a dilute core (at the center of planet) and/or helium-depleted hydrogen near the H/He phase transition region (shallower).
SSLs have been shown to significantly impact the morphology of planetary magnetic fields in three-dimensional magnetohydrodynamic (dynamo) simulations. Specifically, a bottom liquid SSL has been shown to add complexity to magnetic fields, possibly explaining the multipolar natures of the fields of Uranus and Neptune (Stanley & Bloxham, 2004, 2006. Similarly, the possibility of upper stable stratification has been shown to be a plausible mechanism for axisymmetrizing the magnetic field of Saturn (Stanley, 2010;Stevenson, 1982;Yan & Stanley, 2021). However, these simulations were carried out in the Boussinesq approximation, which neglects the effects of compressibility. This approximation may be valid for Earth, where pressure-dependent effects may

10.1029/2022JE007479
4 of 20 only be on order of 10%. However, these effects can be more significant for giant planets like Jupiter, where the density varies by several orders of magnitude throughout the planet (See Figure 2). The capability to model SSLs was recently incorporated into the open-source dynamo code, MagIC (Gastine & Wicht, 2021); that work included a single proof of concept Jovian simulation in the anelastic (compressibility-dependent) approximation, which demonstrated that stable stratification has a strong influence on the dynamics of Jupiter's atmospheric jets and dynamo field. The specific stable layer simulated (∼0.85-0.9 R J ) resulted in a magnetic field that is too dipolar. Here we revisit this problem in more depth by testing a wider variety of layer placements, guided by hypotheses from Jovian interior structure and accretion models.

Overview
In this work, we perform state-of-the-art magnetohydrodynamic (MHD) simulations (dynamo models) of Jupiter's magnetic field morphology to constrain the location and extent of stable stratification within the planet. Specifically, we focus on two potential regions of stable stratification: (a) an upper stable layer (within the bounds of 0.8-0.95 R J ) and (b) a deep layer within 0-0.5 R J . While these features may be present in a variety of dynamical forms-fully convective, layered convection, semi-convection, and/or SSL(s)-we specifically focus on the role of stable stratification in this work, due to its profound implications for Jupiter's thermal evolution. Our goal is to constrain Jupiter's present-day interior dynamics, in order to provide an important basis of comparison with the results of accretion models. Of key interest is determining if a stably stratified dilute core needs to extend beyond 0.2 R J in order to explain Jupiter's magnetic field, which would pose a direct challenge to traditional planetary formation scenarios (Müller et al., 2020).

Methods
We seek to numerically simulate the MHD processes of Jupiter's interior, including the generation of its dynamo magnetic field and zonal winds. In particular, we model the extent of SSLs in the planet's interior, in order to investigate the role of key hypothesized features such as a stratified, inhomogeneous dilute ("fuzzy") core and a hydrogen-helium immiscibility zone. Thus, we choose a cutting-edge, open-source numerical dynamo code, MagIC (Gastine & Wicht, 2012;Jones et al., 2011), which recently added the capability to model the effects of stably stratified fluid layers (Dietrich & Wicht, 2018;Gastine & Wicht, 2021;Gastine et al., 2020). MagIC utilizes a pseudo-spectral approach with a poloidal-toroidal decomposition. The radial grid is determined by zeros of a Chebyshev polynomial. This results in a higher density of grid points near the upper and lower boundaries of the fluid shell, in order to help resolve boundary layers. We describe the numerical setup of the problem and our model parameter choices below.

Background Reference Model
MagIC models the time-evolution of the dimensionless fluid velocity, u, and magnetic field, B, inside a spherical shell with rotation rate Ω. Specifically, it solves the conservation equations for momentum (Navier-Stokes), magnetic induction, and energy. The system is represented in terms of small perturbations around a static, spherically symmetric background reference state, that is, We use the convention of tildes to represent the properties of the background state and primes to represent fluctuations throughout. The background state quantities are in general non-dimensionalized by their values at the outer boundary of the shell.
For most quantities, we use MagIC's built-in Jupiter reference model, JUP (with adjustments to ensure dipolar fields and force stable stratification, to be described). Overall this model is a polynomial fit approximation to the adiabatic Jupiter interior model of French et al. (2012), which features a homogeneous, enriched envelope out to 0.63 R J (see Figure 2). Note that the original background interior model technically has an abrupt, discontinuous jump in composition (Z and helium) and entropy, which would create an infinitesimally thin stratified region. However, these effects are smoothed over by the polynomial fit, and are usually not modeled directly. The radially varying electrical conductivity profile is super-exponential; see Figure 2. The outer boundary of our model is at 0.99 R J due to steep variations in the background profiles which are difficult to resolve numerically beyond this radius. The variation in thermal conductivity in the French et al. (2012) profiles is limited to a factor of about 60 and is thus neglected, following  and Gastine and Wicht (2021). The electrical conductivity (̃( )) profile from French et al. (2012) is shown in Figure 2f. It shows that the transition from a highly conductive metallic region to a semi-conducting fluid region takes place around r/R J ≈ 0.9. To approximate this behavior, we use a formulation introduced in Gómez-Pérez et al. (2010), Thus, the interior is approximated as a combination of a nearly constant conductivity profile in the metallic region followed by an exponential fall off. is the electrical conductivity of the fluid at the transition radius r m . The conductivity is normalized to its value at the inner boundary of the simulation. We use values of σ m = 0.2, r m = 0.9, and ξ = 13 for all our simulations.

Anelastic MHD Equations
Applying the reference state described in Section 2.1.1, we use the anelastic approximation of the MHD equations (Braginsky & Roberts, 1995;Gilman & Glatzmaier, 1981;Lantz & Fan, 1999). In this formulation, the time rate of change of density perturbations relative to the background state is assumed to be small compared to changes in the velocity-that is, sound waves are filtered out. Resolving sound waves in our models would require substantially higher time-resolution, and is not required for the goals of our study. The resulting anelastic equations are: where D/Dt = ∂/∂t + u ⋅∇ is the material derivative. S is the traceless rate-of-strain tensor: and Q ν and Q λ are the viscous and Ohmic heating terms, respectively: In the above equations, length is non-dimensionalized using the spherical shell thickness, d = r o − r i , timescale has been non-dimensionalized by the viscous diffusion time d 2 /ν, the velocity by ν/d, the magnetic field by √ Ω 0 , and the entropy perturbations by d|ds/dr| at the outer boundary of the shell. We also use the non-dimensional parameters: Ekman, Prandtl, the magnetic Prandtl number, and the dissipation number, Di, The final system of equations is solved directly for the 3D fluid velocity and magnetic fields as a function of time, given the initial background profile and a small initial random perturbation.

Stable Layer Implementation
MagIC currently provides the functionality to model planetary dynamos with SSLs. Our description of the methodology used in MagIC closely follows that of the anelastic Jovian study, Gastine and Wicht (2021). Mathematically, MagIC implements SSLs by modifying the background entropy gradient, to allow for non-adiabatic reference states (perturbations are assumed to be small to preserve the validity of the anelastic approximation). An SSL is created by setting the entropy gradient to a constant value relative to the background reference state, where Γ s is the "degree of stratification," and is set as a free parameter. Effectively, this forces stable stratification in regions where the entropy gradient, > 0 . The degree of stratification is related to the Brunt-Väisälä frequency in the stable layer by A constant dimensionless entropy gradient of = −1 is imposed in the surrounding convective regions.
The entropy gradient of the SSL is connected smoothly to the surrounding regions with a hyperbolic tangent function. The current version of MagIC allows for SSLs to be placed either at the outer boundary of a spherical fluid shell or somewhere in the interior of the planet (i.e., sandwiched between two regions of convective fluid; see Figure 3b). For an SSL extending from non-dimensional radii r strat to [r strat + h] (where h is the non-dimensional thickness), the two gradients are smoothly connected as follows: where ζ s is the stiffness of the transition between the stably stratified and convective layers.
In addition to utilizing the existing functionality for interior (sandwiched) stable layers (also referred to as "upper" layers in this paper, due to their simulated location within Jupiter), we added new functionality to the code to enable simulation of SSLs extending all the way to the bottom of the fluid shell (see Figure 3c). We parameterize this new option for a bottom stable layer similarly to the above. Here the planet is assumed to have a stable layer extending between the inner boundary and r strat . Note that whereas in the sandwich-stratification case r strat represented the bottom boundary of the layer, here it represents the upper boundary. This choice was made to minimize modifications to the shared community code.
Just as above, the bottom-stratified layer implementation is defined by a constant positive entropy gradient , and the overlying convective layer has a constant dimensionless entropy gradient of = −1 . The two gradients are smoothly connected as follows: While somewhat ad hoc, forcing stratification is standard in the solar and planetary fluid dynamic literature, and provides a controlled way to examine the effects of stratification on the resulting flows generated inside the planet. (See Takehiro and Lister (2001) and Gastine et al. (2020) for implementation in geodynamo models, and Dietrich and Wicht (2018) for an application to Saturn). We note that this approach also approximates the planet as a single-component fluid, when in reality the stratification would likely be due to chemical gradients of multiple components including H, He, and heavier elements in the planet.

Jovian Numerical Experiments
Using MagIC, we perform several numerical experiments to explore the dynamics of Jupiter's interior. Our experimental design is largely influenced by the nature of dynamo modeling, which is extremely computationally expensive. This is especially true at the resolution levels required to model the pressure-dependent interiors of giant planets. Thus, we perform a tightly focused study.
Since we can only perform a limited number of simulations due to computational resource constraints, we focus our experiments on a small number of parameters. Several previous studies have already explored the effect of various non-dimensional fluid dynamic parameters such as the Rayleigh number (Ra), Ekman number (Ek), Prandtl number (Pr), magnetic Prandtl number (Pm), and Rossby number (Ro); as well as the broader effects of the background interior models (ranging from idealized polytrope number to more "realistic" Jupiter interiors like that of French et al. (2012)); and the role of different electrical conductivity profiles (i.e., different metalization depths/drop-off rates) and heating models (internally heated vs. top-down cooling). Thus, we rely on the important conclusions those authors have derived (see Duarte et al. (2018) for a review), and do not repeat those parameter surveys here. Instead, we focus on the effect of large-scale changes to Jupiter's interior structure: the role of stable stratification due to a possible dilute core and/or helium rainout. See Table 1 for a list of our numerical dynamo simulation cases.

Baseline Control Models
To understand the effect of these new possible interior structures, we first perform a single "dipolar baseline" control model run using a traditional Jupiter interior structure (a solid compact core from 0 to 0. was specifically chosen as the maximum radius considered plausible by traditional accretion models (Müller et al., 2020). Our dipolar baseline model uses a Rayleigh number Ra = 8 × 10 7 , Ek = 10 −4 , Pr = 1, and Pm = 3. This model was based on Case 34 of (Duarte et al., 2018), but with a slightly higher Rayleigh number. We also use a grid resolution of N r = 161, N θ = 384, and N ϕ = 768 grid points in the radial, meridional, and azimuthal directions respectively, leading to a maximum spherical harmonic degree l max = 256. All of our simulations also used these same simulation control properties, with the exception of the Rayleigh number for our multipolar baseline model. We chose to keep most values the same in each simulation, and focused on varying the core radius, and the locations and thicknesses of SSLs.
We acknowledge the common concern that computational restrictions require an unphysically large Ekman number for planetary studies (i.e., too viscous). Some works have hinted at the possibility of different dominant force balances only seen at lower Ekman numbers (Yadav et al., 2016); whether this is consistently the case remains to be seen. However, stably stratified models require significantly higher grid resolutions than traditional compact core models, and so reaching these low Ekman numbers is an even greater challenge here. We choose instead to perform a broader suite of models with a higher Ekman number (Ek = 10 −4 ) that is commonly used in the literature. This represents a balanced approach, and allows us to easily compare the effects of our stably stratified models with both our baseline models and previous models produced by the scientific community.
Our baseline simulation was consistently dipolar for the entire run time period. This baseline model was then used as the initial checkpoint for our simulations with SSLs. We also perform a "multipolar baseline" control model Note. Here we provide a list of dynamo simulations in our case study. The models are of four classes (aligned with Figure 3): fully convective models ("baseline"), models featuring an upper SSL ("Upper SSL"), models with a stably stratified dilute core ("Dilute core"), and combination models with both an upper and lower stably stratified regions ("Combo DC/UP"). While the table highlights parameters that differ between models, many parameters are the same across simulations. For example, all of our models use non-dimensional fluid parameter values of Ra = 8 × 10 7 (except the multipolar baseline model which uses Ra = 1.8 × 10 8 ), Ek = 10 −4 , Pr = 1, Pm = 3, and a simulation upper cutoff radius of 0.99 R J . All simulations use a grid resolution of N r = 161, N θ = 384, and N ϕ = 768 grid points in the radial, meridional, and azimuthal axes respectively. Models including a stably stratified layer use a stratification slope of ζ s = 130. The degree of stratification, Γ s , was fixed to give a maximum N/Ω ratio of approximately 10 in the stratified layer. Here RMSLE stands for root mean square log error.

Table 1
Overview of Dynamo Simulations using identical parameters as the dipolar baseline, with the exception of a higher Rayleigh number (1.8 × 10 8 ) to yield a multipolar model.

Dynamo Models With Stratified Interiors
After creating our baseline comparison models, we perform a series of dynamo simulations featuring a stably stratified dilute core (referred to here as "Dilute Core models" throughout the text, although in reality dilute cores need not be stratified). Our models are simulated as a solid inner core from 0 to 0.2 R J , with a fluid stably stratified region overlying the solid inner core (see Figure 3c); the presence of the inner core provides our models with numerical stability. Specifically, we simulate models with total non-convective core radii (solid inner core plus liquid SSL) out to 0.25, 0.3, 0.4, and 0.5 R J . We performed exploratory simulations with bottom SSLs extending to even higher radii, but the initial results of these simulations were too similar to our other stratified dilute core models to justify the additional expenditure of computational time within our limited parameter study. Overall, our dilute core cases provide a strong contrast with our baseline cases, which only have a total core radius of 0.2 R J -the maximum radius permitted by traditional accretion models (Müller et al., 2020). As a control, we also perform a simulation with a large, entirely solid core (i.e., no SSL) ranging from 0 to 0.4 R J . We note that while we only test the case of a dilute core that is entirely non-convective here (i.e., solid or stably stratified across the entire volumetric region), additional geometries such as multilayered or fully convective dilute cores may be possible (Militzer et al., 2022).
For all of our dilute core models, the extent of the SSL was chosen first, and then the degree of stratification was set inside that region to yield a maximum Brunt-Väisälä frequency N/Ω ≈ 10 in the stratified region. This value was used so as to be directly comparable to the case by Gastine and Wicht (2021). While the true value inside Jupiter may be closer to 1-3 for upper SSLs depending upon the interior model chosen (see Miguel et al. (2016), Debras and Chabrier (2019), and Gastine and Wicht (2021) for details), the increased stiffness of the SSL helps balance out the effect of the artificially enhanced viscosity in dynamo models due to numerical limitations. We did perform initial probing studies using a more weakly or strongly stratified dilute core (N/Ω ≈ 0.5 or 100, respectively). Models with higher values of N/Ω appear to have stronger propagation of internal gravity waves in the core, which could have implications for observations of the planet's interior modes by the Juno spacecraft or telescopic observatories. However, the degree of stratification in the dilute core did not significantly impact the primary concern of our study-the morphology of the magnetic field at the planetary surface-so we do not consider it further here.
We also performed simulations featuring an upper SSL ("Upper SSL" models). Our experiments explored layers with thicknesses of 0.05, 0.1, and 0.15 R J (between 3,500 and 11,000 km in thickness), and a variety of depths (such that the layer was fully encompassed within the radial region between 0.8 and 0.95 R J ; see Figure 3b). An upper SSL may be caused by helium rainout near the outer boundary of the dynamo region. In such case, the layer thickness and depth would be directly determined by the planet's temperature, pressure, and electrical conductivity profiles. However, the H/He mixing phase diagram has sufficient uncertainty that we treat these values as free parameters, both for exploring Jupiter's magnetic field as well as gaining general insight about the effects of SSLs in giant planets. We note that an outer extent of 0.95 R J for He rainout is shallower than predicted by all current H/He phase diagrams, but we include it here as an outer bound due to suggestions from previous magnetic field analyses that it is a radius of interest (Christensen et al., 2020;Moore et al., 2019). We discuss the interesting contrast between the locations of possible stratified layers suggested by high pressure physics models of hydrogen-helium gradients versus evidence from other disciplines like seismology and magnetic field analysis more in Section 4; it is possible that multiple processes causing stable stratification may be at play.
Finally, we performed combination simulations of models including both an extended dilute core and an upper SSL (possibly due to helium rain); see Figure 3d. Such models are denoted as "Combo DC/UP" models throughout the text. The upper layers in these cases were simulated identically to those described above, while the dilute core region was modeled as a large solid inner core for simplicity. This choice is justified by our results (See Section 3), which find no significant differences in field morphology between dilute core models with a purely solid inner core versus those with a bottom SSL.
We note that we have chosen to idealize the possible dilute core within Jupiter as a stratified layer overlying a small compact core (or, for comparison, as a large and entirely solid core region), and the hypothesized helium immiscibility zone as an upper stratified layer. In reality, the dynamics of these features may be further complicated by alternative processes such as layered convection (Moore et al., 2018), marginal stability, or semiconvection/double-diffusive convection (Leconte & Chabrier, 2013;Nettelmann et al., 2015), that is, a region that is compositionally stable but thermally unstable. Consideration of these processes is beyond the scope of our current study, but would be a fruitful and important direction for future research in the field.

Results
In this study, we explore whether adding complexity to Jupiter's interior structure (i.e., a stably stratified dilute core and/or a possible helium rainout zone) could help explain the planet's observed magnetic field morphology. That is to say, we seek a "Goldilocks" or intermediate dynamo model-whereby starting with a dipolar baseline case and adding new interior features, we manage to add a small amount of necessary complexity to bridge the two end-member regimes.
We primarily focus on the broad implications of time-averaged properties (details are described below). Planetary dynamo models are turbulent and may vary strongly in time; similarly, observational data has also found evidence of time-variation in Jupiter's actual magnetic field (Moore et al., 2019). Working with averaged quantities helps to prevent biased conclusions from possibly unrepresentative "snapshots" of the field at isolated moments in time.
Resolving the time-dependence of features such as isolated flux spots like the Great Blue Spot in a statistically meaningful manner requires longer simulation run times and is computationally expensive; we hope to follow up on our promising candidate models in a future study.

Magnetic Field Morphology
First, we examine the magnetic field morphology of our numerical simulations. A planet's dynamo magnetic field at the surface can be represented as a potential field, V, with degree l and order m spherical harmonic coefficients, and : where V is a function of position relative to a planetary reference radius, a (in our case, a = 1 R J = 71,492 km).

1D Power Spectra
Using the field coefficients at the outer boundary of our fluid dynamics simulations (0.99 R J ), we find the time-averaged Lowes power spectrum (Lowes, 1966) to compare the total power, R l , at each spectral degree: The Lowes power spectrum is commonly used to estimate the location of the surface of a planetary dynamo by assuming that the dynamo source emits power as a flat spectrum, and finding the distance to this source region that best-fits the power spectrum of the observed field. This is based on geometric reasoning: the further one is from the source generating region, the more strongly higher order terms in the potential field expansion will decay relative to lower order terms (i.e., "all things resemble a dipole when observed from a far enough distance", see Equation 21). Using this method on Juno-derived magnetic field models, one would estimate that Jupiter's convective dynamo extends to between 0.81 and 0.85 R J (Connerney et al., , 2021Sharan et al., 2022).
This method provides an excellent first approximation, but is largely empirical in nature and primarily informed by studies of Earth's magnetic field. There are no first principles requirements for dynamos to emit at a flat spectrum, and the "dynamo surface" may be poorly defined for fluid planets such as Jupiter which have radially varying conductivity profiles and possibly intricate, multilayered interiors. By performing a suite of 3D MHD simulations, we can compare the data-derived power spectrum to the specific spectra resulting from a range of realistic scenarios for Jupiter's interior.
After initializing each dynamo model from our dipolar baseline checkpoint file, we allowed the simulations to stabilize for 0.9 viscous diffusion times (0.3 magnetic diffusion times for Pm = 3, similar to other works). Then we calculated the time-averaged Lowes spectrum for the remainder of the simulation period, normalized to the power in the main dipole. We note that all dipolar models in our study were relatively steady, while all multipolar models underwent frequent reversals. Figure 4 compares the resulting power spectra of our dynamo simulations (colored lines) to the JRM09 magnetic field model  derived from Juno spacecraft observations (black line), and is summarized in Table 1. We also show Earth's power spectrum using the IGRF-13 model for comparison (Alken et al., 2021). Earth's power spectrum has a much faster falloff with degree than Jupiter's due to the deeper depth of its dynamo region relative to planetary radius (∼0.55 R E vs. ∼0.85-0.95 R J ). Our dipolar and multipolar baseline models reproduce the two end-member regimes observed in the literature for compact core models, as expected. In particular, our choice of dipolar baseline model is a surprisingly reasonable fit to the power spectrum of JRM09, although it has a significant mismatch in the power of the quadrupolar term.

Dilute Core Models
Like traditional compact core models, models including SSLs still appear to fall into the end-member extremes of highly dipolar or highly multipolar morphologies. Past dynamo modeling studies of Uranus and Neptune's interiors in the Boussinesq limit (Stanley & Bloxham, 2006) had found that larger inner core sizes lead to increased field complexity; models with intermediate shell thickness somewhat resembled Jupiter's field. Thus, we had initially expected that the inclusion of a non-convective dilute core region of the correct radius could lead to a case with intermediate morphology, but our results do not support that hypothesis. Instead, our non-convective dilute core models unanimously yield multipolar magnetic field morphologies (red/orange/yellow lines in Figure 4). This was true for all radial extents of the non-convective dilute core examined here, with larger dilute cores generally having worse root mean square log errors to the JRM09 power spectrum than smaller ones (See Table 1). Note that we did not consider interior structures with convective dynamics such as multi-layered dilute cores or double diffusive convection.
While previous studies found that dynamo models in the Boussinesq approximation are sensitive to the difference between a liquid stably stratified core and a solid inner core of the same radius (Stanley & Bloxham, 2004, 2006, no significant difference in field morphology was observed in our anelastic simulations. This is likely because the large number of density scale heights traversed in our anelastic Jupiter simulations reduces the models' overall sensitivity to the specific nature of the deep bottom boundary, aside from the general aspect ratio of the convective region. Overall we find that a fully stably stratified dilute core alone does not appear to explain Jupiter's unusual magnetic field (Moore et al., 2018).  . All spectra have been normalized to the power in the dipolar term. We also show the power spectrum of Earth's dynamo field using the IGRF-13 model for comparison (Alken et al., 2021).

Upper SSL Models
In contrast to the multipolar dilute core models, our models implementing an upper stably stratified region occurring somewhere between 0.8 and 0.95 R J (Figure 3b) are primarily dipolar (Figure 4b). We find that the Lowes power spectrum appears to be strongly correlated with the location of the stable layer: shallower layers have stronger contributions from higher order terms. This trend can be explained by the stable layer's ability to change the effective "depth of the dynamo surface." In MHD simulations of Jupiter with a fully convective interior, dipolar dynamos are confined to deeper regions of the planet due to competition with the zonal winds (while multipolar dynamos are generated in more shallow regions) (Duarte et al., 2018). Here we find the stable layer provides a buffering separation zone between Jupiter's convective dynamo and its strong zonal winds, allowing for a dipolar dynamo to extend to shallower depths than might otherwise occur in the presence of strong zonal winds.
Unlike the layer placement, changing the thickness of an upper stable layer due to helium rain did not appear to affect the Lowes power spectrum significantly (see Figure 4b). Layers from 0.8 to 0.85, 0.8 to 0.9, and 0.8 to 0.95 R J all have nearly identical power spectra, as do the layers from 0.85 to 0.9 and 0.85 to 0.95 R J . Overall we find that the location of the bottom of the stratified layer-and not the layer thickness-appears to be the primary determining factor of the dynamo's power content. This is consistent with interpreting the power spectrum in terms of geometric distance from the primary source of field generation: the convective dynamo region. We also note that from the perspective of the magnetic Reynolds number, the bottom of an upper stable layer is likely where advective shearing of the electromagnetic field would be strongest due to the increased electrical conductivity, although this effect would also depend upon the rate of decay of fluid velocity with depth.

Combination Models
While all of our original dilute core models are multipolar, combination models that have both an extended dilute core region and an upper stable layer (Figure 3d) produce both dipolar and multipolar morphologies, depending upon layer placement (Figure 4c). As a reminder, our previous results found no significant difference between a large, extended core that is entirely solid (0-0.4 R J ) versus a core consisting of a smaller solid inner core with an overlying liquid stably stratified core (0-0.2 and 0.2-0.4 R J , respectively). Thus for numerical convenience, our combination models simulate the dilute core region as a larger solid core extending from 0 to 0.4 R J . We find that combination models with dipolar magnetic fields closely reproduced the original power spectrum of these models when only an upper stable layer was included (see Figures 4b-4c, models with stable layers at 0.85-0.9 and 0.8-0.95 R J ). However, other combination models were entirely multipolar. This likely represents that there is a close tipping point between dipolar and multipolar dynamo modes (as found in other studies, e.g., Duarte et al. (2018)), and that the exact placement and thickness of the stable layer(s) can push this balance one way or the other, depending upon the resulting convective aspect ratios. While bistability was not explicitly observed in any of our combination model simulations, we cannot fully rule it out either.

Best-Fitting Model
Our best-fitting model has an upper SSL extending between 0.9 and 0.95 R J , and a traditional compact core from 0 to 0.2 R J (Figure 3b, not to scale; see also the medium blue solid line in Figure 4b). It provides an excellent fit to the normalized JRM09 power spectrum, with near exact fits to the quadrupole and octupole terms. The fit of our upper SSL model is significantly better than our baseline dipolar compact core model (root mean square log error, or RMSLE, of 0.41 vs. 0.70). More generally, all upper SSL models with a bottom boundary starting at 0.85 R J or shallower provide a better fit to JRM09 than the dipolar baseline model. The best-fitting model's power spectrum is smoother than JRM09, which can be explained by the time-averaged nature of our model (since JRM09 only represents the Jovian field at a specific snapshot in time).
Interestingly, our best-fitting model is also bistable. That is, it stochastically switches between a non-reversing dipolar and a reversing multipolar regime (see Figures 4 and 5). While atypical for Earth-like models, bistability is a commonly observed phenomenon in anelastic dynamo models of giant planets (Duarte et al., 2018). This is usually thought to be due to competition between the planet's dynamo magnetic field and the zonal winds. We find here, however, that this bistable competition is still observed on rare occasions even when a SSL is added as a buffer between the dynamo and the zonal winds.

Zonal Wind Profiles
Explaining Jupiter's atmospheric zonal wind profile alone has proved to be a decades long challenge in geophysical fluid dynamics and atmospheric science. Producing a Jovian model that has both reasonable winds and a dipolar dynamo simultaneously is an even larger challenge, due to competition between the winds and the dynamo. While the primary goal of our study is to understand the role of stable layers on the morphology of Jupiter's magnetic field, our fluid dynamic simulations also provide important insights on this secondary indicator, since our models encompass almost the full radial extent of the planet.
Previous dynamo studies implementing a compact core have typically found it challenging to simultaneously produce a dipole dominated surface field and a Jupiter-like profile of alternating zonal winds (Dietrich & Jones, 2018;Duarte et al., 2013Duarte et al., , 2018Heimpel & Gómez Pérez, 2011;Jones, 2014). Our present study is in agreement with recent studies which have suggested that having a SSL between the highly conducting and non-conducting regions of a simulation helps in overcoming this obstacle (Christensen et al., 2020;Gastine & Wicht, 2021).
We find that models featuring thin, upper stable layers (possibly due to mechanisms such as helium rain) provide largely improved fits to Jupiter's time-averaged zonal winds (See Figures 6h-6k). The width of the prograde jet To separate these two modes more quantitatively, we apply a cutoff dipolarity value of 0.5 (orange line), where high dipolarity signifies a dipolar state. As with our other models, the first 0.9 viscous diffusion times of the simulation are not used in our analysis in order to allow the model to adjust from its dipolar baseline initialization state. (b) The transitions between dipolar and multipolar states are clearly demonstrated by the dynamo's pattern of polarity reversals. Our dipolar Jovian models consistently have steady magnetic polarity, while the multipolar models undergo frequent reversals.
shrinks slightly in latitude to more closely match Jupiter's observed winds, and some models even start to replicate the weak retrograde jets near 20-30° latitude. While the wind morphology of combination models (featuring both a dilute core and upper stable layer) is slightly more varied across models, they provide similarly good fits to the observed winds (see Figures 6n-6s). This again suggests that an upper stable layer helps offset the issues caused by a stratified dilute core.
When the upper SSL is made too thick, however (e.g., 0.1-0.15 R J ; Figures 6k-6m), an unrealistically strong high latitude jet appears to dominate. At the extremes, the models become highly time-varying (large standard deviation), and the equatorial prograde jet is lost entirely. It is possible that these effects are an artifact of treating the stable layer thickness as a free parameter, when in reality it is a function of the conductivity profile. While computational limitations enforced a small parameter suite on our study, this is a rich topic for a follow-up exploration future work.

Discussion
In this work, we performed advanced numerical MHD simulations of Jupiter's interior, investigating the existence and extent of stably stratified regions within the planet. Specifically, we investigated stratified regions representing the effects of a possible dilute core and/or an upper stably stratified region-features postulated Figure 6. Time-averaged surface zonal winds. We compare the normalized, time-averaged zonal winds of our dynamo models at the outer boundary of the fluid shell, 0.99 R J (colored lines). We also show the standard deviation of the normalized winds (gray shaded region), and cloud-tracking observations of Jupiter's surface zonal winds from December 2016 (approximately coincident with Juno perijove 3; thick black line) (Tollefson et al., 2017).
by gravity observations and high pressure physics studies, and key to testing hypotheses of Jupiter's formation. Previous fully convective dynamo studies with only a compact core display a strong dichotomy of too-dipolar or too-multipolar magnetic field morphologies. Our hope was that by adding the correct combination of layer placement and thicknesses to Jupiter's interior, we would identify a dynamo model that better explains field models produced from Juno spacecraft observations. Indeed, we find that an upper SSL does help to explain Jupiter's observed features. Such layers produced dipolar field morphologies; layers with bottoms at 0.85 R J or shallower produced better fits to Jupiter's observed power spectrum than our compact core baseline model, with the shallowest layers having the best fit of all. These models also provide slightly enhanced fits to the planet's zonal winds. Our results are consistent with trends found from dynamo explorations of Jupiter's possible conductivity profile with a compact core (Duarte et al., 2018), which also favor shallower metalization depths. Since the dynamics of hydrogen-helium immiscibility and resulting density gradients/stratified layers are intimately linked to the metalization of hydrogen, this paints a cohesive picture of the planet's interior structure and dynamics.
In contrast, we find that an extended, stably stratified dilute core alone does not appear to explain Jupiter's field. Our stably stratified dilute core models only produced multipolar magnetic fields, and yielded zonal wind profiles with too thick of a prograde jet. Thus, we find that stratified dilute core models (without any additional layering) are outmatched by fully convective/compact core models in all respects, which are usually able to somewhat match Jupiter's observed winds or observed field but not both. We also find that whether the "stratified" dilute core is fully solid or a stably stratified liquid does not appear to make a difference in the anelastic dynamo limit; while we present stratified dilute core simulations out to 0.5 R J here as part of our tightly focused case study, we also performed shorter exploratory studies out to higher radii and still found the results to be consistently multipolar.
When a stratified dilute core is present in combination with an upper SSL, however, these problems are partially resolved. While most combination models were multipolar, dipolar cases could still be found. Additionally, the planet's zonal wind profiles improved as well, yielding thinner equatorial jets more consistent with Jupiter's ground-truth.
We caution that while our evidence is strongly against the specific interior model of a stably stratified dilute core that is entirely non-convective and with no upper stable layer, more research will be required to fully rule this option out. Our study only explored a limited parameter space in conductivity and non-dimensional fluid dynamics parameters due to computational limits. Some studies have suggested that the thickness of the prograde jet may decrease at smaller Ekman numbers . This is unlikely to resolve the mismatch in magnetic field morphology, but it could improve the fit to Jupiter's zonal winds.
Additionally, while we assumed that Jupiter's observed surface zonal winds from cloud tracking extended deeper into the planet (i.e., to at least 0.99 R J where our fluid shell began), this may not be the case. Historically, fluid dynamical theories of Jupiter's zonal winds have fallen into two main camps-(a) that the winds extend geostrophically deep into the planet along so-called Taylor Columns (Busse, 1976), or (b) that the surface winds could be considered a shallower "weather layer" due to effects such as variable solar insolation or latent heat from ortho/para hydrogen phase transitions (see Vasavada and Showman (2005) for a review). Models of Jupiter's interior and zonal winds from Juno gravity observations suggest that the winds do extend deep (Guillot et al., 2018;Kaspi et al., 2018), rendering this latter hypothesis unlikely.
Finally, while our results cast doubt on the considered the case of a dilute core that is stably stratified across the entire region, many other flavors of a dilute core may still be possible. One possibility is a marginally stable dilute core. While we ran initial exploratory simulations and found the field morphology to be much more sensitive to the overall aspect ratio of the convective region than the stability of the bottom boundary, we cannot rule out the possibility of a marginally stable dilute core entirely and recommend it worthy of additional consideration. Other possibilities include a fully convective (i.e., homogeneous) dilute core; layered convection (e.g., featuring an outer stable region where the gradient is sharpest overlying a second interior convective region Militzer et al. (2022)), or semiconvection (Leconte & Chabrier, 2012;Moll et al., 2017). Indeed, while semi-convection traditionally operates as a staircase of hundreds to thousands of thin layers for Earth-like regimes, the dynamics of semi-convection in astrophysical regimes are poorly understood and may deviate substantially from this model (Garaud, 2018). It is clear that an investigation of the possibility of magnetic field generation under semi-convection is required. Despite these caveats, however, our initial evidence from Jupiter's magnetic field and zonal winds is still starkly against an interior profile that has only a dilute core and no upper stratified layer.

Comparisons With Current Jupiter Interior Models
A central question of this work was to determine if a stably stratified dilute core could extend beyond 0.2 R J in Jupiter. This is a critical metric in planetary formation modeling. Müller et al. (2020) found that traditional accretion models cannot produce models with dilute cores that have radial extents greater than 0.2 R J , even when the models were intentionally pushed to possibly nonphysical limits in order to support the formation of an extended dilute core. That radius has proved to be a key cutoff for Jupiter's dynamo magnetic field as well; stratified dilute core radii larger than 0.2 R J did not produce any dipolar (i.e., Jupiter-like) magnetic fields, except for when an upper SSL was also present. Thus, our results are consistent with current planetary accretion models, including new work suggestive of an interior with both a fuzzy core with an upper stable layer (Shibata & Helled, 2022).
Closer to the planet's surface, interior models constrained by Juno gravity data often require extremely low concentrations of atmospheric heavy elements (Z). One way to solve this problem is to artificially enhance the planet's surface temperature (Miguel et al., 2022;Nettelmann et al., 2021), or to assume a steeper-than-adiabatic T-gradient in an inhomogeneous stably stratified region at ∼0.95 R J (Debras & Chabrier, 2019). Our study reports evidence from simulations of Jupiter's magnetic field morphology that strongly support an upper stratified layer in the planet's interior at this depth. The synergistic results from gravity and magnetic field analyses reported here is extremely promising. However, the exact mechanism of such a shallow SSL (0.9-0.95 R J ) is currently unknown.
One possibility worthy of consideration is helium rainout. Outstanding questions remain about the existence and extent of helium rainout within the giant planets. For example, while Saturn is currently believed to be experiencing helium rain, there is significant uncertainty regarding the extent of this process. Evolutionary models predict that Saturn should be colder (cooled more quickly since formation) than its present day temperature, and so a process such as helium rain can help to explain Saturn's slower cooling and excess luminosity (Fortney & Hubbard, 2003). Analysis of Saturnian upper atmospheric helium abundance observations vary depending on the adopted methods, from a slightly depleted state (Conrath & Gautier, 2000;Koskinen & Guerlet, 2018) (minimal helium rainout) to severely depleted (Achterberg & Flasar, 2020;Conrath et al., 1984;Sromovsky et al., 2016) (extensive helium rainout). The case for helium rain is less certain for Jupiter based on the planet's pressure temperature profile (Brygoo et al., 2021;Morales et al., 2013;Schöttler & Redmer, 2018), although atmospheric measurements of outer atmospheric helium and neon depletion from the Galileo probe strongly suggest the presence of helium rain (Niemann et al., 1998).
But our preferred depth for the upper SSL, 0.9-0.95 R J , is much shallower than current phase diagrams of H/ He mixtures suggest helium rain should occur. While there are large uncertainties in the exact location of H/ He demixing, models consistently find it should be in the vicinity of 0.6-0.86 R J , or 0.9-10 Mbar (Hubbard & Militzer, 2016;Nettelmann et al., 2015); see Figure 1. What's more, even if helium rain is present, it may not lead to stable stratification in the rainout region. Rather, it could be fully convective, or lead to more complicated dynamics. For example, Debras and Chabrier (2019) hypothesize helium rain may cause a double diffusive instability in the vicinity of the upper SSL of our best-fitting model. That is to say, the existence of a stable layer at 0.95 R J would imply that either our current understanding of the thermodynamics of hydrogen-helium at high pressures would need revision, or that an effect other than the helium gradient may be responsible for the layer. Thus, the mechanism for a stable layer at approximately 0.95 R J is still an open question.

Conclusions
Numerical dynamo models of giant planet interiors provide insights not only into a planet's present day interior structure and dynamics, but also its past formation and evolution. Our new, state-of-the-art simulations featuring SSLs suggest that an upper SSL helps to explain Jupiter's dynamo and winds; the best-fitting model used featured a SSL from 0.90 to 0.95 R J . In contrast, our models with stably stratified dilute cores provided worse fits to Jupiter's observed field and winds, unless an upper SSL was also included. Our work highlights the ongoing challenges in explaining the planet's origins, and reinforces Jupiter's interdisciplinary nature. As the Juno Mission continues its extended mission, we hope that future modeling work can benefit additional constraints such as maps of the magnetic field's time-variation, and the geometry of possible wind-sheared field features.