# Juno Spacecraft Measurements of Jupiter's Gravity
## Imply a Dilute Core

Burkhard Militzer, William B. Hubbard, Sean Wahl, Jonathan I. Lunine, Eli Galanti, Yohai Kaspi, Yamila Miguel, Tristan Guillot, Kimberly M. Moore, Marzia Parisi, et al.

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

► To cite this version:

Burkhard Militzer, William B. Hubbard, Sean Wahl, Jonathan I. Lunine, Eli Galanti, et al., Juno Spacecraft Measurements of Jupiter's Gravity Imply a Dilute Core. The Planetary Science Journal, 2022, 3, 10.3847/PSJ/ac7ec8. insu-03869585

HAL Id: insu-03869585

https://insu.hal.science/insu-03869585

Submitted on 24 Nov 2022

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Juno Spacecraft Measurements of Jupiter's Gravity Imply a Dilute Core

Burkhard Militzer<sup>1,2</sup> <sup>①</sup>, William B. Hubbard<sup>3</sup> <sup>①</sup>, Sean Wahl<sup>1</sup> <sup>①</sup>, Jonathan I. Lunine<sup>4</sup> <sup>①</sup>, Eli Galanti<sup>5</sup> <sup>①</sup>, Yohai Kaspi<sup>5</sup>, Yamila Miguel<sup>6,7</sup>, Tristan Guillot<sup>8</sup> <sup>①</sup>, Kimberly M. Moore<sup>9</sup> <sup>①</sup>, Marzia Parisi<sup>10</sup> <sup>①</sup>, John E. P. Connerney<sup>11,12</sup>, Ravid Helled<sup>13</sup> <sup>①</sup>, Hao Cao<sup>14</sup> <sup>①</sup>, Christopher Mankovich<sup>9</sup> <sup>①</sup>, David J. Stevenson<sup>9</sup> <sup>①</sup>, Ryan S. Park<sup>10</sup> <sup>①</sup>, Mike Wong<sup>15,16</sup>, Sushil K. Atreya<sup>17</sup>, John Anderson<sup>10</sup>, and Scott J. Bolton<sup>18</sup>

<sup>1</sup>Department of Earth and Planetary Science, University of California, Berkeley, CA 94720, USA; militzer@berkeley.edu

<sup>2</sup>Department of Astronomy, University of California, Berkeley, CA 94720, USA

<sup>3</sup>Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA

<sup>4</sup>Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

<sup>5</sup>Department of Earth and Planetary Sciences, Weizmann Institute of Science, Rehovot 76100, Israel

<sup>6</sup>SRON Netherlands Institute for Space Research, Sorbonnelaan 2, NL-3584 CA Utrecht, The Netherlands

<sup>7</sup>Leiden Observatory, University of Leiden, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

<sup>8</sup>Universite Cote d Azur, OCA, Lagrange CNRS, F-06304 Nice, France

<sup>9</sup>Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, CA 91125, USA

<sup>10</sup>Jet Propulsion Laboratory, 4800 Oak Grove Drive, Pasadena, CA 91011, USA

<sup>11</sup>Space Research Corporation, Annapolis, MD 21403, USA

<sup>12</sup>NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA

<sup>13</sup>Institute for Computational Science, University of Zurich, Winterthurerstr. 190, CH-8057 Zurich, Switzerland

<sup>14</sup>Department of Earth and Planetary Sciences, Harvard University, 20 Oxford Street, Cambridge, MA 02138, USA

<sup>15</sup>Center for Integrative Planetary Science, University of California, Berkeley, CA 94720, USA

<sup>16</sup>SETI Institute, 189 Bernardo Avenue 200, Mountain View, CA 94043, USA

<sup>17</sup>University of Michigan, Climate and Space Sciences and Engineering, Ann Arbor, MI 48109, USA

<sup>18</sup>Southwest Research Institute, San Antonio, TX 78238, USA

Received 2022 May 31; revised 2022 July 1; accepted 2022 July 1; published 2022 August 4

## Abstract

The Juno spacecraft measured Jupiter's gravity field and determined the even and odd zonal harmonics, $J_n$, with unprecedented precision. However, interpreting these observations has been a challenge because it is difficult to reconcile the unexpectedly small magnitudes of the moments $J_4$ and $J_6$ with conventional interior models that assume a large, distinct core of rock and ice. Here we show that the entire set of gravity harmonics can be matched with models that assume an ab initio equation of state, wind profiles, and a dilute core of heavy elements that are distributed as far out as 63% of the planet's radius. In the core region, heavy elements are predicted to be distributed uniformly and make up only 18% by mass because of dilution with hydrogen and helium. Our models are consistent with the existence of primary and secondary dynamo layers that will help explain the complexity of the observed magnetic field.

*Unified Astronomy Thesaurus concepts: Solar system gas giant planets (1191)*

*Supporting material: machine-readable table*

## 1. Introduction

Conventional models for giant planet interiors are constructed with a compact core of rock and ice, atop which is a hydrogen–helium envelope. Since hydrogen and helium are predicted to become immiscible at megabar pressures (Stevenson & Salpeter 1977), one typically separates this envelope into an upper helium-depleted layer of molecular hydrogen, an intermediate helium rain layer, and a deep helium-enriched layer of metallic hydrogen. There is indeed good evidence that helium rain has occurred on Jupiter because the Galileo entry probe measured a helium mass fraction of $\bar{Y} \equiv Y/(X+Y) = 0.238 \pm 0.005$ (von Zahn et al. 1998) that is well below the protosolar value of 0.2777 (Lodders 2010). Furthermore, neon in Jupiter's atmosphere was measured to be ninefold depleted relative to solar, and this can be attributed to efficient dissolution in helium droplets (Roulston & Stevenson 1995; Wilson & Militzer 2010). However, many details of

layering in giant planet interiors have remained uncertain. Since a detailed experimental characterization of H–He phase separation is still outstanding, one relies instead on predictions from ab initio computer simulations (Morales et al. 2013) to constrain the thickness of the helium rain layer. Also, the abundance of elements heavier than hydrogen and helium is poorly constrained in all but the uppermost layer of Jupiter's atmosphere. The Galileo entry probe measured the heavy-element abundances up to a pressure of 22 bars. There, noble gases and several other heavy elements were found to have approximately three times the protosolar concentration (Mahaffy et al. 2000; Wong et al. 2004). In contrast, the probe measured subsolar concentrations of oxygen, the element with the largest mass contribution due to its high relative solar abundance. Prior to Juno, it was debated whether the subsolar oxygen is representative of the average envelope or if it reflects inhomogeneities associated with dynamical processes localized to the 5 µm hot spot into which the probe fell.

The Juno spacecraft, in orbit about Jupiter since 2016, improved our knowledge of the gravity field through multiple close flybys. These measurements are summarized in terms of

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
<table><caption>**Table 1**
Comparison of Juno Measurements (Durante et al. 2020) and Predictions of Model A for the Even and Odd Gravity Coefficients</caption><thead><tr><th></th><th>J<sub>2</sub></th><th>J<sub>4</sub></th><th>J<sub>6</sub></th><th>J<sub>8</sub></th><th>J<sub>10</sub></th></tr></thead><tbody><tr><td>Even J<sub>n</sub> × 10<sup>6</sup></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Interior model</td><td>14,696.4484</td><td>-586.8463</td><td>34.4692</td><td>-2.4983</td><td>0.2067</td></tr><tr><td>Wind model (i)</td><td>0.0579</td><td>0.2377</td><td>-0.2684</td><td>0.0763</td><td>-0.0231</td></tr><tr><td>Interior+wind</td><td>14,696.5063</td><td>-586.6086</td><td>34.2008</td><td>-2.4220</td><td>0.1837</td></tr><tr><td>Juno measurement</td><td>14,696.5063</td><td>-586.6085</td><td>34.2007</td><td>-2.422</td><td>0.181</td></tr><tr><td>3σ error bar</td><td>0.0017</td><td>0.0024</td><td>0.0067</td><td>0.021</td><td>0.065</td></tr><tr><td>Deviation model measurement</td><td>0.0000</td><td>-0.0001</td><td>0.0001</td><td>0.0000</td><td>0.0027</td></tr><tr><td>Odd J<sub>n</sub> × 10<sup>6</sup></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Wind model (i)</td><td>-0.0569</td><td>-0.0750</td><td>0.1354</td><td>-0.1157</td><td></td></tr><tr><td>Juno measurement</td><td>-0.0450</td><td>-0.0723</td><td>0.120</td><td>-0.113</td><td></td></tr><tr><td>3σ error bar</td><td>0.0033</td><td>0.0042</td><td>0.012</td><td>0.036</td><td></td></tr><tr><td>Deviation between model and measurement</td><td>-0.0119</td><td>0.0027</td><td>0.015</td><td>0.003</td><td></td></tr></tbody></table>

Note. The interior and wind models both make contributions to the even $J_n$, while for the odd $J_n$, only the wind contributions matter. The results of wind model (i) are shown. With wind approach (ii), the Juno measurements can be matched to four significant digits. Machine-readable data files for model A are included in the supplemental material.

zonal gravity coefficients, $J_n$, which are integrals over all mass,

$$J_n = -\frac{2\pi}{Ma^n} \int_{-1}^{+1} d\mu \int_{0}^{r_{\max}(\mu)} dr \ r^{n+2} P_n(\mu) \rho(r, \mu), \quad (1)$$

where $M$ and $a$ are the planet's mass and equatorial radius, $P_n$ are the Legendre polynomials, and $\rho$ represents the planet's density at radius $r$ and colatitude $\theta$ with $\mu = \cos\theta$. Throughout this work, we use $GM = 1.266865341 \times 10^{17} \text{ km}^3 \text{ s}^{-2}$ from Durante et al. (2020), a 1 bar radius of $a = 71492$ km from Lindal et al. (1981), and a rotation period of 9:55:29.711 hr, or $870.536^\circ \text{ day}^{-1}$, from Archinal et al. (2010). The dimensionless rotational parameter then becomes $q_{\text{rot}} = \frac{\Omega^2 a^3}{GM} = 0.08919543238$.

Matching the Juno measurements (Durante et al. 2020) in Table 1 with conventional interior models has been a challenge. Models typically predict values for $J_4$ and $J_6$ that are larger in magnitude than was measured, as illustrated in Figure 1. This discrepancy has made it difficult to draw conclusions from the gravity measurements about Jupiter's interior structure and evolution. Earlier work has demonstrated that it is possible to bridge the $J_4-J_6$ discrepancy by reducing the density of hydrogen, helium, and heavy elements in the planet's interior in a number of ways. One can, e.g., assume a subsolar value of the heavy-element fraction, $Z_1$ (Hubbard & Militzer 2016; Wahl et al. 2017). (For the heavy-element abundance of the protosolar nebula, we assume the value of $Z_{\text{solar}} = 0.0153$ from the composition model of Lodders 2010.) One can also make the interior hotter by raising the interior entropy (Wahl et al. 2017), but the resulting models are no longer compatible with the adiabatic temperature profile starting from $T_{\text{1bar}} = 166.1$ K that the Galileo probe measured (Seiff et al. 1997). Nettelmann et al. (2012) showed that one can also reduce the density in the deep interior by adopting an equation of state (EOS) that has a different entropy than is predicted by ab initio simulations (Militzer & Hubbard 2013). Different EOSs and 1 bar temperatures higher than 166.1 K were recently employed by Miguel et al. (2022) to demonstrate that the heavy-element abundance cannot be constant throughout Jupiter's envelope.

Recently, Debras et al. (2021) proposed that the density in Jupiter's deep interior could be reduced by invoking a stably stratified layer at intermediate pressures. To match the Juno gravity data, such models assume that the temperature profile

of a deep adiabat of composition $Z > Z_{\text{solar}}$ is higher, but its density is lower than an adiabat with $Z \approx Z_{\text{solar}}$. So far, this assumption is not supported by experiments or ab initio simulations (Soubiran & Militzer 2015). Earlier, Debras & Chabrier (2018) addressed the $J_4-J_6$ challenge by reducing the density in an intermediate layer from ~1 to 5 Mbar by adopting a higher entropy (or temperature) and/or a subsolar heavy-element abundance.

In this paper, we match $J_4$ and $J_6$ by adopting a dilute core that extends to ~63% of the planet's radius. It makes the density in the deep envelope higher than in models that assume a compact core. We explain why a dilute core allows us to match $J_4$ and $J_6$. Here "dilute core" refers to a still-evolving state in which the Z component has been dissolved and greatly diluted into a hydrogen-rich envelope. Wahl et al. (2017) showed that approximately one-third of the $J_4-J_6$ discrepancy can be ameliorated by assuming Jupiter has a dilute core. Here we develop this approach further by combining assumptions for the planet's interior with models for winds. We optimize the model parameters for the interior and winds simultaneously, which enables us to match the entire set of even and odd gravity harmonics under one self-consistent set of assumptions.

Our results differ from the predictions in a recent paper on Jupiter's interior by Nettelmann et al. (2021). This paper does not model the effects of winds on zonal harmonics or use our optimization procedure to find models consistent with a given EOS but rather modifies the EOS to obtain a match with the Juno data. Nettelmann et al. (2021) computed zonal harmonics using a seventh-order theory of figures, which brings the predictions into better agreement with concentric Maclaurin spheroid (CMS) results than earlier lower-order calculations. With knowledge of the predictions in this paper, Idini & Stevenson (2022) recently invoked a dilute core to relate large values of the tidal Love number $k_{42}$ to internal waves that are trapped in the core.

## 2. Methods and Model Assumptions

All interior models are derived with the CMS method (Hubbard 2013; Militzer et al. 2019) with 1025 or 2049 spheroids that enables one to model a rotating planet in hydrostatic equilibrium without invoking perturbative methods. For a given pressure and entropy, the density of the mixture of
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Comparison of Juno measurements (diamonds) and model predictions for the even and odd gravity coefficients. The large open symbols show results from the dilute core model A in Table 1, while the small symbols represent ensembles of A- and I-type models. The stars show earlier model predictions from Hubbard & Militzer (2016) and Wahl et al. (2017). The large triangles show models that combine dilute and compact cores. The labels indicate the masses of their compact cores in Earth masses.</figcaption></figure>

hydrogen, helium, and heavier elements is obtained by combining the EOSs of Saumon et al. (1995) at low pressure with the ab initio results (Militzer & Hubbard 2013; Militzer 2013) at high pressure. Heavy elements are incorporated according to Hubbard & Militzer (2016).

In this paper, we report results from models of types A, B, and I that invoke slightly differing assumptions for the interior but all rely on a dilute core and a core transition layer where the heavy-element abundance, Z(P), changes (see Figure 2). In model A, we keep the hydrogen–helium mass ratio constant as Z(P) changes across this layer, $\tilde{Y} \equiv Y/(X + \tilde{Y}) = \text{constant}$. In our initial models of type I, we kept the helium mass fraction Y instead of $\tilde{Y}$ constant across the core transition layer. The

impact on the computed gravity field is small, but we consider keeping $\tilde{Y}$ constant to be more plausible because the hydrogen–helium ratio probably remained constant as heavy elements were added to this layer. The other difference between models A and I is that in models I, we represent the heavy-element abundance by a flexible piecewise linear function of log(P) with 11 knots. As we learned that so much flexibility was not needed, we constructed models of type A for which we assume the heavy-element abundance is constant at low pressure, $Z_1$, as well as at high pressure inside the dilute core, $Z_2$. It changes linearly between these two values across the core transition layer (see Figure 2). In Table 1, we summarize all of the parameters and provide values for a reference model of
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Predictions from models I and A for the mass fractions of heavy elements, Z, and helium, $Y/(X+Y)$, are shown as a function of equatorial radius. Instead of a compact core of $Z \approx 1$, our models include a dilute core that reaches to 63% of the planet's radius ($R = 0.63$). In its inner part, the composition is uniform, with Z values ranging from 0.16 to 0.19. This inner region is surrounded by a stably stratified layer where Z gradually decreases until it reaches a constant, approximately solar value for $R \ge 0.63$. The lower panel shows the helium distribution of our models that are constrained to match measurements of the Galileo entry probe near $R = 1$ and, on average, to agree with the protosolar helium abundance.</figcaption></figure>

type A. For models of type B, we make the same assumptions as for the A models, but we change the EOS for hydrogen-helium mixtures over selected pressure intervals (see Table 3) in order to determine the impact on the inferred heavy-element abundance.

We assume Jupiter's outer molecular layer to be homogeneous and convective. Its entropy, $S_1$, is set to match $T_{1bar} = 166.1$ K. Its helium mass fraction, $Y_1$, is set to match the observed value of $Y/(X+Y) = 0.238$. The only adjustable parameter in this layer is the heavy-element fraction, $Z_1$. We introduce two pressures, $P_{rain,1}$ and $P_{rain,2}$, that mark the boundaries of the helium rain layer. Following Militzer et al. (2019), their values are adjustable but constrained to remain

close to the predictions of the ab initio simulations by Morales et al. (2013) so that the entropy $S_2$ is constrained to be between $S_1$ and a maximum entropy consistent with the helium immiscibility curve. Throughout this helium rain layer, we gradually change the entropy and the helium fraction between the value of the layers above and below.

We set the helium mass fraction of the metallic layer, $Y_2$, so that the planet overall has a protosolar fraction of $Y/(X+Y) = 0.2777$ (Lodders 2010). In the helium rain layer from pressure $P_{rain,1}$ to $P_{rain,2}$, we gradually switch from the exterior adiabat of entropy $S_1$ and helium fraction $\tilde{Y}_1$ to an interior adiabat characterized by $S_2$ and $\tilde{Y}_2$. The mass fraction of heavy elements, $Z_1$, remains constant. We employ the algebraic
switching function,

$$f(P) = x^{\alpha} \quad \text{with} \quad x = \frac{\log(P/P_{\text{rain,1}})}{\log(P_{\text{rain,2}}/P_{\text{rain,1}})}, \qquad (2)$$

to control the entropy, $S(P) = S_1 + f(P) \times (S_2 - S_1)$, and helium fraction, $\tilde{Y}(P) = \tilde{Y}_1 + f(P) \times (\tilde{Y}_2 - \tilde{Y}_1)$, at intermediate pressures. The exponent $\alpha$ is a positive, adjustable parameter. If it is chosen to be larger than 1, more helium has been sequestered from the upper region of the rain layer. Figure 2 shows that models of type A favor this scenario. Conversely, if $\alpha$ were set to a value smaller than 1, less helium would have been sequestered.

We employ two numerical methods to generate models that match the Juno gravity data: the downhill simplex method (Press et al. 2001) and Markov Chain Monte Carlo (MCMC) calculations (Goodman & Weare 2010). The cost function is dominated by the $\chi^2$ deviation between the measured values of the even and odd gravity harmonics, $J_n$, and the model predictions. But then we include additional penalty terms as explained in Militzer et al. (2019) that, for example, ensure that the pressure values for the rain boundary, $P_{\text{rain,1}}$ and $P_{\text{rain,2}}$, are compatible with the H-He immiscibility curve derived by Morales et al. (2013). The downhill simplex method allows one to optimize a single model. The disadvantages of this method are that it tends to get stuck in local minima, and, more importantly, it does not provide a practical way to determine whether a minimum is global or local. This makes it difficult to decide whether model assumptions need to be modified in cases when model predictions do not match the spacecraft measurements well. Furthermore, if there are degeneracies among the best models —if, for example, unnecessary model parameters are included —the downhill simplex will not help to identify them. Conversely, MCMC methods are very efficient in mapping out the allowed parameter region and identifying parameter degeneracies. On the other hand, assessing the quality of the generated ensemble may be a challenge. For example, in the MCMC ensemble of Guillot et al. (2018), there were models that matched the gravity data in $J_4-J_6$ space. There were also models that matched in $J_6-J_8$ and $J_8-J_{10}$ spaces, but there was no model at the time that matched all $J_n$ simultaneously. The simplex method is very good at reoptimizing selected models from the MCMC ensemble to assess their quality. For these reasons, we combine simplex and MCMC methods here, which enables us to provide the reference model in Table 2 for further use in addition to generating model ensembles.

## 2.1. Thermal Wind Models with Latitude-dependent Depth

Guillot et al. (2018) and Kaspi et al. (2018) demonstrated that the winds on Jupiter reach a depth of approximately 3000 km. Here we model their gravity effects by solving the thermal wind equation (Kaspi et al. 2016) for a rapidly rotating planet. We adopt the time-averaged wind profiles that were observed by tracking cloud motion (Tollefson et al. 2017) and assume they remain initially constant in the direction parallel to the planet spin axis but then decay at some depth $H$. Then we construct two ensembles of wind models under slightly different assumptions. Both fit the Juno data when combined with our interior models. In this section, we describe our wind approach (i) that directly employs the observed cloud-level wind profile (Tollefson et al. 2017) without modifications.

Matching the even harmonics $\Delta J_n$ requires us to make the wind depth latitude-dependent, as shown in Figure 3. In Section 2.2, we introduce our wind approach (ii) that assumes the wind depth and decay function to be independent of latitude but instead allows the wind profiles to deviate from the cloud-level observations (Figure 3).

In addition to providing the $J_n$ in Equation (1), the CMS method also gives access to the surfaces of constant potential (gravity plus centrifugal terms) throughout the planet's interior. We use these surfaces and the density structure to solve the thermal wind equation (Kaspi et al. 2016) for a rotating, oblate planet (Cao & Stevenson 2017a) in geostrophic balance. On an equipotential surface, we construct paths, $s$, from the equatorial plane to the poles and integrate the dynamical part of the density, $\rho'(s)$, using

$$\frac{\partial \rho'}{\partial s} = \frac{2\Omega}{g} \frac{\partial}{\partial z} [\rho u], \qquad (3)$$

where $z$ is the vertical coordinate that is parallel to the axis of rotation; $\rho$ is the static background density that we derive with the CMS method; $u$ is the differential flow velocity with respect to the uniform rotation rate, $\Omega$, that one obtains from the planet's rotation period; and $g$ is the acceleration that we derive from the gravitational-centrifugal potential in our CMS calculations. We represent the flow field $u$ as a product of the surface winds, $u_s$, and a decay function, $D$,

$$u = u_s(\theta) \times D(d, H(\theta)), \qquad (4)$$

where $d$ presents the distance from the surface, and $H$ represents the wind depth that we allow to vary between 1000 and 5000 km with colatitude $\theta$. We assume that $u$ remains initially constant but then decays to zero over a depth interval from $d_1 = H[1 - w/2]$ to $d_2 = H[1 + w/2]$ according to

$$D(d, H) = D\left(x = \frac{d - d_1}{d_2 - d_1}\right) = \begin{cases} 1 & \text{for } x \le 0 \\ \sin^2\left((1 - x)\frac{\pi}{2}\right) & \text{else} \\ 0 & \text{for } x \ge 1 \end{cases}, \qquad (5)$$

where $w$ is the width of the decay interval. We obtain good models by setting $w$ to 0.1 and 0.2.

## 2.2. Thermal Wind Model with Modified Cloud-level Wind Profiles

Following Kaspi et al. (2020) and Galanti & Kaspi (2021), we take the approach of allowing the cloud-level wind to be modified from the observed values. Unlike in methodology (i), in which the wind depth varies with latitude, for each solution here, we find an optimal wind depth and decay profile that we apply to all latitudes. For an ensemble of interior models of type A, we fit exactly the odd gravity harmonics $J_3$, $J_5$, $J_7$, and $J_9$ and the residual even harmonics $\Delta J_2$, $\Delta J_4$, $\Delta J_6$, $\Delta J_8$, and $\Delta J_{10}$. We start by decomposing the observed cloud-level wind profile into N Legendre polynomials, $P_i$,

$$U^{\text{obs}}(\theta) = \sum_{i=0}^{N} A_i^{\text{obs}} P_i(\cos \theta). \qquad (6)$$
<table><caption>**Table 2**

Description of Model Parameters and Values in Our Preferred Model A</caption><thead><tr><th>Parameter</th><th>Value in Model A</th><th>Description</th></tr></thead><tbody><tr><td>S<sub>1</sub> (k<sub>B</sub>/el.)</td><td>7.078</td><td>Entropy in the outer molecular layer. Fixed to match the temperature at 1 bar of 166.1 K in all models.</td></tr><tr><td>Z<sub>1</sub></td><td>0.0156</td><td>Adjustable parameter that represents the mass fraction of heavy elements in the molecular layer. A penalty is added to models with Z<sub>1</sub> &lt; Z<sub>solar</sub>.</td></tr><tr><td>Y<sub>1</sub></td><td>0.2332</td><td>Helium mass fraction in outer layer. Fixed to match the observed value of Y<sub>1</sub>/(1 - Z<sub>1</sub>) = 0.238.</td></tr><tr><td>P<sub>rain,1</sub> (GPa)</td><td>93.1</td><td>Adjustable parameter for the starting pressure of the helium rain layer. The corresponding temperature, T<sub>1</sub>, follows from S<sub>1</sub> and the EOS. As explained in Militzer et al. (2019), a penalty is introduced if P<sub>rain,1</sub> and T<sub>1</sub> deviate from the H-He immiscibility curve of Morales et al. (2013).</td></tr><tr><td>P<sub>rain,2</sub> (GPa)</td><td>443.2</td><td>Adjustable parameter to represent the high-pressure end of the helium rain layer. T<sub>2</sub> follows from S<sub>2</sub>. Again, a penalty is introduced if P<sub>rain,2</sub> and T<sub>2</sub> deviate from the immiscibility curve.</td></tr><tr><td>α</td><td>9.4</td><td>Exponent in helium rain switching function, also an adjustable parameter.</td></tr><tr><td>S<sub>2</sub> (k<sub>B</sub>/el.)</td><td>7.194</td><td>Adjustable parameter that represents the entropy in the metallic layer. It cannot exceed 7.20 k<sub>B</sub> electron<sup>-1</sup> to be compatible with the H-He immiscibility curve.</td></tr><tr><td>Y<sub>2</sub></td><td>...</td><td>Helium mass fraction below the helium rain layer that is constrained so that the planet overall matches the solar value of Y/(X + Y) = 0.2777 from Lodders (2010). Here Y<sub>2</sub> is only employed in models of type I, where Y<sub>2</sub>(P) is kept constant as Z(P) varies in the core transition layer.</td></tr><tr><td>Ỹ<sub>2</sub> = Y<sub>2</sub>/(1 - Z<sub>2</sub>)</td><td>0.2957</td><td>Employed in models A and B to keep the helium fraction Ỹ<sub>2</sub>(P) below the helium rain layer constant as Z(P) varies. Its value is contained so that the planet overall matches the solar helium fraction.</td></tr><tr><td>P<sub>core,1</sub> (GPa)</td><td>784</td><td>Adjustable parameter in models A and B that marks the outer boundary in pressure of the dilute core. It affects the heavy elements according to Z(P &le; P<sub>core,1</sub>) = Z<sub>1</sub>.</td></tr><tr><td>P<sub>core,2</sub> (GPa)</td><td>2054</td><td>Adjustable parameter in models A and B that marks the inner boundary in pressure where the Z of the core starts to decrease.</td></tr><tr><td>Z<sub>2</sub>(P)</td><td>...</td><td>Piecewise linear function that represents a heavy-element fraction in the metallic layer. It is only used in models of type I. The Z<sub>2</sub> value is adjusted on 11 pressure points.</td></tr><tr><td>Z<sub>2</sub></td><td>0.1830</td><td>Adjustable parameter in models A and B. We set Z(P &ge; P<sub>core,2</sub>) = Z<sub>2</sub>.</td></tr><tr><td>w</td><td>0.1</td><td>Fractional width of the wind decay interval, typically kept fixed at 0.1 or 0.2, but models that vary w with latitude have also been constructed.</td></tr></tbody></table>

**Note.** Machine-readable data files for model A are included in the supplemental material.

(This table is available in its entirety in machine-readable form.)

The resulting coefficients, $A_i^{\text{obs}}$, represent the latitudinal shape of the observed wind profile from Tollefson et al. (2017). Then we introduce a second set of coefficients, $A_i^{\text{sol}}$, to represent a modified cloud-level wind profile,

$$U^{\text{sol}}(\theta) = \sum_{i=0}^{N} A_i^{\text{sol}} P_i(\cos \theta), \quad (7)$$

to represent solutions that may deviate from the observations. The coefficients $A_i^{\text{sol}}$ are optimized for the wind-induced gravity harmonics to fit the spacecraft observations perfectly. We employ a very large number of polynomials, $N=99$, so that the emerging wind solution follows the observed wind profile as close as possible. The following optimization procedure ensures that the large set of coefficients are well constrained. We perform these calculations in a spherical geometry but otherwise follow the same steps as in method (i) when projecting the cloud-level winds inward, allowing it to decay in the radial direction, calculating the induced anomalous

density field, and calculating the wind-induced gravity harmonics.

Our optimization procedure is based on the methodology of Kaspi & Galanti (2016) and Galanti & Kaspi (2021). The parameters to be optimized, i.e., those defining the depth of the wind and the cloud-level wind latitudinal profile, are defined as a control vector,

$$Y = \left\{ \frac{H_0}{h_{\text{norm}}}, X^{\text{sol}} \right\} \quad \text{with} \\ X^{\text{sol}} = \left\{ \frac{A_0^{\text{sol}}}{u_{\text{norm}}}, \dots, \frac{A_N^{\text{sol}}}{u_{\text{norm}}} \right\}, \quad (8)$$

where the parameter $H_0$ represents the radial wind depth, while $h_{\text{norm}} = 10^7$ m and $u_{\text{norm}} = 10^3$ m s$^{-1}$ are the normalization factors for the depth and coefficients, respectively. The goal is to minimize the difference between the model solution and the gravity observation and keep the parameters from attaining unphysical
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Wind speeds as a function of latitude observed in Jupiter's atmosphere (Tollefson et al. 2017). We assume they decay with depth and incorporate them into our thermal wind calculations. In our thermal wind approach (i), this wind profile is adopted without modification, but the wind depth depends on latitude. The gray profiles in the lower panel show an ensemble of wind depth profiles for I-type interior models. The red curve shows the wind depth profiles for our preferred interior model A. In our thermal wind approach (ii), the reconstructed cloud-level wind profiles (yellow region in upper panel) are permitted to deviate from the observations, while the wind decay function is latitude-independent. The dashed orange curve shows the wind solution for the preferred interior model A.</figcaption></figure>

values. The cost function, $L$, is composed of two terms:

$$L = (J^{\text{sol}} - J^{\text{obs}})W(J^{\text{sol}} - J^{\text{obs}})^T + \epsilon_U(X^{\text{sol}} - X^{\text{obs}})(X^{\text{sol}} - X^{\text{obs}})^T. \quad (9)$$

The first is the difference between the measured and calculated gravity harmonics, and the second ensures that the wind solution does not vary too far from the observed one. Here $J^{\text{sol}}$ is a vector that contains the model solution for $J_3, J_5, J_7, J_9, \Delta J_6, \Delta J_8$, and $\Delta J_{10}$, while the vector $J^{\text{obs}}$ represents the corresponding measured values. The covariance matrix, $W$, represents the uncertainties of the gravity measurements (Durante et al. 2020). We set $\epsilon_U = 5 \times 10^8$. Here $X^{\text{obs}}$ are the normalized observed wind coefficients. Given the value of $\epsilon_U$ and the large number of coefficients, $A_i^{\text{sol}}$, the wind is strongly constrained to the observed cloud-level profile, thus ensuring that deviations from the observed values are only permitted if they result in a significant lowering of the cost function. Given an initial guess for $Y$, a minimal value of $L$ is derived using the Matlab function `fmincon`<fnref n="19" /> while taking

advantage of the cost function gradient that is derived from the adjoint of the dynamical model.

### 3. Results and Discussion

Here we show that all of the zonal gravity harmonics can be matched with one set of plausible assumptions, including a dilute core, a helium rain layer, and a model for the zonal wind speeds and depth. Motivated by the work of Wahl et al. (2017), we assume a flexible profile for the heavy-element abundance in the deep interior $Z(P)$ that we represent by a piecewise linear function of $\log(P)$ (see models I in Figure 2). When we simultaneously optimized our interior and wind parameters under these assumptions, we found that all promising models had no or only a very small compact core. Instead, the heavy elements were distributed throughout the deep interior, extending to 63% of the planet's radius, as illustrated in Figure 2. We thus temporarily removed the compact central core from our models and extended the metallic layer to the planet's center. In Section 3.2, we explain why redistributing mass from the compact core reduces the magnitude of $J_4$.

The most promising models of type I that then emerged had a number of features in common. In the dilute core region, the abundance of heavy elements was constant. Above some radius

<footnote n="19">¹⁹ MathWorks: https://www.mathworks.com/help/optim/ug/fmincon.html.</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Temperature and density profiles of model A. The vertical dotted lines mark the boundaries of the helium rain region. The dashed lines show the boundaries of the transition layer. The bottom panel shows the density deviation from a protosolar adiabat ($Y = 0.2735$ and $Z = 0.0153$; Lodders 2010) for $T_{1\text{bar}} = 166.1$ K. Because of the helium sequestration, the density of model A is lower in the molecular and helium rain layers. Conversely, in the dilute core and the core transition layer, the density of model A is much higher because of the enrichment in heavy elements. In the metallic hydrogen layer, the density of model A is slightly higher than that of the protosolar adiabat because of the enrichment in helium.</figcaption></figure>

(quantified below), it started to decay gradually until it reached the value of the outer layer, $Z_1$. We consequently simplified how we represent the distribution $Z_2(P)$. We kept only two values, $Z_1$ and $Z_2$, but introduced two pressure values, $P_{\text{core},1}$ and $P_{\text{core},2}$, that mark the region of decay from $Z_2$ to $Z_1$ with decreasing pressure, as we illustrate in Figure 2. We require this region and the helium rain layer to be Ledoux stable (Ledoux 1947). Under these assumptions, we obtained model A (see Table 2 and Figure 4), whose even gravity coefficients match the Juno observations exceptionally well (Table 1). We still see some minor deviations, e.g., for the odd gravity harmonics $J_3$ but the magnitude of the deviations is reasonably small so that they can be eliminated with our second set of wind approach.

Summarizing, one can say that the dilute core with a constant $Z$ distribution emerged from our flexible I-type models because we tried to match as closely as possible the measured $J_4$ and $J_6$ values but also because we assume that the core region has a

constant entropy and require that the $Z$ fraction does not decrease as one descends into the planet (Ledoux stability). When we constructed our subsequent models of type A, we assumed $Z$ to be constant inside the dilute core.

For 8572 interior models of type A that we constructed with MCMC calculations, we derived wind solutions with approach (ii) to match the Juno measurements. The resulting wind models are shown in Figure 5. Since the observed wind (Tollefson et al. 2017) already allows a reasonable match to the higher gravity harmonics (Kaspi et al. 2020), only moderate modifications are needed to enable an exact fit to all of the even gravity harmonics (Figure 5(a)). At midlatitudes, the deviation of the optimized wind is mostly within the observed uncertainty of around $15 \text{ m s}^{-1}$ (Tollefson et al. 2017), and in the equatorial region, the deviations are larger, reaching $50 \text{ m s}^{-1}$. With the modified cloud-level winds, a perfect match to all gravity harmonics is reached (Figures 5(b)–(e)). For all solutions, the depth of the wind is around 3000 km.
<figure><img src="image_9.png" /><figcaption>**Figure 5.** (a) Observed cloud-level wind (black) and the range of model solutions (red shaded area) developed with thermal wind approach (ii). (b)–(e) Corresponding model solutions for the even gravity harmonics. Shown are the measurements (black crosses, with the size of the cross representing the measurement uncertainty), the interior model solutions (blue dots), and the model solutions when the wind-induced gravity harmonics are added to the interior model solutions (red dots). The insets show the zoom-in on the measured values.</figcaption></figure>

In Figure 1, we separate the contributions from our interior and wind models. Diluting the core changes the interior contributions to $J_4$ and $J_6$ linearly until the residual discrepancy can be bridged with a wind model that we derive self-consistently with wind approach (i). The green and yellow symbols respectively represent the interior and interior+wind contributions to the different gravity harmonics that we derived from an MCMC ensemble of A-type models. The wind models contribute $\Delta J_4 \approx 0.24 \times 10^{-6}$ and $\Delta J_6 \approx -0.27 \times 10^{-6}$, as well as smaller amounts to all other $J_n$. The magnitude of the wind corrections is broadly consistent with that reported by Guillot et al. (2018). Still, the characterization of the wind structure below the visible cloud desk remains a goal of the Juno extended mission. From the single-spot wind measurements of the Galileo entry probe (Atkinson et al. 1998), one

may expect deviations between the deep flows and those visible at the cloud level.

We also inserted small compact cores composed of a 1:1 rock-ice mixture into our dilute core models and reoptimized all model parameters. We were able to accommodate compact cores of up to 3 $M_⊕$ (1% of Jupiter's mass). For larger compact cores, we cannot fit the gravity data. Already for 4 $M_⊕$, we find a discrepancy of $\Delta J_4 \approx 0.9 \times 10^{-6}$ (or $\sim 10^3\sigma$; Durante et al. 2020) between predictions from interior models and gravity measurements that we cannot bridge with our wind model assumptions. Figure 1 shows that for larger compact cores of up to 8 $M_⊕$, the gap between model predictions and the Juno measurements widens in $J_4-J_6$ space, as mass is being moved from the dilute core region into the compact core.

According to model A, Jupiter has a thick central region extending to $0.41 R_J$ that is fully convective and uniformly
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Inner and outer radii of the core transition layer derived from ensembles of type A models (left panel). The probability densities in the two other panels are normalized so that the integrals under the curves equal 1. For the outer and inner radii, we derived values of 0.63 ± 0.01 and 0.41 ± 0.03, respectively. The uncertainty intervals were chosen so that approximately 80% of all models are represented.</figcaption></figure>

enriched by ~12 $M_⊕$ worth of heavy elements (mass fraction 18%; see Figure 2). Figure 4 shows that the density of model A is higher than that of the protosolar adiabat except in the outer two layers, where the depletion of helium reduced it by ~3%. Note that our dilute core solution differs from other recent Jupiter models with dilute cores that invoke compositional gradients and are stable against convection, and heat is transported via conduction/radiation or semiconvection (Vazan et al. 2018; Debras & Chabrier 2019). This may imply that Jupiter's primordial core was mixed with the envelope during a giant impact shortly after its formation (Liu et al. 2019). Analyses of the possible evolution of a primordial core show that (a) likely constituents of the core are fully soluble in hydrogen under the relevant temperature-pressure conditions (Wilson & Militzer 2012a, 2012b; Wahl et al. 2013; Gonzalez-Cataldo et al. 2014), and (b) the primordial core could be eroded by convective motions (Moll et al. 2017; Müller et al. 2020). The convective central region is expected to be the primary source of Jupiter's magnetic field.

The core region is enveloped in a stably stratified transition layer (0.41–0.63 $R_⊕$) where the heavy-element abundance gradually decreases, contributing ~11 $M_⊕$ of heavy elements. This layer is neither convective nor likely to be capable of generating a dynamo. Such layers might be analogous to those invoked for Saturn to explain seismological observations in its rings (Fuller 2014) and the planet's luminosity anomaly (Leconte & Chabrier 2013).

Atop the transition layer in Jupiter is a thin layer of metallic hydrogen (0.63–0.72 $R_⊕$) that is again homogeneous and convective (inset of Figure 2), raising the possibility of a secondary dynamo operating there. Above that is the helium rain region (0.72–0.83 $R_⊕$) that would again be stably stratified. Finally, the outer molecular layer is again homogeneous and convective. It only contributes ~2 $M_⊕$ of heavy elements, bringing the total amount in the planet to ~25 $M_⊕$, which is a lower bound, as we discuss in Section 3.3.

In Figure 6, we show the inner and outer equatorial radii of the core transition region that derived from an ensemble of models that are constructed with the same assumption as model A. Both radii are anticorrelated so that the mass of the core region is approximately preserved when these radii vary. If the outer boundary of this region is extended to a larger radius, the inner boundary shrinks so that the heavy Z-element abundance falls off more gradually over a larger radius interval. In Figure 6,

<table><caption>**Table 3**
Model Predictions for the Heavy-element Fraction in the Molecular Envelope, $Z_1$</caption><thead><tr><th>Model Type</th><th>Pressure Interval [GPa]</th><th>Density Correction Δρ/ρ</th><th>Z<sub>1</sub></th><th>Z<sub>1</sub>/Z<sub>solar</sub></th></tr></thead><tbody><tr><td>A</td><td>None</td><td>0</td><td>0.0156</td><td>1.00</td></tr><tr><td>B1</td><td>20–100</td><td>-0.02</td><td>0.021</td><td>1.37</td></tr><tr><td>B2</td><td>20–4500</td><td>-0.02</td><td>0.034</td><td>2.25</td></tr><tr><td>B3</td><td>10–100</td><td>-0.03</td><td>0.033</td><td>2.14</td></tr><tr><td>B4</td><td>10–4500</td><td>-0.03</td><td>0.045</td><td>2.92</td></tr></tbody></table>

**Note.** In models B1–B4, the density of our H–He EOS has been lowered over the specific pressure range by a certain fraction compared to our original model A.

we also show the probability distributions of the outer and inner radii in the model ensemble. From these distributions, we derived values of 0.63 ± 0.01 and 0.41 ± 0.03, $R_⊕$ for both radii. The uncertainties reflect approximately 80% of the models in the ensemble.

### 3.1. Dependence on the EOS of Hydrogen–Helium Mixtures

Our models of types A and I predict an approximately protosolar enrichment of $Z_1 = 0.0156$ for the molecular layer, apparently consistent with the lower end of the equatorial water determination at 0.7–30 bars of Li et al. (2020) between one and five times solar, only if the Galileo determination of the carbon enrichment of two to four times solar (Wong et al. 2004) were not included. Assuming that these two low-pressure measurements are representative of the entire molecular envelope, the carbon enrichment and $Z_1 = 0.0156$ would force the abundance of oxygen, and hence water, to be subsolar. Since the Juno spacecraft has so far determined the water abundance only near the equator, measurements at other latitudes are warranted.

Here we probe how sensitively the predicted heavy-element abundance for the molecular layer, $Z_1$, depends on the hydrogen–helium EOS that was derived with ab initio simulations. Models B1–B4 in Table 3 were derived by lowering the density by a certain fraction over the specific range of pressures.

While there are no uncertainties in the EOS at low pressure, where the H–He mixture behaves like an ideal gas, it is
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Model predictions for the amount of heavy elements in the molecular envelope, $Z_1$, if we introduce a density correction, $\Delta\rho/\rho$, for all pressures $P \ge P_c$ to the EOS in our A-type models and then reoptimize them to match the Juno measurements. The blue curves represent the results from simple reoptimizations that only match the planet's mass and $J_2$. The red curves were obtained by matching the mass and all gravity coefficients, which led us to conclude that a density correction only in the metallic region for $P \ge 200$ GPa does not significantly affect $Z_1$. Instead, $Z_1$ was found to be sensitive to a density reduction in the pressure interval from 10 to 200 GPa.</figcaption></figure>

reasonable to ask what type of EOS changes at high pressure would lead to an increase of the $Z_1$ values that are predicted by our models. There are uncertainties in the ab initio EOS (Militzer & Hubbard 2013) that were derived with the density functional theory. While this is a state-of-the-art approach, an experimental confirmation is still outstanding because the uncertainties of the existing density and temperature measurement data at megabar pressures are still too large to constrain the results from ab initio simulations.

For any pressure and temperature, an increase in $Z$ always increases the density. We therefore studied how our model predictions would change if the density of our H-He EOS were decreased by a fraction of 1%, 2%, or 3% for all pressures above a certain threshold, $P_c$. We performed this analysis in two steps and report the results in Figure 7. In the first step, we introduce the density change into model A and only match the planet and $J_2$ by adjusting $Z_1$ and $Z_2$. The blue curves in Figure 7 show a substantial increase in $Z_1$ even if the density is only changed for a high value of $P_c \sim 200$ GPa.

In the second step, we reoptimize the entire model with the goal of matching the planet mass and the entire set of gravity coefficients. The red curves in Figure 7 show that a density change $P_c \ge 100$ GPa does not lead to an increase in $Z_1$ because in our models, the helium rain layer starts at $\sim 100$ GPa; thus, our models are less sensitive to a density change in this layer. From these calculations, we selected models B2 and B4, which respectively incorporated a density reduction of 2% for $P_c \ge 20$ GPa and 3% for $P_c \ge 10$ GPa. The resulting $Z_1$ values are reported in Table 3.

We concluded that our models are fairly sensitive to density change in the pressure from 10 to 100 GPa. We therefore constructed two additional models, B1 and B3, in which we reduced the density by 2% and 3% only in the pressure ranges

from 20 to 100 and 10 to 100 GPa, respectively. Under these assumptions, the heavy-element abundance in the molecular layer, $Z_1$, increases to 1.41 and 2.14 times the protosolar value, which is in better agreement with the noble gas and methane abundances in Jupiter's atmosphere. To remove this uncertainty in the EOS, one would need, e.g., to conduct new laboratory experiments that measure the density of H-He mixtures to better than 2% in the pressure range 10-100 GPa while carefully controlling the temperature.

## 3.2. Relation of Dilute Core and $J_4$

There is an approximate way to understand the indirect evidence that Jupiter has an extended, low density contrast central region of enhanced metallicity but is still dominated by hydrogen-helium, rather than a compact high-density core, as had been expected prior to the Juno mission. Consider a traditional pre-Juno Jupiter model from the set published in Hubbard & Militzer (2016). For such a reference model, we select the model DFT-MD 7.13 (bold in Table 1 in Hubbard & Militzer 2016). For the present purpose, this model differs negligibly from the Juno-measured values of $J_2$ and $M$ and the adopted equatorial radius $a$ and rotational parameter $q_{rot}$, but its predicted $J_4 = -596.05 \times 10^{-6}$ is 1.6% larger in absolute value than the Juno measurement. We will refer to this 2016 Jupiter model as the preliminary Jupiter model (PJM). The PJM has a high-density central core of 12.2 $M_E$ with a radius of 0.15$a$.

Here we report a Jupiter interior model that fits all of the Juno gravitational constraints, including $J_4$. Because Jupiter's interior structure approximately resembles a polytrope of index 1, we can use a perturbation approach to examine the effect of a small density modification on the predicted $J_4$. This is possible because unlike $J_2$, the zonal harmonic $J_4$ and all higher zonal harmonics (in the absence of wind corrections) are excited by
<figure><img src="image_12.png" /><figcaption>**Figure 8.** The density modification and corresponding correction to $J_4$ is shown as a function of equatorial spheroid radius $\lambda_i$.</figcaption></figure>

higher-order responses to $q_{\text{rot}}$. Accordingly, Hubbard (1974) derived, for a polytrope-like planet, the approximate relation

$$J_4 \approx -\frac{4\pi}{35} \left[ J_2 + \frac{1}{3} q_{\text{rot}} \right] \frac{b^3}{M} \rho(b), \quad (10)$$

where $b$ is the planet's polar radius, and the density $\rho(b)$ is evaluated on a level surface whose mean radius is $b$, corresponding to a CMS value of about $\lambda \approx 0.93$ (corresponding to a pressure of about 20 GPa in Jupiter's interior). As argued by Hubbard, the density distribution deeper in the polytrope-like planet contributes little to the total external value of $J_4$.

Define $\delta J_4 = J_4 - J_4^{\text{PJM}}$, where $J_4$ is the value for the present optimized model, fitted to the Juno measurement, and $J_4^{\text{PJM}} = -596.05 \times 10^{-6}$. Similarly, let $\delta\rho = \rho - \rho^{\text{PJM}}$, evaluated at $\lambda \approx 0.93$. Using Equation (10) and holding $J_2$, $M$, $q_{\text{rot}}$, and $b$ constant, we get

$$\frac{\delta J_4}{J_4} \approx \frac{\delta \rho_4}{\rho_4}, \quad (11)$$

where $\rho$ and $\delta\rho$ are evaluated around $\lambda \approx 0.93$. Recognizing that Equations (10) and (11) are crude approximations, they correctly show the order of magnitude of the decrease of $|J_4|$ due to negative $\delta\rho = \rho - \rho^{\text{PJM}}$. Figure 8 shows a detailed comparison of the present optimized Jupiter model with the PJM.

A literal application of Equation (11) to the layer at the equatorial radius of $\lambda \approx 0.93$ would predict $\delta\rho/\rho \approx -0.016$, about three times larger than the gray curve in this region, but the sign and order of magnitude are correct. The bulk of the downward shift in $|J_4|$ arises from spheroids between $\lambda \approx 0.75$ and 0.96; a smaller upward contribution comes from spheroids between $\lambda \approx 0.75$ and 0.6. There is no significant contribution to $|J_4|$ from deeper layers.

As we see from Figure 8, $J_4$ is fitted by shifting the density downward by 1%–4% in the outer layers of Jupiter. Once we have done this, we must readjust the layers for $\lambda < 0.6$ in order

to maintain constant the values of $J_2$ and $M$. With the possible exception of extreme, unrealistic interior models, there is a correlation between $J_2$ and the dimensionless moment of inertia, $C/\text{Ma}^2$. Thus, the value of $C/\text{Ma}^2$ must be kept constant as well, always very close to 0.264. If $\delta\rho/\rho$ is negative for $\lambda > 0.65$, it must be positive for smaller values of $\lambda$. The value of $M$ can be kept constant by inserting a small, high-density central core, but the additional requirement that $C/\text{Ma}^2$ be constant requires the density augmentation to be spread out to a larger central region. The effect of the starting pressure $P_c$ for the density correction $\delta\rho/\rho$, discussed above, can also be understood in the framework of Equation (11), since $\lambda \approx 0.93$ corresponds to about 20 GPa pressure, and $\lambda \approx 0.86$ corresponds to about 100 GPa.

The link between a dilute core and the increase in $J_4$ is also explained by models that keep Jupiter's equatorial radius and $J_2$ value fixed while the density profile is arbitrarily increased over a range of pressures. Figure 5 from Guillot (1999) showed that a 5% increase in density in the 10 Mbar region yields an increase of $J_4$ by about $10^{-5}$, which is in line with the trend that we show in Figure 1.

### 3.3. Superadiabatic Temperature Gradients in Core Transition Layer

In our models, the dilute core is surrounded by a stably stratified transition layer in which the heavy-element abundance increases from $Z \approx 0.0156$ to 0.183. While we did not do so in model A, one might adopt a superadiabatic temperature gradient for this layer (Leconte & Chabrier 2013). This would increase the temperature in the dilute core and thereby enable it to accommodate more heavy elements. We study such an increase under the following assumptions and report the resulting increase in the total Z-element budget of the planet in Figure 9.

While Debras & Chabrier (2019) assumed the dilute core itself to be stably stratified, we assume it to be homogeneous and convective and confine the superadiabatic temperature gradient to the core transition layer. The dilute core may then have an entropy, $S_3$, that is higher than that of the metallic layer, $S_2$. In the core transition layer, we interpolate the entropy linearly from $S_2$ to $S_3$ as a function of $\log(P)$ from $P_{\text{core},1}$ to $P_{\text{core},2}$,

$$S_{\text{tr}}(P) = S_2 + g(P) \times (S_3 - S_2) \\ \geq S_2 \quad \text{with} \quad g(P) = \frac{\log(P/P_{\text{core},1})}{\log(P_{\text{core},2}/P_{\text{core},1})}. \quad (12)$$

We employ the hydrogen–helium EOS (Militzer & Hubbard 2013) to relate the entropy to temperature before we adjust the heavy-element fraction, $Z_{\text{tr}}(P)$, so that the density for a given pressure stays exactly the same as in model A,

$$\rho(S_{\text{tr}}(P), Z_{\text{tr}}(P)) = \rho(S_2, Z_A(P)) \quad \text{with} \\ Z_{\text{tr}}(P) \geq Z_A(P) + g(P) \times (Z_2 - Z_1). \quad (13)$$

Here entropy and the heavy Z fraction trade off against each other while leaving the density unchanged. This means that the gravity coefficients are the same as in model A and that the superadiabatic model is again Ledoux stable (Ledoux 1947). In Figure 9, we plot the revised total budget of heavy elements. For simplicity, we chose the central temperature rather than the core entropy, $S_3$, as the X-axis. We find that the planet's total
<figure><img src="image_13.png" /><figcaption>**Figure 9.** If one assumes a superadiabatic temperature gradient (circles) in the stably stratified transition layer (0.41–0.63 *R*₁), one can construct interior models that accommodate a larger amount of heavy elements (Y-axis) than in model A (diamond) because the temperature in the core region is higher. We represent this temperature change by plotting against the temperature in the planet's center (X-axis).</figcaption></figure>

heavy-element contents increase linearly with the planet's central temperature. If the central temperature doubles in magnitude, the total Z amount increases from 25.0 to 41.7 *M*⊕. In this regard, our adiabatic model A represents a lower bound for the amount of heavy elements stored in the planet. At present, we are not able to place an upper bound on the core entropy because that depends on the planet's formation and how efficiently heat can be transported across the stably stratified layer. The slope in Figure 9 depends on where in the planet the superadiabatic region is introduced and what mass fraction is affected by the resulting temperature increase. While we chose to make the core transition layer superadiabatic, one could place it higher in the planet or put it inside the dilute core region, thus introducing a larger or smaller increase in the amount of heavy elements than we show in Figure 9.

## 4. Conclusions

The Juno mission's measurements of Jupiter's gravity require a thorough revision of models of the interior of the solar system's largest planet. Doing so has proved to be a challenge given other constraints that must be satisfied, including atmospheric temperature and abundances of elements heavier than helium. Here we report models of Jupiter's interior that meet these constraints. Our models have a dilute core—a region substantially enriched in heavier elements—spanning 63% of the planet's radius. This unexpected feature challenges conventional models of the formation and early evolution of Jupiter.

Furthermore, our models have important implications for the dynamo process inside Jupiter. We predict the existence of two separate, fully convective dynamo layers inside Jupiter, the metallic layer and the dilute core, which are separated by a stably stratified layer, as shown in Figure 2. Such a double-dynamo configuration could help explain the planet's

surprisingly complex magnetic field (Connerney et al. 2018) as proposed in Moore et al. (2018). Wind activity in the bottom part of the convective molecular hydrogen layer just above the helium rain layer could further modify the observed magnetic field (Cao & Stevenson 2017b) but is unlikely to determine the main features of the magnetic field. Our interior model implies the simultaneous operation of a deep-seated, thick-shell dynamo and a shallower, thin-shell dynamo, in contrast to the homogeneous dynamo region assumed in most Jovian dynamo studies (Jones 2014; Duarte et al. 2018). The upcoming low-altitude polar passes of the Juno Extended Mission will map the field in the northern hemisphere with increased spatial resolution and test this hypothesis.

Finally, our results beg the question of how a dilute core can form inside Jupiter. Standard models of the formation of Jupiter, whether by core accretion or disk instability, do not predict such a structure (Müller et al. 2020). An obvious possibility is that a discrete primordial core was mixed with the envelope, either over a cosmic period of time by subsequent double-diffusive mixing and solution in the metallic hydrogen (Gonzalez-Cataldo et al. 2014; Moll et al. 2017) or abruptly during a giant impact shortly after its formation (Liu et al. 2019). What such a catastrophic model implies for the other aspects of the Jovian system, such as the presence of four large regular satellites, are intriguing questions stimulated by Jupiter's surprising interior. Among the ways to verify the predictions in this paper are high-pressure experiments on hydrogen, additional measurements by the Juno spacecraft, and Jovian seismology.

## ORCID iDs

Burkhard Militzer https://orcid.org/0000-0002-7092-5629

William B. Hubbard https://orcid.org/0000-0003-3185-4538

Sean Wahl https://orcid.org/0000-0003-2451-7939

Jonathan I. Lunine https://orcid.org/0000-0003-2279-4131
Eli Galanti https://orcid.org/0000-0002-5440-8779

Tristan Guillot https://orcid.org/0000-0002-7188-8428

Kimberly M. Moore https://orcid.org/0000-0003-3162-437X

Marzia Parisi https://orcid.org/0000-0003-4064-6634

Ravid Helled https://orcid.org/0000-0001-5555-2652

Hao Cao https://orcid.org/0000-0002-6917-8363

Christopher Mankovich https://orcid.org/0000-0002-4940-9929

David J. Stevenson https://orcid.org/0000-0001-9432-7159

Ryan S. Park https://orcid.org/0000-0001-9896-4585

## References

Archinal, B., A'Hearn, M., Bowell, E., et al. 2010, CeMDA, 109, 101

Atkinson, D. H., Pollack, J. B., & Seiff, A. 1998, JGR, 103, 22911

Cao, H., & Stevenson, D. J. 2017a, JGRE, 122, 686

Cao, H., & Stevenson, D. J. 2017b, Icar, 296, 59

Connerney, J. E. P., Kotsiaros, S., Oliversen, R. J., et al. 2018, GeoRL, 45, 2590

Debras, F., & Chabrier, G. 2018, A&A, 609, A97

Debras, F., & Chabrier, G. 2019, ApJ, 872, 100

Debras, J., Chabrier, G., & Stevenson, D. J. 2021, ApJL, 913, 21

Duarte, L. D. V., Wicht, J., & Gastine, T. 2018, Icar, 299, 206

Durante, D., Buccino, D. R., Tommei, G., et al. 2020, GeoRL, 47, e86572

Fuller, J. 2014, Icar, 242, 283

Galanti, E., & Kaspi, Y. 2021, MNRAS, 501, 2352

Gonzalez-Cataldo, F., Wilson, H. F., & Militzer, B. 2014, ApJ, 787, 79

Goodman, J., & Weare, J. 2010, Communications in Applied Mathematics and Computational Science, 5, 65

Guillot, T. 1999, P&SS, 47, 1183

Guillot, T., Miguel, Y., Militzer, B., et al. 2018, Natur, 555, 227

Hubbard, W. B. 1974, Icar, 21, 157

Hubbard, W. B. 2013, ApJ, 768, 43

Hubbard, W. B., & Militzer, B. 2016, ApJ, 820, 80

Idini, B., & Stevenson, D. J. 2022, PSJ, 3, 89

Jones, C. A. 2014, Icar, 241, 148

Kaspi, Y., Davighi, J. E., Galanti, E., & Hubbard, W. B. 2016, Icar, 276, 170

Kaspi, Y., & Galanti, E. 2016, ApJ, 820, 91

Kaspi, Y., Galanti, E., Hubbard, W., et al. 2018, Natur, 555, 223

Kaspi, Y., Galanti, E., Showman, A. P., et al. 2020, SSRv, 216, 84

Leconte, J., & Chabrier, G. 2013, NatGe, 6, 347

Ledoux, P. 1947, ApJL, 105, 305

Li, C., Ingersoll, A., Bolton, S., et al. 2020, NatAs, 4, 609

Lindal, G. F., Wood, G. E., Levy, G. S., et al. 1981, JGR, 86, 8721

Liu, S.-F., Hori, Y., Muller, S., et al. 2019, Natur, 572, 355

Lodders, K. 2010, in Principles and Perspectives in Cosmochemistry, Astrophysics and Space Science Proc., ed. A. Goswami & B. E. Reddy (Berlin: Springer), 379

Mahaffy, P. R., Niemann, H. B., Alpert, A., et al. 2000, JGR, 105, 15061

Miguel, Y., Bazot, M., Guillot, T., et al. 2022, A&A, 662, A18

Militzer, B. 2013, PhRvB, 87, 014202

Militzer, B., & Hubbard, W. B. 2013, ApJ, 774, 148

Militzer, B., Wahl, S., & Hubbard, W. B. 2019, ApJ, 879, 78

Moll, R., Garaud, P., Mankovich, C., & Fortney, J. J. 2017, ApJ, 849, 24

Moore, K. M., Yadav, R. K., Kulowski, L., et al. 2018, Natur, 561, 76

Morales, M. A., McMahon, J. M., Pierleonie, C., & Ceperley, D. M. 2013, PhRvL, 110, 065702

Müller, S., Helled, R., & Cumming, A. 2020, A&A, 638, A121

Nettelmann, N., Becker, A., Holst, B., & Redmer, R. 2012, ApJ, 750, 52

Nettelmann, N., Movshovitz, N., Ni, D., et al. 2021, PSJ, 2, 241

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2001, Numerical Recipes in C++ (Cambridge: Cambridge Univ. Press)

Roulston, M., & Stevenson, D. 1995, EOS, 76, 59

Saumon, D., Chabrier, G., & Horn, H. M. V. 1995, ApJS, 99, 713

Seiff, A., Kirk, D. B., Knight, T. C. D., et al. 1997, Sci, 276, 102

Soubiran, F., & Militzer, B. 2015, ApJ, 806, 228

Stevenson, D., & Salpeter, E. 1977, ApJS, 35, 239

Tollefson, J., Wong, M. H., de Pater, I., et al. 2017, Icar, 296, 163

Vazan, A., Helled, R., & Guillot, T. 2018, A&A, 610, id.L14

von Zahn, U., Hunten, D. M., & Lehmacher, G. 1998, JGR, 103, 22815

Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, GeoRL, 44, 4649

Wahl, S. M., Wilson, H. F., & Militzer, B. 2013, ApJ, 773, 95

Wilson, H. F., & Militzer, B. 2010, PhRvL, 104, 121101

Wilson, H. F., & Militzer, B. 2012a, PhRvL, 108, 111101

Wilson, H. F., & Militzer, B. 2012b, ApJ, 745, 54

Wong, M., Mahaffy, P. R., Atreya, S. K., Niemann, H. B., & Owen, T. C. 2004, Icar, 171, 153