# Jupiter's interior from Juno: Equation-of-state uncertainties and dilute core extent*

S. Howard<fnref n="1" />, T. Guillot<fnref n="1" />, M. Bazot<fnref n="2" />,<fnref n="3" />, Y. Miguel<fnref n="4" />,<fnref n="5" />, D. J. Stevenson<fnref n="6" />, E. Galanti<fnref n="7" />, Y. Kaspi<fnref n="7" />, W. B. Hubbard<fnref n="8" />, B. Militzer<fnref n="9" />, R. Helled<fnref n="10" />, N. Nettelmann<fnref n="11" />,<fnref n="12" />, B. Idini<fnref n="6" />, and S. Bolton<fnref n="13" />

Accepted February 6, 2023

## ABSTRACT

**Context.** The Juno mission has provided measurements of Jupiter's gravity field with an outstanding level of accuracy, leading to better constraints on the interior of the planet. Improving our knowledge of the internal structure of Jupiter is key to understanding its formation and evolution but is also important in the framework of exoplanet exploration.

**Aims.** In this study, we investigated the differences between the state-of-the-art equations of state and their impact on the properties of interior models. Accounting for uncertainty on the hydrogen and helium equation of state, we assessed the span of the interior features of Jupiter.

**Methods.** We carried out an extensive exploration of the parameter space and studied a wide range of interior models using Markov chain Monte Carlo (MCMC) simulations. To consider the uncertainty on the equation of state, we allowed for modifications of the equation of state in our calculations.

**Results.** Our models harbour a dilute core and indicate that Jupiter's internal entropy is higher than what is usually assumed from the Galileo probe measurements. We obtain solutions with extended dilute cores, but contrary to other recent interior models of Jupiter, we also obtain models with small dilute cores. The dilute cores in such solutions extend to ~20% of Jupiter's mass, leading to better agreement with formation-evolution models.

**Conclusions.** We conclude that the equations of state used in Jupiter models have a crucial effect on the inferred structure and composition. Further explorations of the behaviour of hydrogen-helium mixtures at the pressure and temperature conditions in Jupiter will help to constrain the interior of the planet, and therefore its origin.

**Key words.** planets and satellites: interiors – planets and satellites: gaseous planets – planets and satellites: composition – planets and satellites: individual: Jupiter – equation of state

## 1. Introduction

Despite the significant improvement to measured gravitational moments provided by Juno (Iess et al. 2018; Durante et al. 2020), the interior of Jupiter remains mysterious. After the two first perijoves of Juno, Wahl et al. (2017) proposed the presence of a dilute core inside the planet: a region above the central compact core where heavy elements are gradually mixed with hydrogen and helium in the envelope. However, most models led to atmospheric abundances that are incompatible with observations. De-

bras & Chabrier (2019) then looked for models compatible with atmospheric abundances and proposed models that require an inward decrease of the abundance of heavy elements (negative Z gradient). While this cannot be ruled out, it seems unlikely from the point of view of long-term stability and formation models. Recently, Militzer et al. (2022) found models with both an atmosphere of protosolar composition of heavy elements and a positive Z gradient, but with large deviations of the gravitational moments requiring specific differential rotation solutions. In Miguel et al. (2022), we found solutions with smaller differential rotation offsets but these required a higher interior entropy than that measured by the Galileo probe.

While all these recent interior models rely on different assumptions, most of them still yield a dilute core inside Jupiter

<footnote n="1">¹ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, France
e-mail: saburo.howard@oca.eu</footnote>

<footnote n="2">² Heidelberg Institute for Theoretical Studies (HITS gGmbH), Schloss-Wolfsbrunnenweg 35,69118 Heidelberg, Germany</footnote>

<footnote n="3">³ CITIES, NYUAD Institute, New York University Abu Dhabi, PO Box 129188, Abu Dhabi, United Arab Emirates,</footnote>

<footnote n="4">⁴ SRON Netherlands Institute for Space Research, Niels Bohrweg 4, 2333 CA Leiden, the Netherlands,</footnote>

<footnote n="5">⁵ Leiden Observatory, University of Leiden, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands,</footnote>

<footnote n="6">⁶ Division of Geological and Planetary Sciences, California Institute of Technology, Pasadena, California 91125, USA,</footnote>

<footnote n="7">⁷ Department of Earth and Planetary Sciences, Weizmann Institute of Science, Rehovot 76100, Israel,</footnote>

<footnote n="8">⁸ Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721, USA,</footnote>

<footnote n="9">⁹ Department of Earth and Planetary Science, University of California, Berkeley, CA 94720, USA,</footnote>

<footnote n="10">¹⁰ Institute for Computational Science, University of Zurich, Winterthurerstr. 190, CH8057 Zurich, Switzerland,</footnote>

<footnote n="11">¹¹ Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA</footnote>

<footnote n="12">¹² Institute of Planetary Research, German Aerospace Center, D-12489 Berlin, Germany</footnote>

<footnote n="13">¹³ Southwest Research Institute, San Antonio, TX 78238, USA</footnote>

<footnote>* Models of Tables E.1. and E.2. are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/cgi-bin/qcat?J/A+A/ and at https://doi.org/10.5281/zenodo.7598377.</footnote>
that is very extended. Debras & Chabrier (2019) and Militzer et al. (2022) find dilute cores that respectively reach 65%-75% and 63% of Jupiter's radius. These values correspond to a dilute core that extends up to respectively ~60% - 75% and ~50% of Jupiter's mass. From the point of view of evolution, Vazan et al. (2018) showed that the outer envelope is mixing efficiently and the outer 60% in mass are of uniform composition after 4.5 Gyr. However, Müller et al. (2020) then modelled the formation of Jupiter using realistic initial entropies and found even more efficient mixing during the planetary evolution where only the inner 20% of the mass is left intact, suggesting that Jupiter's dilute core is not very extended. It has therefore been challenging so far to find agreement over the extent of the dilute core between interior and formation-evolution models of Jupiter, unless an additional process, such as a giant impact, is considered (Liu et al. 2019). In addition, the gravitational imprint of the dilute core in Jupiter's tidal signal registered by Juno (Idini & Stevenson 2022a,b) has lead to increased uncertainty over the extent of the dilute core.

Furthermore, constraining the internal structure of Jupiter requires a good understanding of the behaviour of hydrogen and helium at the pressure and temperature conditions in the planet. Interior models therefore also rely on a key ingredient, which is a hydrogen and helium equation of state (hereafter H-He EOS). Experiments and simulations have been extensively conducted to provide accurate EOSs (see Helled et al. (2020) and references therein for a review) but some uncertainty remains. In addition, Howard & Guillot (2023) recently emphasised the importance of accounting for H-He interactions when calculating EOSs for interior models. The aim of this study is to estimate how extended Jupiter's dilute core is given the current uncertainty on the H-He EOS.

In Section 2, we first explain our methods and the details of our interior models. In Section 3, we compare some H-He EOSs in order to assess the uncertainty that must be accounted for in interior models. Section 4 is devoted to the results, where we present models with original EOSs but also models including a modification of the EOS.

## 2. Methods

Following the method described in Miguel et al. (2022) (see also Guillot et al. 2018), our approach is designed to obtain statistically robust estimates of the properties of Jupiter's interior.

### 2.1. Observed physical characteristics

Interior models satisfy Jupiter's mass, which is obtained through $GM = 1.266865341 \times 10^{17} \text{ m}^3/\text{s}^2$ (Durante et al. 2020). Jupiter's equatorial radius has also been measured $R_{eq} = 71492 \pm 4$ km (Lindal 1992). Jupiter's fast rotation, in 9hr55min29.7s (Davies et al. 1986) flattens the planet and the departure from sphericity can be measured through the gravitational moments, which themselves are measured very accurately by Juno (Durante et al. 2020):

$$J_{2n} = -\frac{1}{MR_{eq}^{2n}} \int \rho(r')^{2n} P_{2n}(\cos\theta) d^3 r' \quad (1)$$

where $R_{eq}$ is Jupiter's equatorial radius, $\rho$ the density, $P_{2n}$ the Legendre polynomials, and $r$ and $\theta$ the radius and the colatitude, respectively. We aim to find models in agreement with

the observed equatorial radius and gravitational moments. As Jupiter exhibits zonal winds in its upper region, we account for the contribution of the latitude-dependent differential rotation to the even gravitational moments (Kaspi et al. 2017, 2018; Guillot et al. 2018). Our interior models must therefore match $J_{2n}^{\text{static}} = J_{2n}^{\text{Juno}} - J_{2n}^{\text{differential}}$ with $J_{2n}^{\text{Juno}}$, that is, the gravitational moments measured by Juno (Durante et al. 2020) and $J_{2n}^{\text{differential}}$, the contribution due to the differential rotation (see Miguel et al. (2022)).

In addition, Jupiter's interior is constrained by observations of its atmosphere. The temperature at the 1 bar level was measured first from Voyager radio occultations (Lindal et al. 1981) and found to be $T_{1bar} = 165 \pm 5$ K. The Galileo entry probe then provided the first and only in situ measurements of Jupiter's atmosphere. The temperature was measured at a single location and found to be $T_{1bar} = 166.1 \pm 0.8$ K (Seiff et al. 1998). Nonetheless, Gupta et al. (2022) questioned the representability of this temperature over a wider range of latitudes and longitudes. After reassessing the Voyager radio occultations, these latter authors found that the temperature at 1 bar could reach a value of $T_{1bar} = 170.3 \pm 3.8$ K. Thus, $T_{1bar}$ is either fixed or allowed to vary in our models.

Galileo also provided measurements of the atmospheric composition. We set the helium abundance in our models through the measurement of $Y_{atm} = Y_1/(X_1 + Y_1) = 0.238$ (von Zahn et al. 1998) and $Y_{proto} = 0.277$ (Serenelli & Basu 2010). The abundances of heavy elements were also measured but the particularly low value of the water abundance from Galileo led to an updated analysis in the equatorial region with the Juno mission (Li et al. 2020). Figure 1 sums up the mass fractions of the most abundant heavy elements in Jupiter's atmosphere. We compute the total mass fraction of ices considering only CH$_4$, NH$_3$, H$_2$O, and H$_2$S and using a root mean square to estimate its uncertainty. As we do not have measurements of the refractory elements in the atmosphere of Jupiter, we show two values of the total abundance of heavy elements $Z_{tot}$, assuming their abundance to be zero and three times the protosolar value. Depending on the amount of rocks we considered, the total abundance of heavy elements is likely to be between 1.6 ($Z_{tot} = 0.024$) and 4.0 ($Z_{tot} = 0.061$) times the protosolar abundance. Our models presented in Section 4 use an abundance of 1.3 times the protosolar abundance ($Z = 0.02$) which is close to the lower limit of $Z_{tot}$.

### 2.2. Details and calculations of the interior models

We keep almost the same parameterisation as in Miguel et al. (2022), with the main difference being that we only consider dilute core models here. Standard three-layer models can be considered as a particular case of dilute core models in which the dilute core is of uniform composition and extends all the way to the helium phase separation region. We note that, as the helium phase separation appears late in the evolution (Mankovich & Fortney 2020), it is unlikely that a very strong change in composition can appear in that region. Here we choose not to consider three-layer models and focus on dilute core models. This latter type of model is composed of a central compact core made exclusively of heavy elements, a dilute core region where the heavy elements are gradually distributed outwards, a metallic hydrogen layer, and an outer molecular hydrogen layer. Figure 2 shows the typical distribution of heavy elements in a model with a dilute core. The mass fractions of heavy elements in the molecular hydrogen layer, the metallic hydrogen layer, and the dilute core are respectively monitored with the parameters $Z_1, Z_2$, and
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Abundances of heavy elements in the atmosphere of Jupiter. Here, 'ices' refers to the metallicity of the atmosphere considering only ices (CH$_4$, NH$_3$, H$_2$O and H$_2$S). Two values of Z$_tot$ are shown, which correspond to the metallicity of the atmosphere considering ices (volatiles) and assuming either no rocks (refractories) or an enrichment of three times the protosolar value for rocks. Abundances of methane and hydrogen sulfide are taken from Wong et al. (2004), and ammonia and water abundances are taken from Li et al. (2020). Protosolar abundances are taken from Asplund et al. (2021). The protosolar mass fraction of heavy elements is 0.0154. For ices, we calculated the protosolar value considering C, N, O, S, and Ne and included additional elements with the value from Lodders et al. (2009). Concerning rocks, we only considered Mg, Si, and Fe.</figcaption></figure>

Z$_dilute$. As a strong change in composition is not expected at the location where helium rain occurs, we fix Z$_1$ = Z$_2$. We stress that this assumption can affect the gravity harmonics, especially the higher ones, which are more sensitive to the parameters of the outer envelope (Z$_1$, Z$_2$). The extent of the dilute core is controlled by m$_dilute$, which is the normalised mass where the dilute core merges with the metallic hydrogen layer. The mass fraction of heavy elements in the dilute core region is defined by:

$$Z(m) = Z_1 + \frac{Z_{\text{dilute}} - Z_1}{2} \left[ 1 - \text{erf} \left( \frac{m - m_{\text{dilute}}}{\delta m_{\text{dil}}} \right) \right], \quad (2)$$

where δm$_dil$ is the slope of the gradual change in the mass fraction of heavy elements, and is set to 0.075. The compact core is only made of heavy elements (rocks) and its mass is defined by M$_core$. We define two useful quantities, M$_Z,dil*$ and M$_Z,env*$, in order to assess how predominant the dilute core is, but also to provide estimates of the amount of heavy elements that needs to be accreted onto the compact core during the formation of the planet.

Finally, we note that we use the code CEPAM (Guillot & Morel 1995) to calculate interior models. We compute the gravitational moments using the theory of figures (Zharkov & Trubitsyn 1978) at fourth order. The gravitational moments are then calibrated (Guillot et al. 2018) using the accuracy of the concentric MacLaurin spheroid (Hubbard 2012, 2013). More details can be found in Miguel et al. (2022).

<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Distribution of heavy elements for a model with a dilute core. The three hashed areas correspond to the mass of the compact core (only made of heavy elements), the mass of heavy elements in the dilute core region excluding the area where Z < Z$_1$, and the mass of heavy elements in the rest of the envelope (where Z < Z$_1$). Z$_1$ is the mass fraction of heavy elements in the outer envelope. Z$_dilute$ is the maximum mass fraction in the dilute core. m$_dilute$ monitors the extent of the dilute core in terms of mass.</figcaption></figure>

## 2.3. The MCMC approach

We aim to explore a great number of interior models to understand which ones best represent Jupiter's internal structure. We followed the method used in Miguel et al. (2022), which is based on Bazot et al. (2012). We took a Bayesian approach based on a Markov chain Monte Carlo (MCMC) code. The equatorial radius as well as the gravitational moments are considered as the data while measurements of the parameters give us a priori information. The MCMC samples models calculated with CEPAM. We computed the posterior density function and therefore the likelihood at various points of the parameter space. We were then able to obtain the joint posterior probability densities. Additional details about the method can be found in the aforementioned papers. Furthermore, in Section 4.1 we discuss the priors chosen for the parameters of our models and a table with all the parameters is provided in Appendix.

# 3. Uncertainty on the H-He equation of state

## 3.1. A diversity of equations of state

Jupiter being mostly composed of hydrogen and helium, interior models of the planet require the use of an appropriate EOS for these chemical elements. The SCvH95 (Saumon et al. 1995) EOS was commonly applied for studying Jupiter and extrasolar giant planets (Thorngren et al. 2016). However, in the last decade, with the improvement of shock Hugoniot data analyses and the development of ab initio simulations, new H-He EOSs have emerged (Militzer & Hubbard 2013; Chabrier et al. 2019; Mazevet et al. 2022). To compare these EOSs, we derived adiabats for pure hydrogen-helium mixtures (Y = 0.245) by simply integrating the adiabatic gradient starting from 1 bar, 166.1 K (i.e. Jupiter's conditions). For the sake of comparison, heavy elements have not been included but they would only slightly affect the adiabats (see Helled (2018)). Table 1 details
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Adiabats obtained from different EOSs and corresponding to a homogeneous model, with no compact core, with $Y = 0.245$. Table 1 lists the details of the EOSs. Below are shown the contribution functions of the gravitational moments of order 0 to 8 centred at their peak. Their extent corresponds to their full width at half maximum (see Guillot (2005) for more details).</figcaption></figure>

the specifics of the EOSs used to calculate these adiabats. The MGF16+MH13 EOS was derived by Miguel et al. (2016) from the ab initio simulation data of Militzer & Hubbard (2013). We built MH13* by fitting an adiabat provided by Burkhard Militzer (private communication), based on the ab initio EOS of Militzer & Hubbard (2013). We used a polynomial function $g$ to fit the residuals between the MGF16+MH13 adiabat and the one provided by BM. We then perturbed the MGF16+MH13 EOS so that $\rho_{MH13*} = (1-g)\rho_{MGF16+MH13}$. The comparison between MGF16+MH13 and MH13* is detailed in Section 3.2. More information about the derivation of CD21 can be found in Debras et al. (2021). The HG23+CMS19 and HG23+MLS22 EOSs both make use of the tables of the non-linear mixing effects provided in Howard & Guillot (2023). We note that all these EOSs, except SCvH95, account for the interactions between hydrogen and helium particles. Still, MGF16+MH13, CD21, and MH13* include the non-ideal mixing effects but they remain fixed and equal to that calculated by Militzer & Hubbard (2013) for a single composition ($Y = 0.245$). On the other hand, HG23+CMS19 and HG23+MLS22 include the H-He interactions and remain valid for any composition of the mixture (see Howard & Guillot (2023)). For convenience, as $P \propto \rho^2$ in Jupiter's interior (Hubbard 1975), Figure 3 compares values of $\rho/\sqrt{P}$. The SCvH95 EOS is much less dense than the other, more recent EOSs, ranging from 0.1 to 10 Mbar. (Quantum Monte Carlo simulations (Mazzola et al. 2018) yield even denser hydrogen (between 0.3 and 2.6 Mbar) than found by EOSs obtained from density functional theory.) We find a maximum relative deviation between the different EOSs (except SCvH95) that amounts to 5.5% (around 0.03 Mbar). The discrepancy between the various EOSs gives us an estimate of the uncertainty on the EOS to be accounted for in Jupiter models.

## 3.2. Interpolation uncertainties

The differences between the EOSs seen in Fig. 3 are surprising, because in that parameter space, with the exception of SCvH95,

they are all based on the results obtained by the same ab initio calculations from Militzer & Hubbard (2013). In order to understand where these differences come from, we must examine the way the EOS tables are constructed. Two tables are available in Militzer & Hubbard (2013). Table 1 provides the thermodynamic quantities obtained from the density functional molecular dynamics (DFT-MD) calculations. This latter was directly used by Miguel et al. (2016) and grafted to the SCvH95 EOS to construct the MGF16+MH13 table. On the other hand, Table 2 from Militzer & Hubbard (2013) provides coefficients for a free-energy fit from which one can calculate all thermodynamic quantities. This free-energy fit is used by Militzer et al. (2022) and forms the basis of the MH13 EOS (similar to our MH13* EOS). Both the CD21 EOS (Chabrier & Debras 2021) and the nonideal mixing tables of Howard & Guillot (2023) use this EOS and have to interpolate the results with the SCvH95 EOS in the molecular regime. To assess the uncertainty due to the choice of table, but also due to the way we interpolate through the table, we derived adiabats from both Table 1 and Table 2 from Militzer & Hubbard (2013). To do so, we used a one-dimensional interpolation to evaluate pressure and temperature at a typical value of entropy for Jupiter ($7.078061 k_B e^{-1}$), for each density value. This procedure is straightforward for Table 1, but for Table 2 we followed the procedure prescribed by Militzer & Hubbard (2013) before deriving the adiabats. We then tried three different types of interpolation (linear, quadratic, and cubic) when calculating the adiabats. Figure 4 shows the different extent in parameter space of the two tables from Militzer & Hubbard (2013) and how different choices, in particular on the order of the interpolation, affect the resulting adiabat. Table 2 is slightly extended compared to Table 1 and provides density and entropy for temperatures between 1000 and 80,000 K, and pressures between 0.1 and 300 Mbar. The maximum deviation between adiabats calculated from Table 1 and Table 2 is of the order of 2%. The order of the interpolation brings a maximum deviation that amounts to 1.3%. These uncertainties lead to the differences between the MGF16+MH13 and the MH13* EOSs. However, at pressures of lower than 0.1 Mbar, there are also discrepancies between both EOSs that are certainly due to the combination of the DFT-MD calculations of Militzer & Hubbard (2013) with the SCvH95 EOS of Saumon et al. (1995). The construction of an EOS is very sensitive to the merging of several tables, particularly around the regions where the tables are connected. This may explain the high values of density around 0.03 Mbar and the slightly lower densities at $P < 0.01$ Mbar for CD21, HG23+CMS19, and HG23+MLS22 (see Fig. 3). Furthermore, we can see that the points of Table 1—displayed in Fig. 4—are sparse, particularly between 0.1 and 1 Mbar, at densities relevant to the region used to derive a Jupiter adiabat, affecting the accuracy of the interpolation through the table.

## 3.3. A thermodynamically consistent modification of the EOS

With the uncertainty on the EOSs in hand, we want to account for it in our interior models. To do so, we need a function perturbing an EOS. Initially, we simply used a Gaussian function to perturb the density profile of our models (similarly to Nettelmann et al. (2021)). Nonetheless, an EOS cannot be perturbed freely. Indeed, any variations of the EOS should satisfy the limits of a thermodynamical potential. The Helmholtz free energy, which is usually where an EOS comes from, is relatively well known at low and high density regimes. Low densities correspond to the regime of an ideal gas and the free energy is known from experimental data and statistical mechanics. High densi-
<table><caption>**Table 1.** H-He EOSs used in interior models. Concerning the heavy elements, the EOSs used for ices and rocks are respectively those of water and dry sand from Lyon & Johnson (1992).</caption><thead><tr><th>Name of the H-He EOS</th><th>H EOS</th><th>He EOS</th><th>H-He interactions</th><th>Notes/References</th></tr></thead><tbody><tr><td>MGF16+MH13</td><td>MGF16-H</td><td>SCvH95-He</td><td>Included in the H EOS but fixed at Y = 0.245</td><td>Derived by Miguel et al. (2016) from Militzer & Hubbard (2013)</td></tr><tr><td>MH13*</td><td>--</td><td>--</td><td>--</td><td>Adjusted from MGF16+MH13 to fit a Y = 0.245 adiabat from Militzer et al. (2022)</td></tr><tr><td>CD21</td><td>CD21-H</td><td>CMS19-He</td><td>Included in the H EOS but fixed at Y = 0.245</td><td>Derived by Chabrier & Debras (2021) from Militzer & Hubbard (2013)</td></tr><tr><td>HG23+CMS19</td><td>CMS19-H</td><td>CMS19-He</td><td>HG23</td><td>H and He EOSs from Chabrier et al. (2019) and non-ideal mixing effects from Howard & Guillot (2023)</td></tr><tr><td>HG23+MLS22</td><td>MLS22-H</td><td>CMS19-He</td><td>HG23</td><td>H EOS from Mazevet et al. (2022) and non-ideal mixing effects from Howard & Guillot (2023)</td></tr><tr><td>SCvH95</td><td>SCvH95-H</td><td>SCvH95-He</td><td>None</td><td>H and He EOSs from Saumon et al. (1995)</td></tr></tbody></table>

ties correspond to the regime well above the metallisation pressure where hydrogen is fully ionised and the free energy here is known from theory and simulations. Here, we use the internal energy because $P = -\left(\frac{dU}{dV}\right)_S$ can lead to an integral constraint on legitimate density changes of the EOS. If we know the internal energy at low and high regimes for a given entropy, that is, $U(\rho_1, S)$ and $U(\rho_2, S)$, respectively, we can obtain an expression of the differences between these two terms, which is to be conserved by perturbations along the adiabat:

$$ \Delta U = U(\rho_2, S) - U(\rho_1, S) = \int_{\rho_1}^{\rho_2} \frac{P}{\rho^2} d\rho. \qquad (3) $$

If $\delta\rho$ corresponds to a slight density modification, then we have

$$ \Delta U_{\text{modif. EOS}} - \Delta U_{\text{orig. EOS}} = \int_{\rho_1}^{\rho_2} \frac{P}{(\rho + \delta\rho)^2} d\rho - \int_{\rho_1}^{\rho_2} \frac{P}{\rho^2} d\rho, \quad (4) $$

Using the approximation $P \propto \rho^2$, we then get

$$ \Delta U_{\text{modif. EOS}} - \Delta U_{\text{orig. EOS}} \propto \int_{\rho_1}^{\rho_2} \frac{\delta\rho}{\rho} d\rho = 0. \qquad (5) $$

This difference must be 0 if $\Delta U$ in Eq. (3) is to be conserved. Equation (5) provides a constraint on how we are allowed to perturb an EOS. Hence, we need to choose an appropriate function to perturb the EOS while verifying this constraint. This function, denoted $f$ here, comes from the equation

$$ \rho_{\text{modif}} = \rho + \delta\rho = \rho(1 + f). \qquad (6) $$

Therefore, we need to find $f$ so that

$$ \int_{\rho_1}^{\rho_2} f d\rho = 0. \qquad (7) $$

Using again $P \propto \rho^2$ and changing variables, we need to choose an $f$ that satisfies

$$ \int_{P_1}^{P_2} f(\log_{10}(P)) \sqrt{P} d\log_{10}(P) = 0. \qquad (8) $$

We naturally define $f$ as

$$ f = K \frac{d\rho}{\sqrt{P}} \exp\left(-\left[\frac{\log_{10}(P/P_{\text{modif}})}{\Delta P}\right]^2\right) \text{erf}\left(\frac{\log_{10}(P/P_{\text{modif}})}{\Delta P}\right). \qquad (9) $$

The function $f$ is composed of a Gaussian and an error function, and includes a division by the square root of the pressure, with $d\rho$ being the amplitude of the density modification. From Eq. (5), we infer that to satisfy this integral constraint, a density reduction at a certain pressure will imply a density increase at another pressure. To modify an EOS and obtain this trend, we use the product of a Gaussian and an error function. And to properly satisfy the integral constraint from Eq. (5) after changing the integral as a function of the density into an integral as a function of $\log_{10}(P)$, we need to multiply by the square root of the pressure. This function depends on three parameters: $P_{\text{modif}}$ corresponds to the pressure (in cgs units) at which the density modification is applied, $\Delta P$ controls the width of the modification; and $d\rho$ is the amplitude of the density change. The constant $K$ is in units of square root of pressure divided by density and is set to $K = 1.04 \times 10^6 \sqrt{\text{dyn.cm}^{-2}\cdot\text{g}^{-1}\cdot\text{cm}^3}$. Figure 5 shows how the integral constraint from Eq. (5) is satisfied for two different models. One model simply uses a Gaussian function to modify the density profile and clearly does not respect the integral constraint. On the other hand, the second model, which uses the function $f$ defined above, satisfies this constraint well, because the $\Delta U$ difference falls close to zero at high pressures. More precisely, this value at high pressures is not exactly zero due to the perturbation theory approximation applied to Eq. (4) where $\delta\rho/\rho$ is assumed small. But overall, there is a significant difference between how the integral constraint is respected between the two models of Figure 5, and the effort of satisfying this constraint must be underlined. We stress that our function $f$ was naturally but arbitrarily chosen; there are certainly other functions that could satisfy Eq. (5).

### 3.4. Priors on the modification of the EOS

For MCMC runs in which we allow modifications of the EOS, we use priors on $d\rho$, $P_{\text{modif}}$, and $\Delta P$, as defined in Eq. (9). The priors are either Gaussian or uniform (a more detailed discussion can be found in Section 4.2) with boundaries set to avoid physically inconsistent modifications. We set $P_{\text{modif}}$ between $10^{11.5}$ and $10^{12.5}$ dyn.cm$^{-2}$ (which correspond to 0.3 and 3 Mbar), as a preliminary study shows us that a density reduction followed by an increase in the density at higher pressure is preferred over a density increase followed by a decrease. We set $\Delta P$ between 0.2 and 0.8 so that low ($< 10^{-2}$ Mbar) and high ($> 10$ Mbar) pres-
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Comparison between the MH13* and MGF16+MH13 EOSs. *Top panel.* Comparison of Table 1 and Table 2 from Militzer & Hubbard (2013). *Bottom panel.* Comparison of adiabats obtained with the MH13* and MGF16+MH13 EOSs and of points obtained by direct interpolation of the above tables. The $(P, \rho/\sqrt{P})$ points of Table 1 from Militzer & Hubbard (2013) are shown as blue crosses. The locations at which the entropy is equal to $7.078061 \ k_B\text{el.}^{-1}$ (a typical value of entropy for Jupiter) according to the interpolations in the two tables are shown with different symbols, as labelled.</figcaption></figure>

sures are not significantly affected by the density modification. Allowing $P_{\text{modif}}$ and $\Delta P$ to pass the chosen boundaries would allow modifications of the EOS that are not consistent with the constraints. The boundaries of $d\rho$ are set to allow a change of amplitude in density that does not exceed 10%. Using a random sampling, Figure 6 shows the possible modifications that can be allowed from the original EOSs with the priors and boundaries we chose. When plotting the density according to the pressure, we can see that the differences are small even when we perturb the adiabats. Table A.1 sums up the priors we chose to run further MCMC simulations.

## 4. Results

Constraints on the interior structure and composition of Jupiter are derived as follows. First, as described in the following subsection, in order to avoid results in which the solutions are too far from the observational constraints, we chose to fix some parameters. On this basis, we ran models for the different EOSs and outer envelope metallicity $Z_1$. In the following subsections, we present two types of results: (a) models with the original EOSs, variable $T_{1\text{bar}}$ and $Z_1 = 0.02$; and (b) models with modified EOSs, a value of $T_{1\text{bar}}$ equal to either the Galileo value of 166.1 K or the upper limit from Gupta et al. (2022), 174.1 K, and values of $Z_1 = 0.02, 0.029$, and $0.035$ (1.3 to 2.3 × protosolar). Our results confirm and extend those of Miguel et al. (2022). When using the same hypotheses, we obtain the same results as
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** Comparisons of the ΔU difference for a model modified with only a Gaussian function (blue) and a model modified with the defined function *f* (see Eq. (9)) composed of a Gaussian and an error function (orange). Top panel. Comparison of the adiabats ($P_{\text{modif}} = 10^{11} \text{ dyn.cm}^{-2}$, $\Delta P=0.5$, $d\rho=-0.06$ for the blue curve & $P_{\text{modif}} = 10^{11.5} \text{ dyn.cm}^{-2}$, $\Delta P=0.6$, $d\rho=-0.05$). Bottom panel. Difference between the ΔU of the modified models and that of the reference model (MH13).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 6.** Possible adiabats (gray lines) that can be obtained when modifying the MGF16+MH13, MH13*, CD21, HG23+CMS19, and HG23+MLS22 EOSs with the chosen priors (see Table A.1).</figcaption></figure>

Militzer et al. (2022) (see Appendix D), and results that are consistent with those presented by Debras & Chabrier (2019).

## 4.1. Choice of priors

The optimisation problem is characterised by two issues: (1) The constraints on Jupiter's gravitational moments are extremely tight, meaning that only a tiny fraction of the parameter space allows for successful models, and (2) the high density of recent

H-He EOSs imply that in most cases, fitting only Jupiter's radius and $J_2$ value would require nonphysical negative $Z_1$ or core mass values. In practice, this means that some parameters can be led to values that are significantly offset from their prior.

To assess the influence of the priors, we ran four different simulations (using the MGF16+MH13 EOS) with our MCMC code and focus on three parameters: $Y_{\text{proto}}$, $Y_{\text{atm}}$, $T_{\text{1bar}}$ and two data: $J_4$, $J_6$. In the first run, we set all parameters to vary freely using a Gaussian prior (mean $\mu$, standard deviation $\sigma$) centred on the value from observations. In the second, third, and fourth runs, we respectively fix one, two, and three of the aforementioned parameters. We note that the accuracy on the $J_{2n}$ is clearly higher compared to that on the other parameters. $\sigma/\mu$ is equal to 0.01% and 0.2% for $J_4$ and $J_6$ (0.0002% for $J_2$), respectively, while it is equal to 2% for the three other parameters. Figure 7 shows the posterior distribution of $Y_{\text{proto}}$, $Y_{\text{atm}}$, $T_{\text{1bar}}$, $J_4$, and $J_6$ for the four different runs. When setting all parameters free, the MCMC code samples models well around the mean values of the $J_{2n}$ measured by Juno (accounting for the influence of differentially rotating winds; see Miguel et al. (2022)). However, $Y_{\text{proto}}$, $Y_{\text{atm}}$, and $T_{\text{1bar}}$ are at $2\sigma$ or $3\sigma$ from the mean value of their respective prior. When fixing $Y_{\text{proto}}$, the sampled values of $J_6$ are now at $1\sigma$ from the observed mean value and sampled values of $Y_{\text{atm}}$ are at $7-8\sigma$ from Galileo's measurement. When fixing $Y_{\text{proto}}$ and $Y_{\text{atm}}$, the $J_6$ fit is slightly poorer than the previous run. However, $T_{\text{1bar}}$ has now a $4\sigma$ difference from the mean value of the prior. Finally, when fixing $Y_{\text{proto}}$, $Y_{\text{atm}}$, and $T_{\text{1bar}}$, $J_4$ differs by more than $20\sigma$ from Juno's measurement and $J_6$ differs by $5\sigma$.

This shows that to satisfy the observational constraints on parameters such as $Y_{\text{proto}}$ and $Y_{\text{atm}}$, we need to fix them instead of using a Gaussian prior, because the uncertainties on the gravitational moments dominate. Otherwise, these parameters will be further than $1\sigma$ from the observed value. As the observational constraint is looser on $T_{\text{1bar}}$ because of the question of latitudinal dependency —which raises the possibility that Jupiter's deep temperature may be higher than the Galileo probe reference value (see Section 2.1)—, only $T_{\text{1bar}}$ will be set as a free parameter in some of the runs presented in this paper.

## 4.2. Runs with a modified EOS

For MCMC runs in which we allow for modifications of the EOS, we first chose Gaussian priors on $P_{\text{modif}}$, $\Delta P$, and $d\rho$ (see Sections 3.3 and 3.4). Our goal was to penalise models with a substantial change of the EOS and favour models with only a slight modification of the EOS; as for the other parameters discussed in the previous section, this led to large deviations of the EOS. Figure 8 compares the adiabats of models including EOS modifications, with a uniform prior on $P_{\text{modif}}$ and either uniform or Gaussian priors on $\Delta P$ and $d\rho$. For the Gaussian priors, the parameters for $\Delta P$ were $\mu = 0.5$ and $\sigma = 0.02$, and the parameters for $d\rho$ were $\mu = 0$. and $\sigma = 0.01$. The difference between $\rho/\sqrt{P} - P$ profiles obtained after modification of the EOS is subtle between runs with Gaussian and uniform priors. This slight difference between the modified adiabats leads to a better fit of the data (equatorial radius, gravitational moments) for models obtained with uniform priors. When using Gaussian priors, most of the models are at $2\sigma$ or $3\sigma$ from the observed equatorial radius and $J_4$ and are at $4\sigma$ to $5\sigma$ from the mean value assumed for the prior on $P_{\text{He}}$ (see Fig. B.1). As the modified adiabats are very similar and the agreement with the observational constraints is slightly better, we present results obtained with uniform priors.
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Posterior distributions of $Y_{\text{proto}}$, $Y_{\text{atm}}$, $T_{\text{lbar}}$, $J_4$, and $J_6$ for four different MCMC runs using the MGF16+MH13 EOS. The blue histograms correspond to a run where the four parameters are free. In the orange run, $Y_{\text{proto}}$ is fixed. In the green run, $Y_{\text{proto}}$ and $Y_{\text{atm}}$ are fixed. In the red run, $Y_{\text{proto}}$, $Y_{\text{atm}}$, and $T_{\text{lbar}}$ are fixed. When histograms show a black vertical solid line, this indicates that a prior was set and was centred at the value from observations ($Y_{\text{proto}} = 0.277 \pm 0.006$, $Y_{\text{atm}} = 0.238 \pm 0.005$, $T_{\text{lbar}} = 165 \pm 4$ K, $J_4 \cdot 10^6 = -586.53 \pm 0.0836$ and $J_6 \cdot 10^6 = 34.18 \pm 0.07682$). The gray areas correspond to the standard deviation of the prior (1 sigma).</figcaption></figure>

As previously mentioned, the modification of the initial MGF16+MH13 EOS is relatively substantial, with a change in density that can reach up to ~ 11% in amplitude (see Fig. 8). Figure 9 shows the modifications of the EOS for models (using uniform priors) at respectively $T_{\text{lbar}} = 166.1$ K and 174.1 K. At higher $Z_1$, the modifications of the EOS occur at higher pressures but in any case, the amplitude remains significant (between 6 and 11%). In addition, these changes to the EOS are likely to be incompatible with Hugoniot data (Knudson & Desjarlais 2017). We therefore provide the results with modified EOS as a way to test the robustness of the solutions, but we generally focus on results using the original EOSs.

### 4.3. Surface temperature $T_{\text{lbar}}$ and helium transition pressure $P_{\text{He}}$

As mentioned in the preamble of Section 4, we present two sets of models: with original EOSs and with a modification of the EOS. Here, we focus on the 1 bar temperature, which prescribes the entropy inside Jupiter, and on the pressure where helium rain occurs (Stevenson & Salpeter 1977). The latter sets the limit between the molecular hydrogen (helium-poor) and the metallic hydrogen (helium-rich) layers. Figure 10 shows the values of these two parameters for the two types of interior models sampled by our MCMC code. Interior models with the original EOSs (and $Z_{\text{atm}} = 0.02$) all yield a 1 bar temperature, which is higher than the value measured by Galileo, which ranges from 171 to 188 K. In particular, with the MGF16+MH13 EOS, we obtain a $T_{\text{lbar}}$ of between 180 and 188 K, while we were obtaining a $T_{\text{lbar}}$ of between 175 and 183 K for $Z_{\text{atm}} = 0.0153$ (protosolar value) in Miguel et al. (2022). Therefore, a Z abundance of $1.3 \times \text{protosolar}$ instead of $1 \times \text{protosolar}$ in the outer envelope leads to a 5 K increase in $T_{\text{lbar}}$ to fit the $J_{2n}$ (see Section 4.4). Only models using MH13* or HG23+MLS22 seem to be in line with the upper end of the temperature provided by Gupta et al.

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Adiabats obtained for models with a modification of the EOS. The black solid line corresponds to the original MGF16+MH13 EOS. Red shows results obtained with Gaussian priors on $\Delta P$ ($\mu = 0.5$ and $\sigma = 0.02$) and $d\rho$ ($\mu = 0$. and $\sigma = 0.01$). Blue shows results obtained with uniform priors. The prior on $P_{\text{modif}}$ remains uniform in both cases. The dashed lines correspond to the adiabat obtained with the mean values of $P_{\text{modif}}$, $\Delta P$, and $d\rho$ (see Section 3.3) of a subsample of 100 models randomly drawn from the MCMC output. We compute the standard deviation ($\sigma$) of the 100 adiabats and the envelopes show the adiabats of the 1 $\sigma$ spread from the mean modified adiabat (dashed line). Here, $T_{\text{lbar}}$ is fixed at 166.1 K.</figcaption></figure>

(2022). Overall, all models present a high 1 bar temperature, which could correspond to a deep entropy in line with a hotter interior due to a potential superadiabaticity (Guillot 1995; Leconte & Chabrier 2012). The models using original EOSs exhibit a helium transition pressure of between 0.8 and 4.5 Mbar, which is in agreement with the values obtained by simulations and experiments (see Lorenzen et al. (2011); Morales et al. (2013); Schöttler & Redmer (2018); Brygoo et al. (2021)). For all EOSs, we
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Adiabats obtained for models with modified EOS. Top Panel. $T_{1\text{bar}}$ is fixed at 166.1 K. The light blue area shows results for $Z_1 = 0.02$ (1.3 × the protosolar value) while dark blue shows results for $Z_1 = 0.0286$ (1.9 × the protosolar value). Bottom Panel. $T_{1\text{bar}}$ is fixed at 174.1 K. The yellow area shows results for $Z_1 = 0.02$ (1.3 × the protosolar value) while gray shows results for $Z_1 = 0.035$ (2.3 × the protosolar value). The black solid line corresponds to the original MGF16+MH13 EOS. Other details of the figure can be found in the caption of Fig. 8.</figcaption></figure>

can distinguish two ensembles of solutions: one with high $P_{\text{He}}$ corresponding to models with a compact core of a few earth masses (1-6 $M_\oplus$) and a highly extended dilute core ($m_{\text{dilute}}$ between 0.4 and 0.6) and a second one with lower $P_{\text{He}}$ corresponding to models with almost no compact core and a less extended dilute core ($m_{\text{dilute}}$ between 0.15 and 0.45) (more details in Section 4.6). The difference in $T_{1\text{bar}}$ between the two ensembles of solutions is only of 2-3 K. Therefore, accurate constraint of the pressure at which helium rain occurs could help to characterize the dilute core and to determine the atmospheric entropy of Jupiter.

Concerning models with a modification of the EOS, we calculate two subsets of models where we fix $T_{1\text{bar}}$ at 166.1 or 174.1 K. For each temperature, we present results for two values of the abundance of heavy elements in the outer envelope $Z_1$. At $T_{1\text{bar}}$ = 166.1 K, models with $Z_1$ = 0.02 have $P_{\text{He}}$ values concentrated between 3 and 4.5 Mbar. With $Z_1$ = 0.029, models tend to high transition pressures, around 6 Mbar, and are far from fitting the observed equatorial radius and gravitational moments (see Section 4.4). At $T_{1\text{bar}}$ = 174.1 K, we obtain values of $P_{\text{He}}$ of between 1.5 and 3.5 Mbar for $Z_1$ = 0.02 and of between 3 and 4.5 Mbar for $Z_1$ = 0.035, which are both close to what is expected.

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Temperature at 1 bar vs. helium transition pressures for two types of models. Left panel. Models using original EOSs. $Z_1$ = 0.02. Right panel. Models with a modification of the EOS. The initial EOS that has been modified is MH13. We present two subsets of models: with $T_{1\text{bar}}$ = 166.1 K and $T_{1\text{bar}}$ = 174.1 K. The black dotted line corresponds to the 1 bar temperature measured by the Galileo probe.</figcaption></figure>

### 4.4. Equatorial radius $R_{\text{eq}}$ and gravitational moments $J_{2n}$

Here, we examine the fit of our models to the gravitational moments measured by Juno and accounting for differential rotation. Figure 11 shows the equatorial radius and the gravitational moments obtained with our interior models. All models with original EOSs can reproduce the equatorial radius and all the gravitational moments except $J_6$. We find solutions for MGF16+MH13 that can match $J_6$ corrected by differential rotation. For the four other EOSs, the sampled values of $J_6$ are in the 2 - 3 $\sigma$ range. We can see a correlation between $T_{1\text{bar}}$ and $J_6$: models using an EOS that yields higher $T_{1\text{bar}}$ present lower values of $J_6$. With MGF16+MH13, $J_6 \times 10^6$ is between 34.1 and 34.3. With MH13*, CD21 and HG23+CMS19, $J_6 \times 10^6$ is between 34.2 and 34.4. With HG23+MLS22, $J_6 \times 10^6$ is between 34.3 and 34.5. Militzer et al. (2022) found a value of $J_6 \times 10^6$ of 34.47 for $T_{1\text{bar}}$ = 166.1 K using their EOS from Militzer & Hubbard (2013).

Concerning models allowing for a modification of the EOS, at $T_{1\text{bar}}$ = 166.1 K, we manage to find models matching $R_{\text{eq}}$ and all $J_{2n}$ when $Z_1$ = 0.02. However, at $Z_1$ = 0.029, we can no longer fit $R_{\text{eq}}$ or $J_{2n}$. These models have an equatorial radius that is several sigma below the observational constraint but we still retain them to test the robustness of our results. We then set $T_{1\text{bar}}$ to 174.1 K and find models reproducing $R_{\text{eq}}$ and all $J_{2n}$, even for $Z_1$ = 0.035 (2.3 × protosolar).

### 4.5. Heavy-element distribution

We now compare the distribution of heavy elements in our models. Figure 12 shows the heavy elements masses defined in Section 2.2. Models with original EOSs have a total mass of heavy elements of between 18 and 33 $M_\oplus$ and a compact core of less than 6 $M_\oplus$. These results are in line with those obtained by Miguel et al. (2022), showing that increasing the Z abundance in the outer envelope from 1 × protosolar to 1.3 × protosolar does not lead to a drastic change in the distribution of heavy elements. $M_{Z,\text{dil}*}$ is between 10 and 25 $M_\oplus$, which is larger than $M_{Z,\text{env}*}$ by up to a factor of 4. Hence, models with no modification of the
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Equatorial radius and gravitational moments for two types of models. Left panels. Models using original EOSs. $Z_1 = 0.02$. Right panels. Models with a modification of the EOS. The initial EOS that has been modified is MH13. We present two subsets of models: with $T_{1bar} = 166.1$ K and $T_{1bar} = 174.1$ K. The circles with a J correspond to the measurements of the gravitational moments by Juno (Durante et al. 2020). The black error bars correspond to the gravitational moments corrected by differential rotation (see 2.1).</figcaption></figure>

EOS have most of their heavy elements in the dilute core region rather than in the rest of the envelope.

Allowing for modifications of the EOS generally leads to a lower total mass of heavy elements, mostly between 12 and 20 $M_\oplus$. This is due to modifications of the EOS that make the adiabats (hence the H-He mixture) denser at depth. The mass of the compact core does not exceed 8 $M_\oplus$ for these models. $M_{Z,dil*}$ is similar in all of our four cases: models are concentrated around a region where $M_{Z,dil*} \sim 5-7 M_\oplus$. However, $M_{Z,env*}$ clearly depends on the value of $Z_1$. For $Z_1 = 0.02$, $M_{Z,env*}$ is around 6 $M_\oplus$, for $Z_1 = 0.029$, $M_{Z,env*}$ is around 9 $M_\oplus$, and for $Z_1 = 0.035$, $M_{Z,env*}$ is around 11 $M_\oplus$. Therefore, our models with a modified EOS do not lead to a dilute core that is predominant in heavy elements compared to the rest of the envelope, contrary to what we find for models with original EOSs.

We stress that the total masses of heavy elements inferred here are lower limits: The presence of compositional gradients implies that parts of the interior may be super-adiabatic because it is Ledoux-stable, double-diffusive (see Leconte & Chabrier 2012), or stable to moist convection (Guillot 1995; Leconte et al. 2017), meaning that the interior could be warmer and thus retain more heavy elements than calculated here. For example, Militzer et al. (2022) estimated that a doubling of the central temperature of Jupiter would increase the mass of heavy elements from 25 to 42 $M_\oplus$.

<figure><img src="image_12.png" /><figcaption>**Fig. 12.** Masses of heavy elements in our interior models. Top panels. Mass of the compact core vs. total mass of heavy elements in Jupiter. Bottom panels. $M_{Z,dil*}$ vs. $M_{Z,env*}$ (see Section 2.2).</figcaption></figure>

## 4.6. Dilute core characteristics

A key question is the extent of this dilute core, which connects interior models of Jupiter and formation and evolution models. Recent interior models from Militzer et al. (2022) and Debras & Chabrier (2019) suggest considerably extended dilute cores that respectively reach 63% and 65%-75% of Jupiter's radius. These values correspond to ~50% and ~60%-75% of Jupiter's mass, respectively, which can be compared to the value of our $m_{dilute}$ parameter. We note that the comparison is not exact because it is affected by the functional form chosen for the dilute core (see Section 2.2), but the effect is minor: In our case, the added heavy element mass fraction in the dilute core drops from 50% of its maximal value at $m = m_{dilute}$ to only 8% at $m_{dilute} + \delta m_{dil}$ with $\delta m_{dil} = 0.075$.

The preferred model from Militzer et al. (2022) would have a value of $m_{dilute}$ ~ 0.36 in our parameterisation. Figure 13 shows the values of $m_{dilute}$ found for our models. Using original EOSs, we obtain $m_{dilute}$ of between 0.15 and 0.6 (the dilute core extends from ~15% to ~60% of Jupiter's total mass). We confirm that we find models with very extended dilute cores as in Debras & Chabrier (2019) and Militzer et al. (2022). But we are also finding models with relatively narrow dilute cores (down to ~15% of Jupiter's mass). These solutions with relatively small, dilute cores are in better agreement with the mixing and evolution calculations of Müller et al. (2020), which yield a dilute core that does not exceed 20% of Jupiter's mass, resulting in a dilute core that extends only up to Jupiter's inner ~60 $M_\oplus$. This leads to a formation scenario that is consistent for both Jupiter and Saturn (see Guillot et al. 2022), because Saturn is likely to harbour a dilute core that extends to 52–60 $M_\oplus$ of the planet's total mass (Mankovich & Fuller 2021).

Concerning models with a modification of the EOS, we find a few models with very extended dilute cores but the majority yield $m_{dilute} < 0.15$. We suspect that these results are spurious. The changes in the H-He EOS lead to an increase in density at high pressures that can mimic the effect of a dilute core. As the H-He mixture is denser in the dilute core region, less heavy
<figure><img src="image_13.png" /><figcaption>**Fig. 13.** Z<sub>dilute</sub> vs. m<sub>dilute</sub> for two types of models. Z<sub>dilute</sub> is the maximum mass fraction of heavy elements in the dilute core region while m<sub>dilute</sub> controls the extent of the dilute core in terms of mass (see Section 2.2). Left panel. Models using original EOSs. Z<sub>1</sub> = 0.02. Right panel. Models with a modification of the EOS. The initial EOS that has been modified is MH13. We present two subsets of models: with T<sub>1bar</sub> = 166.1 K and T<sub>1bar</sub> = 174.1 K. The thin gray dotted line corresponds to m<sub>dilute</sub> = 0.2, which corresponds approximately to the extent of the dilute core of 20% of Jupiter's mass predicted by formation models from Müller et al. (2020).</figcaption></figure>

elements can be added, which leads to these low values of m<sub>dilute</sub>
(and also M<sub>Z,dil*</sub>, see Section 4.5).

**5. Conclusion**

We explored a wide variety of interior models of Jupiter con-
strained by all available observations and using available EOSs.
Our models assume the presence of a central dense compact
core, a dilute core of variable extent and heavy-element com-
position, and an outer envelope of uniform Z composition. The
helium phase separation is modelled as a jump in helium abun-
dance in the Mbar pressure range.

While high-pressure experiments and ab initio calculations
have led to significant improvements, H-He EOSs remain a
source of uncertainty when modelling Jupiter's interior. We ob-
serve a range in the adiabatic density profiles of up to 5% at pres-
sures ranging from 10 kbar to 10 Mbar. We interpret these varia-
tions as resulting from changes between different EOS tables and
from different interpolation choices in relatively sparsely popu-
lated tables.

An important source of uncertainty results from our poor
knowledge of Jupiter's complex atmosphere and the possibil-
ity of a higher entropy than generally assumed. By allowing the
T<sub>1bar</sub> parameter to vary (see Miguel et al. 2022), we obtain mod-
els that fit all constraints for all EOSs used. The values of T<sub>1bar</sub>
obtained range from 171 K to 188 K, significantly higher than
166.1 K from the Galileo probe (von Zahn et al. 1998), but within
164K - 174K, which is the range of values obtained by Gupta
et al. (2022) from a reanalysis of Voyager's radio occultations.
Interestingly, MH13* and HG23+MLS22 (see Table 1), the two
EOSs that lead to the lowest T<sub>1bar</sub> values, yield the highest values
of J<sub>6</sub>, which are slightly outside the range expected from differ-
ential rotation. Conversely, the MGF16+MH13 EOS leads to the
smallest values of J<sub>6</sub>, well within expectations, but the largest
values of T<sub>1bar</sub>.

In all cases, we obtain a dilute core of heavy-element mass of
between 10 and 25 M<sub>⊕</sub>, confirming the result obtained by Miguel
et al. (2022) that Jupiter's envelope is inhomogeneous. The range
of values is also fully compatible with the results obtained by
Militzer et al. (2022) and Debras & Chabrier (2019). The mass of
the compact core ranges between 0 and 6 M<sub>⊕</sub>. The total mass of
heavy elements that we find ranges from 18 to 33 M<sub>⊕</sub>. However,
we must stress that, given the possibility of (perhaps significant)
superadiabatic regions (see Guillot 1995; Leconte & Chabrier
2012; Leconte et al. 2017, and Section 4.5), these masses are
lower limits.

Our dilute cores are characterised by a global mass fraction
of heavy elements of between 0.02 and 0.27 (in addition to the
envelope heavy element mass fraction Z₁ ~ 0.02) and extend
from ~ 15% to ~ 60% of Jupiter's total mass. We reiterate that
the exact extent of the dilute core will depend on the shape of
its compositional gradient. These solutions therefore encompass
those of Debras & Chabrier (2019) and Militzer et al. (2022), but
also allow for small, dilute cores. These solutions with small, di-
lute cores are compatible with the formation–evolution models
of Müller et al. (2020), which suggest that the outer 80% in mass
should be fully mixed by convection, leaving a primordial dilute
core extending only up to Jupiter’s inner ~ 60 M⊕. This result
is also promising, in light of interior models for Saturn, which
indicate that a dilute core extends to 52 – 60 M⊕ of the planet’s
total mass (Mankovich & Fuller 2021). This could lead to a for-
mation scenario that is consistent for both Jupiter and Saturn (see
Guillot et al. 2022).

However, there is an important caveat to consider: As pointed
out by Wahl et al. (2017) and Debras & Chabrier (2019), and
confirmed by all further modelling efforts, interior models of
Jupiter constrained by Juno's gravitational moments favour so-
lutions with small values of Z in the outer envelope. While our
solutions are calculated by imposing Z₁ = 0.02, this represents
a bare minimum, given all spectroscopic constraints (see Fig. 1).
When imposing higher values of Z₁, we find an increasingly
more difficult situation, with solutions departing from the con-
straints on R<sub>eq</sub> and J₂ and/or modifications to the EOSs that were
too important and most likely incompatible with the experimen-
tal constraints. Interestingly, a similar situation arises for Saturn
(Mankovich & Fuller 2021), indicating that we may be missing
an important physical ingredient.

Several directions of research could lead to significant im-
provements in our understanding of Jupiter's interior structure
and composition. For example, progress in the analysis of Juno
microwave radiometer data to infer abundances of ammonia and
water as well as temperatures as a function of depth and altitude
(Li et al. 2020) will help us to understand heat transport and the
composition of the deep atmosphere (see also Stevenson et al.
2022). Future radio occultations with Juno should further test
observational constraints on temperature and shape. Improve-
ments on EOSs, both experimentally and numerically, with par-
ticular emphasis on the hydrogen–helium mixture at pressures
between 10 kbar and 10 Mbar and near Jupiter’s adiabat (temper-
ature from 1000 K to 20,000 K on this pressure range) would
be extremely valuable. Finally, while indications of the presence
of normal modes of Jupiter exist (Gaulme et al. 2011; Durante
et al. 2022), their identification from dedicated observational ef-
forts (Gonçalves et al. 2019; Shaw et al. 2022) would be an ex-
tremely powerful tool for fully constraining the interior structure
and composition.

***Acknowledgements.*** The authors thank the Juno Interior Working Group for useful discussions and comments. This research was carried out at the Observatoire
de la Côte d'Azur under the sponsorship of the Centre National d'Etudes Spatiales.

## References

Asplund, M., Amarsi, A. M., & Grevesse, N. 2021, A&A, 653, A141

Bazot, M., Bourguignon, S., & Christensen-Dalsgaard, J. 2012, MNRAS, 427, 1847

Brygoo, S., Loubeyre, P., Millot, M., et al. 2021, Nature, 593, 517

Chabrier, G. & Debras, F. 2021, ApJ, 917, 4

Chabrier, G., Mazevet, S., & Soubiran, F. 2019, ApJ, 872, 51

Davies, M. E., Abalakin, V. K., Bursa, M., et al. 1986, Celestial Mechanics, 39, 103

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

Debras, F., Chabrier, G., & Stevenson, D. J. 2021, ApJ, 913, L21

Durante, D., Guillot, T., Iess, L., et al. 2022, Nature Communications, 13, 4632

Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e86572

Gaulme, P., Schmider, F. X., Gay, J., Guillot, T., & Jacob, C. 2011, A&A, 531, A104

Gonçalves, I., Schmider, F. X., Gaulme, P., et al. 2019, Icarus, 319, 795

Guillot, T. 1995, Science, 269, 1697

Guillot, T. 2005, Annual Review of Earth and Planetary Sciences, 33, 493

Guillot, T., Fletcher, L. N., Helled, R., et al. 2022, arXiv e-prints, arXiv:2205.04100

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

Guillot, T. & Morel, P. 1995, A&AS, 109, 109

Gupta, P., Atreya, S. K., Steffes, P. G., et al. 2022, PSJ, 3, 159

Helled, R. 2018, in Oxford Research Encyclopedia of Planetary Science, 175

Helled, R., Mazzola, G., & Redmer, R. 2020, Nature Reviews Physics, 2, 562

Howard, S. & Guillot, T. 2023, arXiv e-prints, arXiv:2302.07902

Hubbard, W. B. 1975, Soviet Ast., 18, 621

Hubbard, W. B. 2012, ApJ, 756, L15

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

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

Idini, B. & Stevenson, D. J. 2022b, PSJ, 3, 11

Iess, L., Folkner, W. M., Durante, D., et al. 2018, Nature, 555, 220

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

Kaspi, Y., Guillot, T., Galanti, E., et al. 2017, Geophys. Res. Lett., 44, 5960

Knudson, M. D. & Desjarlais, M. P. 2017, Phys. Rev. Lett., 118, 035501

Leconte, J. & Chabrier, G. 2012, A&A, 540, A20

Leconte, J., Selsis, F., Hersant, F., & Guillot, T. 2017, A&A, 598, A98

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

Lindal, G. F. 1992, AJ, 103, 967

Lindal, G. F., Wood, G. E., Levy, G. S., et al. 1981, J. Geophys. Res., 86, 8721

Liu, S.-F., Hori, Y., Müller, S., et al. 2019, Nature, 572, 355

Lodders, K., Palme, H., & Gail, H. P. 2009, Landolt B&ouml;l;rnstein, 4B, 712

Lorenzen, W., Holst, B., & Redmer, R. 2011, Phys. Rev. B, 84, 235109

Lyon, S. P. & Johnson, J. D. 1992, LANL Report, LA-UR-92-3407

Mankovich, C. R. & Fortney, J. J. 2020, ApJ, 889, 51

Mankovich, C. R. & Fuller, J. 2021, Nature Astronomy, 5, 1103

Mazevet, S., Licari, A., & Soubiran, F. 2022, A&A, 664, A112

Mazzola, G., Helled, R., & Sorella, S. 2018, Phys. Rev. Lett., 120, 025701

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

Miguel, Y., Guillot, T., & Fayon, L. 2016, A&A, 596, A114

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

Militzer, B., Hubbard, W. B., Wahl, S., et al. 2022, PSJ, 3, 185

Morales, M. A., Hamel, S., Caspersen, K., & Schwegler, E. 2013, Phys. Rev. B, 87, 174105

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

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

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

Schöttler, M. & Redmer, R. 2018, Phys. Rev. Lett., 120, 115703

Seiff, A., Kirk, D. B., Knight, T. C. D., et al. 1998, J. Geophys. Res., 103, 22857

Serenelli, A. M. & Basu, S. 2010, ApJ, 719, 865

Shaw, C. L., Gulledge, D. J., Swindle, R., Jefferies, S. M., & Murphy, N. 2022, Frontiers in Astronomy and Space Sciences, 9, 768452

Stevenson, D. J., Bodenheimer, P., Lissauer, J. J., & D'Angelo, G. 2022, PSJ, 3, 74

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

Thorngren, D. P., Fortney, J. J., Murray-Clay, R. A., & Lopez, E. D. 2016, ApJ, 831, 64

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

von Zahn, U., Hunten, D. M., & Lehmacher, G. 1998, J. Geophys. Res., 103, 22815

Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649

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

Zharkov, V. N. & Trubitsyn, V. P. 1978, Physics of planetary interiors
# Appendix A: Priors in MCMC simulations

Table A.1 lists the priors used for the parameters of our MCMC calculations.

# Appendix B: Comparison between runs with Gaussian and uniform priors

Figure B.1 compares the posterior distributions of two MCMC simulations: using Gaussian or uniform priors on the parameters to modify the EOS (see Section 4.2).

# Appendix C: Corner plots of models with and without modification of the EOS

Figures C.1, C.2, C.3, C.4, and C.5 show the posterior distributions of the MCMC simulations using original EOSs, respectively, MGF16+MH13, MH13*, CD21, HG23+CMS19, and HG23+MLS22.

Figures C.6 and C.7 show the posterior distributions of the MCMC simulations using modified EOS, with $T_{1\text{bar}} = 166.1$ K and respectively for $Z_1 = 0.02$ and $Z_1 = 0.0286$. Figures C.8 and C.9 show the posterior distributions of the MCMC simulations using modified EOS, with $T_{1\text{bar}} = 174.1$ K and respectively for $Z_1 = 0.02$ and $Z_1 = 0.035$.

# Appendix D: Comparison with Militzer et al. (2022)

We ran MCMC simulations to reproduce the results obtained by Militzer et al. (2022). To do so, we changed the values of the gravitational moments around which the MCMC is sampling models. We used the gravitational moments of the interior model of Militzer et al. (2022) (see their Table 1). We used the same properties: no compact core, $Z_1 = 0.0153$, $T_{1\text{bar}} = 166.1$ K, and the MH13* EOS. Figure D.1 shows the posterior distributions we obtain. We find models with similar properties to Model A from Militzer et al. (2022): the same gravitational moments, the same $P_{\text{He}}$, and comparable characteristics for the dilute core.

# Appendix E: Subsample of models

We extracted one single model from each MCMC simulation using original EOSs (see Fig. C.1, C.2, C.3, C.4, C.5) and list them in Table E.1 and Table E.2.
<table><caption>**Table A.1.** Parameters explored in our MCMC calculations for dilute core models.</caption><thead><tr><th>Parameter</th><th>Distribution</th><th>Lower bound</th><th>Upper bound</th><th>μ</th><th>σ</th></tr></thead><tbody><tr><td><i>M</i><sub>core</sub> (M<sub>⊕</sub>)</td><td>Uniform</td><td>0</td><td>24</td><td>—</td><td>—</td></tr><tr><td><i>P</i><sub>He</sub> (Mbar)</td><td>Normal</td><td>0.8</td><td>9</td><td>3</td><td>0.5</td></tr><tr><td><i>T</i><sub>jump</sub> (K)</td><td>Uniform</td><td>0</td><td>2000</td><td>—</td><td>—</td></tr><tr><td><i>Z</i><sub>dilute</sub><sup>rock</sup></td><td>Uniform</td><td>0</td><td>0.5</td><td>—</td><td>—</td></tr><tr><td><i>Z</i><sub>dilute</sub><sup>ice</sup></td><td>Uniform</td><td>0</td><td>0.5</td><td>—</td><td>—</td></tr><tr><td><i>m</i><sub>dilute</sub></td><td>Uniform</td><td>0</td><td>0.6</td><td>—</td><td>—</td></tr><tr><td><i>T</i><sub>1bar</sub> (K)</td><td>Normal</td><td>135</td><td>215</td><td>165</td><td>4</td></tr><tr><td><i>P</i><sub>modif</sub> (dyn.cm<sup>-2</sup>)</td><td>Uniform</td><td>10<sup>11.5</sup></td><td>10<sup>12.5</sup></td><td>—</td><td>—</td></tr><tr><td>Δ<i>P</i></td><td>Uniform</td><td>0.2</td><td>0.8</td><td>—</td><td>—</td></tr><tr><td><i>dρ</i></td><td>Uniform</td><td>-0.1</td><td>0.1</td><td>—</td><td>—</td></tr></tbody></table>

**Notes.** The parameter is given in the first column, the corresponding distribution in the second, and the lower and upper bounds in the third and fourth. When relevant, the mean and the standard deviation of the truncated normal are given in columns five and six. The prior on $T_{1bar}$ is used in models where $T_{1bar}$ is not fixed. The priors on $P_{modif}$, $\Delta P$, and $d\rho$ are used in models allowing for a modification of the EOS. We recall that $Y_{atm} = 0.238$ and $Y_{proto} = 0.277$.

<table><caption>**Table E.1.** Comparison of the parameters of selected models extracted from MCMC simulations using original EOSs.</caption><thead><tr><th>EOS</th><th><i>M</i><sub>core</sub> (M<sub>⊕</sub>)</th><th><i>P</i><sub>He</sub> (Mbar)</th><th><i>m</i><sub>dilute</sub></th><th><i>Z</i><sub>dilute</sub></th><th><i>T</i><sub>jump</sub> (K)</th><th><i>T</i><sub>1bar</sub> (K)</th><th><i>M</i><sub>Z env*</sub> (M<sub>⊕</sub>)</th><th><i>M</i><sub>Z dil*</sub> (M<sub>⊕</sub>)</th><th><i>M</i><sub>Z tot</sub> (M<sub>⊕</sub>)</th></tr></thead><tbody><tr><td>MGF16+MH13</td><td>0.0305</td><td>2.12</td><td>0.237</td><td>0.214</td><td>60.3</td><td>185.3</td><td>6.52</td><td>14.46</td><td>20.98</td></tr><tr><td>MH13*</td><td>0.0288</td><td>2.05</td><td>0.317</td><td>0.175</td><td>82.3</td><td>176.1</td><td>6.52</td><td>15.51</td><td>22.03</td></tr><tr><td>CD21</td><td>0.0115</td><td>1.64</td><td>0.280</td><td>0.186</td><td>16.1</td><td>183.3</td><td>6.52</td><td>14.57</td><td>21.08</td></tr><tr><td>HG23+CMS19</td><td>0.231</td><td>1.66</td><td>0.324</td><td>0.177</td><td>9.10</td><td>182.1</td><td>6.52</td><td>16.05</td><td>22.57</td></tr><tr><td>HG23+MLS22</td><td>0.466</td><td>2.06</td><td>0.441</td><td>0.180</td><td>60.0</td><td>174.6</td><td>6.52</td><td>22.19</td><td>28.71</td></tr></tbody></table>

**Notes.** Models are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/cgi-bin/qcat?J/A+A/ and at https://doi.org/10.5281/zenodo.7598377.

<table><caption>**Table E.2.** Comparison of the equatorial radius and the gravitational moments of the same models as in Table E.1.</caption><thead><tr><th>EOS</th><th><i>R</i><sub>eq</sub> (km)</th><th><i>J</i><sub>2</sub> × 10<sup>6</sup></th><th><i>J</i><sub>4</sub> × 10<sup>6</sup></th><th><i>J</i><sub>6</sub> × 10<sup>6</sup></th><th><i>J</i><sub>8</sub> × 10<sup>6</sup></th><th><i>J</i><sub>10</sub> × 10<sup>6</sup></th></tr></thead><tbody><tr><td>MGF16+MH13</td><td>71487.6</td><td>14695.42</td><td>-586.649</td><td>34.211</td><td>-2.4548</td><td>0.2011</td></tr><tr><td>MH13*</td><td>71492.1</td><td>14695.62</td><td>-586.622</td><td>34.339</td><td>-2.4750</td><td>0.2035</td></tr><tr><td>CD21</td><td>71491.8</td><td>14695.57</td><td>-586.611</td><td>34.291</td><td>-2.4676</td><td>0.2027</td></tr><tr><td>HG23+CMS19</td><td>71491.4</td><td>14695.53</td><td>-586.559</td><td>34.309</td><td>-2.4704</td><td>0.2030</td></tr><tr><td>HG23+MLS22</td><td>71491.2</td><td>14695.65</td><td>-586.625</td><td>34.436</td><td>-2.4904</td><td>0.2054</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Fig. B.1.** Posterior distributions obtained with a modification of the EOS. Red shows results obtained with Gaussian priors. Blue shows results obtained with uniform priors. $T_{1\text{bar}}$ is fixed at 166.1 K and $Z_1 = 0.02$. The black points correspond to the measured $J_{2n}$ by Juno. The black error bars correspond to Juno's measurements accounting for differential rotation for the $J_{2n}$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. C.1.** Posterior distributions obtained with the MGF16+MH13 EOS, where $T_{1\text{bar}}$ is a free parameter, $Z_1 = 0.02$ (1.3 × the protosolar value). The black points correspond to the measured $J_{2n}$ by Juno. The black error bars correspond to Juno's measurements accounting for differential rotation for the $J_{2n}$ and Galileo's measurement for $T_{1\text{bar}}$.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. C.2. Same as Fig. C.1 but with the MH13* EOS. The red star shows the Militzer et al. (2022) preferred (static) model.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Fig. C.3. Same as Fig. C.1 but with the CD21 EOS.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. C.4. Same as Fig. C.1 but with the HG23+CMS19 EOS.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. C.5. Same as Fig. C.1 but with the HG23+MLS22 EOS.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. C.6.** Posterior distributions obtained with a modification of the EOS, where $T_{1\text{bar}}$ is fixed at 166.1 K, $Z_1 = 0.02$ (1.3 × the protosolar value). The black points correspond to the measured $J_{2n}$ by Juno. The black error bars correspond to Juno's measurements accounting for differential rotation for the $J_{2n}$.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>Fig. C.7. Same as Fig. C.6 but with $Z_1 = 0.0286$ (1.9 × the protosolar value)</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Fig. C.8.** Posterior distributions obtained with a modification of the EOS, where $T_{1\text{bar}}$ is fixed at 174.1 K, $Z_1 = 0.02$ (1.3 × the protosolar value). The black points correspond to the measured $J_{2n}$ by Juno. The black error bars correspond to Juno's measurements accounting for differential rotation for the $J_{2n}$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Fig. C.9. Same as Fig. C.8 but with $Z_1 = 0.035$ (2.3 × the protosolar value)</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. D.1.** Posterior distributions obtained with MH13*. The red star shows the Militzer et al. (2022) preferred (static) model. *M*<sub>core</sub> = 0, *T*<sub>1bar</sub> is fixed at 166.1 K and *Z*<sub>1</sub> = 0.0153. The black error bars correspond to Juno's measurements accounting for differential rotation for the *J*<sub>2n</sub>.</figcaption></figure>