© The Authors 2023

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. Debras & 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 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 EOS s (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 EOS s 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 Sect. 2, we first explain our methods and the details of our interior models. In Sect. 3, we compare some H-He EOS s 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 EOS s 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 × 10¹⁷ m³ s⁻² (Durante et al. 2020). Jupiter’s equatorial radius has also been measured R_eq = 71492 ± 4 km (Lindal 1992). Jupiter’s fast rotation, in 9 h 55 min 29.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}{M{R}_{\text{eq}}^{2n}}{\displaystyle \int \rho {\left(r^{\prime }\right)}^{2n}{P}_{2n}\left(\cos \theta \right){\text{d}}^3{r}^{\prime }},$$(1)

where R_eq is Jupiter’s equatorial radius, ρ the density, Р_2n the Legendre polynomials, and r and θ 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 ± 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 Т_1bar = 166.1 ± 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 ± 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₁/(X₁ + Y₁) = 0.238 (von Zahn et al. 1998) and Y_proyo = 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₄, NH₃, H₂O, and H₂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 Sect. 4 use an abundance of 1.3 times the protosolar abundance (Z = 0.02) which is close to the lower limit of Z_tot.

Fig. 1 Abundances of heavy elements in the atmosphere of Jupiter. Here, ‘ices’ refers to the metallicity of the atmosphere considering only ices (CH4, NH3, H2O and H2S). Two values of Ztot 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.

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 < Z1, and the mass of heavy elements in the rest of the envelope (where Z < Z1). Z1 is the mass fraction of heavy elements in the outer envelope. Zdilute is the maximum mass fraction in the dilute core, mdilute monitors the extent of the dilute core in terms of mass.

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₁, Z₂, and Z_dilute. As a strong change in composition is not expected at the location where helium rain occurs, we fix Z₁ = Z₂. 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₁, Z₂). 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\left(m\right)=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],$$(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).

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 Sect. 4.1 we discuss the priors chosen for the parameters of our models and a table with all the parameters is provided in Appendix A.

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 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 (priv. comm.), based on the ab initio EOS of Militzer & Hubbard (2013). We used a polynomial function ɡ to fit the residuals between the MGF16+MH13 adiabat and the one provided by BM. We then perturbed the MGF16+MH13 EOS so that ρ_MH13* = (1 – g)ρ_MGF16+MH13. The comparison between MGF16+MH13 and MH13^* is detailed in Sect. 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 ∝ ρ² in Jupiter’s interior (Hubbard 1975), Fig. 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 Tables 1 and 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_bel.⁻¹), 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 Tables 1 and 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.

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).

Fig. 4 Comparison between the MH13* and MGF16+MH13 EOSs. Top panel: comparison of Tables 1 and 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 $\left(P,\rho /\sqrt{P}\right)$ points of Table 1 from Militzer & Hubbard (2013) are shown as blue crosses. The locations at which the entropy is equal to 7.078061 kbel.−1 (a typical value of entropy for Jupiter) according to the interpolations in the two tables are shown with different symbols, as labelled.

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 Nettehnann et al. 2021). Nonetheless, an EOS cannot be perturbed freely. Indeed, any variations of the EOS should satisfy the limits of a thermo-dynamical 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 densities 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{\text{d}U}{\text{d}V}\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(ρ₁,S) and U(ρ₂,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\left({\rho }_2,S\right)-U\left({\rho }_1,S\right)={\displaystyle {\int }_{{\rho }_1}^{{\rho }_2}\frac{P}{{\rho }^2}\text{d}\rho .}$$(3)

If δρ corresponds to a slight density modification, then we have $$\Delta U_{\text{modif.EOS}}-\Delta U_{\text{orig.EOS}}={\displaystyle {\int }_{{\rho }_1}^{{\rho }_2}\frac{P}{{\left(\rho +\delta \rho \right)}^2}\text{d}\rho -{\displaystyle {\int }_{{\rho }_1}^{{\rho }_2}\frac{P}{{\rho }^2}\text{d}\rho .}}$$(4)

Using the approximation P ∝ ρ², we then get $$\text{Δ}{U}_{\text{modif.EOS}}-\text{Δ}{U}_{\text{orig.EOS}}\propto {\displaystyle {\int }_{{\rho }_1}^{{\rho }_2}\frac{\delta P}{\rho }\text{d}\rho =0.}$$(5)

This difference must be 0 if Δ 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 \left(1+f\right).$$(6)

Therefore, we need to find f so that $${\displaystyle {\int }_{{\rho }_1}^{{\rho }_2}f\text{d}\rho =0.}$$(7)

Using again P ∝ ρ² and changing variables, we need to choose an f that satisfies $${\displaystyle {\int }_{{\rho }_1}^{{\rho }_2}f\left({\log }_{10}\left(P\right)\right)\sqrt{P}\text{d}{\log }_{10}\left(P\right)=0.}$$(8)

We naturally define f as $$f=K\frac{\text{d}\rho }{\sqrt{P}}\exp \left(-{\left[\frac{{\log }_{10}\left(P/P_{\text{modif}}\right)}{\text{Δ}P}\right]}^2\right)\text{erf}\left(\frac{{\log }_{10}\left(P/P_{\text{modif}}\right)}{\text{Δ}P}\right).$$(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ρ 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₁₀(P), we need to multiply by the square root of the pressure. This function depends on three parameters: P_modif corresponds to the pressure (in cgs units) at which the density modification is applied, ΔP controls the width of the modification; and dρ 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}{\text{.cm}}^{-2}{\text{g}}^{-1}{\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 Δ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 δρ/ρ is assumed small. But overall, there is a significant difference between how the integral constraint is respected between the two models of Fig. 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).

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 (Pmodif = 1011 dyn cm−2, ΔP = 0.5, dρ = −0.06 for the blue curve & Pmodif = 1011.5 dyn cnr−2, ΔP = 0.6, dρ = −0.05). Bottom panel: difference between the ΔU of the modified models and that of the reference model (MH13).

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ρ, P_modif, and ΔP, as defined in Eq. (9). The priors are either Gaussian or uniform (a more detailed discussion can be found in Sect. 4.2) with boundaries set to avoid physically inconsistent modifications. We set P_modif between 10^11.5 and 10^12.5 dyn. cm⁻² (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 ΔP between 0.2 and 0.8 so that low (<10⁻² Mbar) and high (>10 Mbar) pressures are not significantly affected by the density modification. Allowing P_modif and ΔP to pass the chosen boundaries would allow modifications of the EOS that are not consistent with the constraints. The boundaries of dρ are set to allow a change of amplitude in density that does not exceed 10%. Using a random sampling, Fig. 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.

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).

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₁. In the following subsections, we present two types of results: (a) models with the original EOSs, variable T_1bar and Z₁ = 0.02; and (b) models with modified EOSs, a value of T_1bar 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₁ = 0.02, 0.029, and 0.035 (1.3–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 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₂ value would require nonphysical negative Z₁ 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_proto, Y_atm, T_1bar and two data: J₄, J₆. In the first run, we set all parameters to vary freely using a Gaussian prior (mean µ, standard deviation σ) 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, σ/µ is equal to 0.01 and 0.2% for J₄ and J₆ (0.0002% for J₂), respectively, while it is equal to 2% for the three other parameters. Figure 7 shows the posterior distribution of Y_proto, Y_atm, T_1bar, J₄, and J₆ 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_proto, Y_atm, and T_1bar are at 2σ or 3σ from the mean value of their respective prior. When fixing Y_proto, the sampled values of J₆ are now at 1σ from the observed mean value and sampled values of Y_atm are at 7–8σ from Galileo’s measurement. When fixing Y_proto and Y_atm, the J₆ fit is slightly poorer than the previous run. However, T_1bar has now a 4σ difference from the mean value of the prior. Finally, when fixing Y_рrоto, Y_atm, and T_1bar, J₄ differs by more than 20σ from Juno’s measurement and J₆ differs by 5σ.

This shows that to satisfy the observational constraints on parameters such as Y_proto and Y_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σ from the observed value. As the observational constraint is looser on T_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 Sect. 2.1) –, only T_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_modif ΔP, and dρ (see Sects. 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_modif and either uniform or Gaussian priors on ΔP and dρ. For the Gaussian priors, the parameters for ΔP were µ = 0.5 and σ = 0.02, and the parameters for dρ were µ = 0. and σ = 0.01. The difference between $\rho /\sqrt{\left(P\right)}-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σ or 3σ from the observed equatorial radius and J₄ and are at 4σ to 5σ from the mean value assumed for the prior on P_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.

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_1bar = 166.1 and 174.1 K. At higher Z₁, 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.

Fig. 7 Posterior distributions of Yproto, Yatm. T1bar, J4, and J6 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, Yproto is fixed. In the green run, Yproto and Yatm are fixed. In the red run, Yproto, Yatm, and T1bar 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 (Yproto = 0.277 ± 0.006, Yatm = 0.238 ± 0.005, T1bar = 165 ± 4K, J4.106 = −586.53 ± 0.0836 and J6.106 = 34.18 ± 0.07682). The gray areas correspond to the standard deviation of the prior (1 sigma).

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 ΔP (µ = 0.5 and σ = 0.02) and dρ (µ = 0. and σ = 0.01). Blue shows results obtained with uniform priors. The prior on Pmodif remains uniform in both cases. The dashed lines correspond to the adiabat obtained with the mean values of Pmodif, ΔP, and dρ (see Sect. 3.3) of a subsample of 100 models randomly drawn from the MCMC output. We compute the standard deviation (σ) of the 100 adiabats and the envelopes show the adiabats of the 1σ spread from the mean modified adiabat (dashed line). Here, T1bar is fixed at 166.1 K.

4.3 Surface temperature T_1bar and helium transition pressure P_He

As mentioned in the preamble of Sect. 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_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_1bar of between 180 and 188 K, while we were obtaining a T_1bar of between 175 and 183 К for Z_atm = 0.0153 (protosolar value) in Miguel et al. (2022). Therefore, a Z abundance of 1.3 × protosolar instead of 1 × protosolar in the outer envelope leads to a 5 K increase in T_1bar to fit the J_2n (see Sect. 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. (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 can distinguish two ensembles of solutions: one with high P_He corresponding to models with a compact core of a few earth masses (1–6 М_⊕) and a highly extended dilute core (m_dilute between 0.4 and 0.6) and a second one with lower P_He corresponding to models with almost no compact core and a less extended dilute core (m_dilute between 0.15 and 0.45; more details in Sect. 4.6). The difference in T_1bar 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 characterise 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_1bar 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₁. At T_1bar = 166.1 K, models with Z₁ = 0.02 have P_не values concentrated between 3 and 4.5 Mbar. With Z₁ = 0.029, models tend to high transition pressures, around 6 Mbar, and are far from fitting the observed equatorial radius and gravitational moments (see Sect. 4.4). At T_1bar = 174.1 K, we obtain values of P_не of between 1.5 and 3.5 Mbar for Z₁ = 0.02 and of between 3 and 4.5 Mbar for Z₁ = 0.035, which are both close to what is expected.

Fig. 9 Adiabats obtained for models with modified EOS. Top panel: T1bar is fixed at 166.1 K. The light blue area shows results for Z1 = 0.02 (1.3 × the protosolar value) while dark blue shows results for Z1 = 0.0286 (1.9 × the protosolar value). Bottom panel: T1bar is fixed at 174.1 K. The yellow area shows results for Z1 = 0.02 (1.3 × the protosolar value) while gray shows results for Z1 = 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.

Fig. 10 Temperature at 1 bar vs. helium transition pressures for two types of models. Left panel: models using original EOSs. Z1 = 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 T1bar = 166.1 and T1bar = 174.1 K. The black dotted line corresponds to the 1 bar temperature measured by the Galileo probe.

Fig. 11 Equatorial radius and gravitational moments for two types of models. Left panels: models using original EOSs. Z1 = 0.02. Right panels: models with a modification of the EOS. The inital EOS that has been modified is MH13. We present two subsets of models: with T1bar = 166.1 and T1bar = 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 Sect. 2.1).

4.4 Equatorial radius R_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₆. We find solutions for MGF16+MH13 that can match J₆ corrected by differential rotation. For the four other EOSs, the sampled values of J₆ are in the 2–3σ range. We can see a correlation between T_1bar and J₆: models using an EOS that yields higher T_1bar present lower values of J₆. With MGF16+MH13, J₆ × 10⁶ is between 34.1 and 34.3. With MH13*, CD21 and HG23+CMS19, J₆ × 10⁶ is between 34.2 and 34.4. With HG23+MLS22, J₆ × 10⁶ is between 34.3 and 34.5. Militzer et al. (2022) found a value of J₆ × 10⁶ of 34.47 for T_1bar = 166.1 K using their EOS from Militzer & Hubbard (2013).

Concerning models allowing for a modification of the EOS, at T_1bar = 166.1 K, we manage to find models matching R_eq and all J_2n when Z₁ = 0.02. However, at Z₁ = 0.029, we can no longer fit R_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_1bar to 174.1 K and find models reproducing R_eq and all J_2n even for Z₁ = 0.035 (2.3 × protosolar).

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: MZ,dil* vs. МZ,env* (see Sect. 2.2).

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 Sect. 2.2. Models with original EOSs have a total mass of heavy elements of between 18 and 33 M_⊕ and a compact core of less than 6 M_⊕. 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,dil* is between 10 and 25 M_⊕, which is larger than M_Z,env*. by up to a factor of 4. Hence, models with no modification of the 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_⊕. 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_⊕ for these models. M_Z,dil* is similar in all of our four cases: models are concentrated around a region where M_Z,dil* ~ 5–7 M_⊕. However, M_Z,env* clearly depends on the value of Z₁. For Z₁ = 0.02, M_Z,env* is around 6 M_⊕, for Z₁ = 0.029, M_Z,env* is around 9 M_⊕, and for Z₁ = 0.035, M_Z,env* is around 11 M_⊕. 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_⊕.

Fig. 13 Zdilute vs. mdilute for two types of models. Zdilute is the maximum mass fraction of heavy elements in the dilute core region while mdilute controls the extent of the dilute core in terms of mass (see Sect. 2.2). Left panel: models using original EOSs. Z1 = 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 T1bar = 166.1 K and T1bar = 174.1 K. The thin gray dotted line corresponds to mdilute = 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).

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 Sect. 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 + δm_dil with δ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_⊕. 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_⊕ 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 elements can be added, which leads to these low values of m_dilute (and also M_Z,_dil* see Sect. 4.5).

5 Conclusion

We explored a wide variety of interior models of Jupiter constrained 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 composition, and an outer envelope of uniform Z composition. The helium phase separation is modelled as a jump in helium abundance 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 observe a range in the adiabatic density profiles of up to 5% at pressures ranging from 10 kbar to 10 Mbar. We interpret these variations as resulting from changes between different EOS tables and from different interpolation choices in relatively sparsely populated tables.

An important source of uncertainty results from our poor knowledge of Jupiter’s complex atmosphere and the possibility of a higher entropy than generally assumed. By allowing the T_{1 bar} parameter to vary (see Miguel et al. 2022), we obtain models that fit all constraints for all EOSs used. The values of T_{1 bar} obtained range from 171 to 188 K, significantly higher than 166.1 K from the Galileo probe (von Zahn et al. 1998), but within 164–174 K, 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_{1 bar} values, yield the highest values of J₆, which are slightly outside the range expected from differential rotation. Conversely, the MGF16+MH13 EOS leads to the smallest values of J₆, well within expectations, but the largest values of T_{1 bar}.

In all cases, we obtain a dilute core of heavy-element mass of between 10 and 25 Μ_⊕, 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 Μ_⊕. The total mass of heavy elements that we find ranges from 18 to 33 Μ_⊕. However, we must stress that, given the possibility of (perhaps significant) superadiabatic regions (see Guillot 1995; Leconte & Chabrier 2012; Leconte et al. 2017, and Sect. 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, dilute 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 Μ_⊕. This result is also promising, in light of interior models for Saturn, which indicate that a dilute core extends to 52–60 Μ_⊕ of the planet’s total mass (Mankovich & Fuller 2021). This could lead to a formation 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 solutions 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 constraints on R_eq and J₂ and/or modifications to the EOSs that were too important and most likely incompatible with the experimental 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 improvements 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. Improvements on EOSs, both experimentally and numerically, with particular emphasis on the hydrogen–helium mixture at pressures between 10 kbar and 10 Mbar and near Jupiter’s adi-abat (temperatures from 1000 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 efforts (Gonçalves et al. 2019; Shaw et al. 2022) would be an extremely 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.

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).

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. T1ba! is fixed at 166.1 K and Z1 = 0.02. The black points correspond to the measured J2n by Juno. The black error bars correspond to Juno’s measurements accounting for differential rotation for the J2n.

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_1bar = 166.1 K and respectively for Z₁ = 0.02 and Z₁ = 0.0286. Figures C.8 and C.9 show the posterior distributions of the MCMC simulations using modified EOS, with T_1bar = 174.1 K and respectively for Z₁ = 0.02 and Z₁ = 0.035.

Fig. C.1 Posterior distributions obtained with the MGF16+MH13 EOS, where T1bar is a free parameter, Z1 = 0.02 (1.3 × the protosolar value). The black points correspond to the measured J2n by Juno. The black error bars correspond to Juno’s measurements accounting for differential rotation for the J2n and Galileo’s measurement for T1bar.

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.

Fig. C.3 Same as Fig. C.1 but with the CD21 EOS.

Fig. C.4 Same as Fig. C.1 but with the HG23+CMS19 EOS.

Fig. C.5 Same as Fig. C.1 but with the HG23+MLS22 EOS.

Fig. C.6 Posterior distributions obtained with a modification of the EOS, where T1bar is fixed at 166.1 K, Z1 = 0.02 (1.3 × the protosolar value). The black points correspond to the measured J2n by Juno. The black error bars correspond to Juno’s measurements accounting for differential rotation for the J2n.

Fig. C.7 Same as Fig. C.6 but with Z1 = 0.0286 (1.9 × the protosolar value)

Fig. C.8 Posterior distributions obtained with a modification of the EOS, where T1bar is fixed at 174.1 K, Z1 = 0.02 (1.3 × the protosolar value). The black points correspond to the measured J2n by Juno. The black error bars correspond to Juno’s measurements accounting for differential rotation for the J2n.

Fig. C.9 Same as Fig. C.8 but with Z1 = 0.035 (2.3 × the protosolar value)

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₁ = 0.0153, T_1bar = 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_He, and comparable characteristics for the dilute core.

Fig. D.1 Posterior distributions obtained with MH13*. The red star shows the Militzer et al. (2022) preferred (static) model. Mcore = 0, T1bar is fixed at 166.1 K and Z1 = 0.0153. The black error bars correspond to Juno's measurements accounting for differential rotation for the J2n.

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.

References

All Tables

All Figures

Fig. 1 Abundances of heavy elements in the atmosphere of Jupiter. Here, ‘ices’ refers to the metallicity of the atmosphere considering only ices (CH4, NH3, H2O and H2S). Two values of Ztot 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.

In the text

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 < Z1, and the mass of heavy elements in the rest of the envelope (where Z < Z1). Z1 is the mass fraction of heavy elements in the outer envelope. Zdilute is the maximum mass fraction in the dilute core, mdilute monitors the extent of the dilute core in terms of mass.

In the text

	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).
In the text

Fig. 4 Comparison between the MH13* and MGF16+MH13 EOSs. Top panel: comparison of Tables 1 and 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 $\left(P,\rho /\sqrt{P}\right)$ points of Table 1 from Militzer & Hubbard (2013) are shown as blue crosses. The locations at which the entropy is equal to 7.078061 kbel.−1 (a typical value of entropy for Jupiter) according to the interpolations in the two tables are shown with different symbols, as labelled.

In the text

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 (Pmodif = 1011 dyn cm−2, ΔP = 0.5, dρ = −0.06 for the blue curve & Pmodif = 1011.5 dyn cnr−2, ΔP = 0.6, dρ = −0.05). Bottom panel: difference between the ΔU of the modified models and that of the reference model (MH13).

In the text

	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).
In the text

Fig. 7 Posterior distributions of Yproto, Yatm. T1bar, J4, and J6 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, Yproto is fixed. In the green run, Yproto and Yatm are fixed. In the red run, Yproto, Yatm, and T1bar 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 (Yproto = 0.277 ± 0.006, Yatm = 0.238 ± 0.005, T1bar = 165 ± 4K, J4.106 = −586.53 ± 0.0836 and J6.106 = 34.18 ± 0.07682). The gray areas correspond to the standard deviation of the prior (1 sigma).

In the text

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 ΔP (µ = 0.5 and σ = 0.02) and dρ (µ = 0. and σ = 0.01). Blue shows results obtained with uniform priors. The prior on Pmodif remains uniform in both cases. The dashed lines correspond to the adiabat obtained with the mean values of Pmodif, ΔP, and dρ (see Sect. 3.3) of a subsample of 100 models randomly drawn from the MCMC output. We compute the standard deviation (σ) of the 100 adiabats and the envelopes show the adiabats of the 1σ spread from the mean modified adiabat (dashed line). Here, T1bar is fixed at 166.1 K.

In the text

Fig. 9 Adiabats obtained for models with modified EOS. Top panel: T1bar is fixed at 166.1 K. The light blue area shows results for Z1 = 0.02 (1.3 × the protosolar value) while dark blue shows results for Z1 = 0.0286 (1.9 × the protosolar value). Bottom panel: T1bar is fixed at 174.1 K. The yellow area shows results for Z1 = 0.02 (1.3 × the protosolar value) while gray shows results for Z1 = 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.

In the text

Fig. 10 Temperature at 1 bar vs. helium transition pressures for two types of models. Left panel: models using original EOSs. Z1 = 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 T1bar = 166.1 and T1bar = 174.1 K. The black dotted line corresponds to the 1 bar temperature measured by the Galileo probe.

In the text

Fig. 11 Equatorial radius and gravitational moments for two types of models. Left panels: models using original EOSs. Z1 = 0.02. Right panels: models with a modification of the EOS. The inital EOS that has been modified is MH13. We present two subsets of models: with T1bar = 166.1 and T1bar = 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 Sect. 2.1).

In the text

	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: MZ,dil* vs. МZ,env* (see Sect. 2.2).
In the text

Fig. 13 Zdilute vs. mdilute for two types of models. Zdilute is the maximum mass fraction of heavy elements in the dilute core region while mdilute controls the extent of the dilute core in terms of mass (see Sect. 2.2). Left panel: models using original EOSs. Z1 = 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 T1bar = 166.1 K and T1bar = 174.1 K. The thin gray dotted line corresponds to mdilute = 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).

In the text

	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. T1ba! is fixed at 166.1 K and Z1 = 0.02. The black points correspond to the measured J2n by Juno. The black error bars correspond to Juno’s measurements accounting for differential rotation for the J2n.
In the text

	Fig. C.1 Posterior distributions obtained with the MGF16+MH13 EOS, where T1bar is a free parameter, Z1 = 0.02 (1.3 × the protosolar value). The black points correspond to the measured J2n by Juno. The black error bars correspond to Juno’s measurements accounting for differential rotation for the J2n and Galileo’s measurement for T1bar.
In the text

	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.
In the text

	Fig. C.3 Same as Fig. C.1 but with the CD21 EOS.
In the text

	Fig. C.4 Same as Fig. C.1 but with the HG23+CMS19 EOS.
In the text

	Fig. C.5 Same as Fig. C.1 but with the HG23+MLS22 EOS.
In the text

	Fig. C.6 Posterior distributions obtained with a modification of the EOS, where T1bar is fixed at 166.1 K, Z1 = 0.02 (1.3 × the protosolar value). The black points correspond to the measured J2n by Juno. The black error bars correspond to Juno’s measurements accounting for differential rotation for the J2n.
In the text

	Fig. C.7 Same as Fig. C.6 but with Z1 = 0.0286 (1.9 × the protosolar value)
In the text

	Fig. C.8 Posterior distributions obtained with a modification of the EOS, where T1bar is fixed at 174.1 K, Z1 = 0.02 (1.3 × the protosolar value). The black points correspond to the measured J2n by Juno. The black error bars correspond to Juno’s measurements accounting for differential rotation for the J2n.
In the text

	Fig. C.9 Same as Fig. C.8 but with Z1 = 0.035 (2.3 × the protosolar value)
In the text

	Fig. D.1 Posterior distributions obtained with MH13*. The red star shows the Militzer et al. (2022) preferred (static) model. Mcore = 0, T1bar is fixed at 166.1 K and Z1 = 0.0153. The black error bars correspond to Juno's measurements accounting for differential rotation for the J2n.
In the text