Issue 
A&A
Volume 662, June 2022



Article Number  A18  
Number of page(s)  16  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202243207  
Published online  08 June 2022 
Jupiter’s inhomogeneous envelope
^{1}
SRON Netherlands Institute for Space Research,
Niels Bohrweg 4,
2333 CA
Leiden,
The Netherlands
email: ymiguel@strw.leidenuniv.nl
^{2}
Leiden Observatory, University of Leiden,
Niels Bohrweg 2,
2333 CA
Leiden,
The Netherlands
^{3}
Heidelberg Institute for Theoretical Studies (HITS gGmbH),
SchlossWolfsbrunnenweg 35,
69118
Heidelberg,
Germany
^{4}
CITIES, NYUAD Institute, New York University Abu Dhabi,
PO Box 129188,
Abu Dhabi,
UAE
^{5}
Université Côte d'Azur, OCA, Lagrange CNRS,
06304
Nice,
France
^{6}
Department of Earth and Planetary Sciences, Weizmann Institute of Science,
Rehovot
76100,
Israel
^{7}
Lunar and Planetary Laboratory, University of Arizona,
Tucson,
AZ
85721,
USA
^{8}
Department of Earth and Planetary Science, University of California,
Berkeley,
CA
94720,
USA
^{9}
Institute for Computational Science, University of Zurich,
Winterthurerstr. 190,
8057
Zurich,
Switzerland
^{10}
University of Michigan, Climate and Space Sciences and Engineering,
Ann Arbor,
MI
48109,
USA
^{11}
Space Research Corporation,
Annapolis,
MD
21403,
USA
^{12}
NASA Goddard Space Flight Center,
Greenbelt,
MD
20771,
USA
^{13}
Department of Mechanical and Aerospace Engineering, Sapienza University of Rome,
Italy
^{14}
Department of Earth and Planetary Sciences, Harvard University,
Cambridge,
MA,
USA
^{15}
Department of Astronomy, Cornell University,
122 Sciences Drive,
Ithaca,
NY
14853,
USA
^{16}
Carl Sagan Institute, Cornell University,
122 Sciences Drive,
Ithaca,
NY
14853,
USA
^{17}
Division of Geological and Planetary Sciences, California Institute of Technology,
Pasadena,
CA
91125
USA
^{18}
Southwest Research Institute,
San Antonio,
TX
78238,
USA
Received:
27
January
2022
Accepted:
2
March
2022
Context. While Jupiter’s massive gas envelope consists mainly of hydrogen and helium, the key to understanding Jupiter’s formation and evolution lies in the distribution of the remaining (heavy) elements. Before the Juno mission, the lack of highprecision gravity harmonics precluded the use of statistical analyses in a robust determination of the heavyelement distribution in Jupiter’s envelope.
Aims. In this paper, we assemble the most comprehensive and diverse collection of Jupiter interior models to date and use it to study the distribution of heavy elements in the planet’s envelope.
Methods. We apply a Bayesian statistical approach to our interior model calculations, reproducing the Juno gravitational and atmospheric measurements and constraints from the deep zonal flows.
Results. Our results show that the gravity constraints lead to a deep entropy of Jupiter corresponding to a 1 bar temperature that is 515 K higher than traditionally assumed. We also find that uncertainties in the equation of state are crucial when determining the amount of heavy elements in Jupiter’s interior. Our models put an upper limit to the inner compact core of Jupiter of 7 M_{Earth}, independently of the structure model (with or without a dilute core) and the equation of state considered. Furthermore, we robustly demonstrate that Jupiter’s envelope is inhomogeneous, with a heavyelement enrichment in the interior relative to the outer envelope. This implies that heavyelement enrichment continued through the gas accretion phase, with important implications for the formation of giant planets in our Solar System and beyond.
Key words: planets and satellites: interiors / planets and satellites: gaseous planets / planets and satellites: formation / planets and satellites: composition
© ESO 2022
1 Introduction
The standard model of Jupiter formation starts with the accretion of solids followed by a rapid gas accretion phase wherein Jupiter’s hydrogen and helium envelope is captured from the primitive solar nebula. In this scenario, the dominant size of the accreted solids could be either kilometresized planetesimals (Pollack et al. 1996) or centimetresized pebbles (Lambrechts et al. 2014), with dramatic implications for formation timescales and the distribution of heavy elements in Jupiter’s envelope (Venturini & Helled 2020; Vazan et al. 2018). In the planetesimaldriven scenario, the accretion of solids continues through the gas accretion phase and stops when all of the planetesimals in the planet’s vicinity have been accreted. The fragmentation and ablation of these solid planetesimals cause a nonhomogenous distribution of heavy elements in the envelope (Alibert et al. 2018). In contrast, in the pebbledriven scenario, the fast orbital decay of pebbles caused by gas drag provides a continuous resupply of solid material that enriches the growing planet (Ormel et al. 2021). Nevertheless, the supply stops once the socalled pebbleisolation mass is reached (Lambrechts et al. 2014), after which only gas accretion continues unless or until pebbles grow to planetesimal size. Here we seek to distinguish between these scenarios, analysing the internal structure of Jupiter today and determining whether Jupiter’s envelope harbours a nonhomogeneous distribution of heavy elements. We explore the possibility of differentiation of the heavy elements in the envelope, testing whether this differentiation is a natural outcome of the set of models that reproduce all observational constraints, including those provided by the recent Juno observations.
2 Methods
2.1 Jupiter’s structure: Threelayer and dilutecore models
We constructed two contrasting sets of Jupiter interior models and explored the largest possible ensemble of realistic models to date that represent the interior of Jupiter. The first set of models was constructed using a threelayer Jupiter model, where Jupiter’s interior consists of an outer H_{2}dominated region, an intermediate H_{metallic}dominated region (e.g. Hubbard & Militzer 2016), and a core made of 100% heavy elements (Fig. 1a). The second set of models was built using three layers and adding a dilute core: a region above the inner core where the H and He of the envelope are gradually mixed with a greater abundance of heavy elements (Wahl et al. 2017; see our Fig. 1b). For the latter models, we characterised the dilute core with three parameters: the maximum mass mixing ratio of heavy elements in the diluted core region (Z_{dilute}), a parameter that controls the extent of the dilute core in terms of mass (m_{dilute}), and a parameter controlling how steep the heavyelement gradient (δm_{dil}) is. In the end, the ratio of the heavyelement mass mixing between the inner core outer boundary and the helium transition is given by (1)
with δm_{dil} = 0.075. We tested possible variations of this parameter and found that the results of the calculation are not sensitive to the choice of δm_{dil}. In addition, for all the models, we allowed for an increase in the interior entropy in the heliumdemixing region by including a temperature jump (ΔT_{He}), following Hubbard & Militzer (2016). An a posteriori examination of the models shows that this region is convectively stable for values of this parameter lower than about 2000 K (Guillot et al. 2018). The outermost layer of the envelope in both sets of models is separated from the inner layer of the envelope by the transition to a region of immiscibility of helium in hydrogen ("helium rain"). The precise location of this region is not well known, but numerical calculations (Morales et al. 2013; Schöttler et al. 2018) suggest that the immiscibility should be located between ~0.8 and ~3 Mbar, and lab experiments indicate that it may extend even deeper into the planet (Brygoo et al. 2021). In our models, we assume a Gaussian distribution, taking these constraints into account. We ran Markov chain Monte Carlo (MCMC) calculations for the two scenarios using a Bayesian statistical model (Bazot et al. 2012), wherein calculations are concentrated around the relevant regions of the parameter space. This allows millions of solutions (models that fit Jupiter’s observational constraints) to be evaluated with a feasible computational cost. More details on the Bayesian approach are provided in the appendix.
2.2 Details on each Jupiter interior model
We calculated each Jupiter interior structure model with the code CEPAM (Guillot & Morel 1995; Guillot et al. 2018), which solves the equations of hydrostatic equilibrium, mass and energy conservation, and energy transport. We performed multiple sets of MCMC runs using different equations of state. Because Jupiter is primarily composed of hydrogen and helium, the equations of state of these two elements are fundamental when modelling its interior structure (Miguel et al. 2016). Several equations of state have been published in the past few years, and many uncertainties still remain (Mazevet et al. 2022). Therefore, we repeated our calculations using several proposed equations of state to avoid a bias in our results. We used MH13H (Militzer & Hubbard 2013; Miguel et al. 2016), CMS19H (Chabrier et al. 2019), and MLS21H (Mazevet et al. 2022) for the equation of state of hydrogen, and for helium we used SCVH95He (Saumon et al. 1995) and CMS19He (Chabrier et al. 2019). For the heavy elements, we used the equations of state of a mixture of silicates (dry sand) and the one of water in SESAME (Lyon & Johnson 1992).
Jupiter’s adiabat is prescribed by its entropy or, equivalently, by its temperature at the 1 bar level (T_{1bar}). Voyager radio occultations (Lindal 1981) initially indicated T_{1bar} = 165 ± 5 K. In situ measurements from the Galileo probe measured, in a socalled hot spot, a 1 bar temperature of 166.1 K (Seiff et al. 1998). However, a reassessment of the Voyager radio occultations after the corrected helium abundance as measured by the Galileo probe led to temperatures at the occultation latitudes that are 35 K higher (Gupta et al. 2022), raising the possibility that the Galileo probe temperature measurement may not be strictly appropriate to the entire planet. We used a Gaussian distribution for T_{1bar} centred at 165 K and with a dispersion in agreement with these observations. Measurements by the Galileo probe and recent results by the Juno mission provide constraints for the metallicity and helium abundance in a shallow region within Jupiter’s H_{2}dominated layer. We took these constraints into account and used a He abundance of Y_{1} = 0.238 (von Zahn et al. 1998) and metallicity Ζ_{1} = 1 Z_{solar} in our calculations (Li et al. 2020). For the H_{metalic} dominated region, we required that the overall abundance of He be consistent with the protosolar abundance (Y_{proto} = 0.277, Serenelli et al. 2010) and allowed either different metallicity from the H_{2}dominated region or the same (the latter option for the dilutecore models).
Fig. 1 Schematic view of the Jupiter structure models used in this paper. Panel a: threelayer model. Panel b: dilutecore model. The different regions are indicated with arrows. 
Fig. 2 Offsets between ToF and CMS for two subsamples (MCMC models and optimised models). Left panels: offsets on the gravitational moments. The dashed black line shows the linear fit of our models. The green dot shows the origin. The yellow star corresponds to the linear fit of the optimised model median of ΔJ_{2}. The pink error bars show the uncertainty of the Juno measurements that account for the differential rotation. The purple error bars show the uncertainty of the Juno measurements. The green cross shows the former offsets from Guillot et al. (2018). Right panels: residuals. 
2.3 Gravity harmonics calculation
The modelled density profiles were used to calculate the even gravity harmonics (J_{2n},n=1,…,5), which were ultimately compared with the Junoderived values to find the models that are consistent with the Juno observations. The even gravity harmonics measured by Juno include contributions from both the static background state (calculated with our interior models that assume rigid body rotation) and dynamical processes (mostly from the deep atmospheric zonal winds; Kaspi et al. 2017, 2018; Kulowski et al. 2020). We can calculate the even gravity harmonics produced from the interior density profile by subtracting the computed dynamical contributions from the Juno measurements. Therefore, the effective gravity harmonics that our interior model calculations must match are (n = 1,…_{;} 5), where are the values provided by Juno mission measurements (Iess et al. 2018; Durante et al. 2020) and is the contribution to the gravity harmonics due to differential rotation. We calculated the static component of the gravity harmonics using the theory of figures (ToF) of fourth order (Zharkov & Trubitsyn 1978; Nettelmann 2017) combined with an integration of the reconstructed density structure in two dimensions using a GaussLegendre quadrature (Guillot et al. 2018) or using the concentric MacLaurin spheroid (CMS) theory (Hubbard 2013). Because results calculated with the CMS method are more accurate but more computationally demanding, in this paper the majority of the runs are made using the first method, which allows us to perform many more runs in a shorter time; however, we use the accuracy of the CMS method by calibrating the gravity harmonics obtained from the ToF to the CMS values (Guillot et al. 2018). Because the gravity harmonics measured by Juno have reached a very high accuracy, estimating these offsets is essential.
This calibration was done by assessing offsets (δJ_{2n} = J_{2n}(ToF)  J_{2n}(CMS)), between the gravity harmonics from the ToF and the ones from the CMS method, using a random sample of 100 of our preferred models. We then took these models and performed an optimisation procedure, modifying M_{core} and Z_{1} to perfectly fit the observational measurements of the equatorial radius R_{eq} and J_{2}  to get offsets from the most accurate models  and then computed the offsets for this second set of models using both methods. Figure 2 displays the offsets for both subsamples (models with and without the optimisation procedure). Our results show a correlation between the offsets, where δJ_{4} and δJ_{6} depend strongly on δJ_{2}. Higher order gravitational moments (δJ_{8}, δJ_{10}) are less dependent on δJ_{2}. Thanks to these linear relationships, we can quantify the impact of an offset on J_{2} on the offsets of the higher order gravitational moments. We also note that the residuals are very small compared to the Juno error bars that consider differential rotation (Fig. 3). We calculated the median value of δJ_{2} among the 100 optimised models and calculated the higher order offsets using the linear relationships. Table 1 lists the offsets found; we note that previous offsets from Guillot et al. (2018) lie almost on the linear regression curves, showing that our offsets have changed very slightly and that our calculations are thus robust.
2.4 Dynamical contribution to the gravity harmonics
In order to provide the plausible range of dynamical contribution to the even gravity harmonics , we considered the widest possible range of flow profiles that still matched the gravity harmonics. The Juno gravity measurements revealed significant northsouth asymmetries in Jupiter’s gravity field (Iess et al. 2018). Such hemispherical asymmetries (represented through the measurement of J_{3}, J_{5}, J_{7}, and J_{9}) can only be due to interior flows. The resemblance of the measured gravity signature of the flow to that obtained by extending the eastwest (zonal) observed cloudlevel winds inwards (Kaspi 2013) suggested that not only do the flows reach deep (~3000 km, where the pressure is ~100 kbar), but that the flow pattern is generally similar to the flow profile at the cloud level (Kaspi et al. 2018, 2020). Although variations from this profile at depth are possible, it is statistically unlikely that the deep flow profile varies significantly, and it will still match all four measured odd gravity harmonics (Tollefson et al. 2017; Duer et al. 2020). Any deep flow profile will also have a signature in the even gravity harmonics, providing a dynamical contribution to the measurements beyond that coming from the interior density profile.
To determine such a dynamical contribution to the even gravity harmonics, we retrieved the wind field from the gravity field data using an adjoint inverse technique (Galanti et al. 2016, 2017). We allowed the zonal wind structure to vary both in its latitudinal profile (from the observed cloudlevel wind structure) and its vertical profile (Galanti et al. 2019). This was achieved through randomly varying , n=1,5 within a range of 30σ of the Juno values (we chose a large range because it is not known what fraction of the even harmonic measurements comes from the dynamics; Kaspi et al. 2017). Every profile was extended inwards along the direction of the spin axis, assuming a radial decay profile; then, given that the flow must be in thermal wind balance with the density field, the balancing density field was integrated to give the dynamical gravity harmonics (Kaspi et al. 2016). For each of the 3000 cases, we calculated the best fitting wind profile that is also consistent with the odd gravity harmonics, thus matching seven gravity values for each case (odd harmonics to and even differential contributions # to ). Once all solutions were obtained (grey profiles in Figs. 4a and b for the meridional and vertical profiles, respectively), we selected as plausible solutions those where the meridional profile of the zonal wind was within 20 m s^{−1} of the observed cloudlevel profile (black profiles in Figs. 4a, b). This uncertainty is consistent with the measurement error (Kaspi et al. 2020). The resulting range of the physically plausible dynamical contribution to the even gravity harmonics is: , , , , and ) (where ).
Fig. 3 Results for four sets of 100 models: random models from our preferred MCMC runs, optimised models with new offsets, optimised models with former offsets from Guillot et al. (2018), and optimised models with CMS calculation. Top left panel: equatorial radius vs. J_{2}. Top right panel: core mass vs. fraction of ices in the molecular hydrogen layer. Grey lines show the pairing between each model and its optimised version. Middle and bottom panels: gravitational moments. The black error bars show the uncertainty of the Juno measurements that account for the differential rotation. 
Offsets between ToF and CMS.
Fig. 4 Contribution to the gravity harmonics from the deep atmospheric flow. Shown are the solutions for the cloudlevel wind (a), the radial decay profile (b), and the windinduced even gravity harmonics (cf). The grey lines and dots show all of the 3000 cases, and the black lines and dots show the plausible solutions. The red line in (a) shows the observed cloudlevel wind, and the green line in (b) shows the solution in Kaspi et al. (2018). The red crosses in (ef) show the expected mean value of the differential contribution. 
3 Results
3.1 Envelope inhomogeneity
All of our models reproduce J_{2} and Jupiter’s radius, and Figs. 5 (panels ac) and 6 (panels ac) show the values of J_{4} and J_{6} for the MCMC solutions in the case of a threelayer model and with a dilutecore model, respectively (the other J_{2n} are shown in Methods). We also note that all our models are consistent with the metallicities observed by Juno and the helium abundance in Jupiter’s atmosphere as measured by the Galileo probe. When analysing the distribution of heavy elements in the envelope, we considered the envelope to consist of the H_{2} and H_{metallic} dominated layers for the threelayer models and to consist of the H_{2}  and H_{metallic}dominated layers and the dilute core for the dilutecore models. In Fig. 5d we show results of the threelayer models, and we see that the difference in heavy elements between the H_{2} and H_{metallic}dominated layers (ΔΖ = Z_{2}  Z_{1}) is larger than zero for all models calculated with the MH13 and MLS20 equations of state. For the CMS19 equation of state, we find that the probability of finding models with ΔΖ > 0 is 97.6%.
Figure 6d shows the difference in heavy elements within the envelope for models with a dilute core. In this case, we calculated the difference in heavy elements between the H2dominated and dilutecore regions (ΔΖ = Z_{dilute}  Z_{1}). We see that all models have ΔZ > 0, independently of the equation of state used in the calculations. Our results robustly demonstrate that Jupiter’s envelope is not homogeneous: the external layer of the envelope is depleted of heavy elements compared to the inner parts of the envelope. This result is independent of both the models adopted for the interior of the planet and the equation of state used. We note, however, that different equations of state lead to different distributions of the heavy elements in the interior of the planet.
3.2 Distribution and mass of the heavy elements
An analysis of the mass of heavy elements in the different layers for a random sample of 1000 of our models is shown in Fig. 7. We see that the total mass of heavy elements varies between 11 and 30 M_{Earth}, with differences resulting from the choice of the equation of state. Calculations done with the MH13 equation of state have , models with MLS21 have , and models with CMS19 have , independently of the model of Jupiter adopted (three layers or dilute core). For the threelayer models, the differences mostly arise from discrepancies in the mass of heavy elements in the H_{metallic}dominated region , which varies between 2 and 23 M_{Earth} depending on the equation of state. For models with a dilute core, the differences are mostly due to differences in the mass of heavy elements in the dilutecore region , which is found to vary between 1 and 25 M_{Earth} depending on the equation of state adopted in the calculation. The mass of heavy elements in the H_{2}dominated region is quite similar for models with different equations of state, independently of the model adopted for Jupiter. This is expected given the prior used to match the observational constraints from the Juno and Galileo missions. Regarding the inner core, we find that it varies between 0 and 7 M_{Earth} for all models, with its exact value depending on the model adopted for Jupiter’s interior. The threelayer models have M_{core}  6 M_{Earth}, independently of the equation of state adopted. Conversely, we see two groups of solutions for models with a dilute core: a group with inner core masses between 3 and 6 M_{Earth} and another group with small masses of up to 2 M_{Earth}.
Interestingly, models tend to have temperatures at 1 bar going from ~169 to ~188 K for the threelayer models to values of between ~165 and ~182K for the dilutecore models, depending on the equation of state adopted (see Sect. 4). While these values are high compared to observations by the Galileo probe (166.1 K; Seiff et al. 1998), it is not clear whether those in situ measurements represent the typical 1 bar temperature on Jupiter because of latitudinal variability, which may exist to a limited extent (Fletcher et al. 2020). Furthermore, a reassessment of the Galileo probe data led to an increase in the temperature of ≃4K (Gupta et al. 2022), and, additionally, the possibility of superadiabaticity in the interior (Leconte et al. 2017; Guillot et al. 2020) could yield a deep entropy corresponding to a temperature a few degrees higher than the measured value at 1 bar. In all cases, dilutecore models have temperatures more in agreement with the expectations given the values observed and the uncertainties in this parameter. Regarding the separation pressure for the immiscibility of helium in hydrogen, our values are always close to 3 Mbar (see the appendix), which is in agreement with the higher limit of numerical calculations (Morales et al. 2013; Schöttler et al. 2018) and more in agreement with recent laboratory experiments estimations (Brygoo et al. 2021).
Fig. 5 Details of models with three layers. Panels a, b, and c: J_{4} vs. J_{6} for threelayer models calculated using different equations of state for H, indicated in the figure. Cyan squares show the Juno measurements , and gravitational harmonics corrected by differential rotation  the ones to match with our interior models  are shown in red. Panel d: histogram showing the difference in heavy elements (ΔZ = Z_{2}  Z_{1}) between the H_{metallic} and H_{2}dominated layers. 
Fig. 6 Specifics of models with dilute cores. Panels a, b, and c: J_{4} vs. J_{6} obtained considering different equations of state for H. The Juno measurement is indicated with the cyan square. The effective gravity harmonics to be matched by our interior models are shown in red. Panel d: histogram that shows the difference in heavy elements (ΔZ = Z_{dilute}  Z_{1}) between the dilutecore region and the H_{2}dominated layer. 
3.3 Comparison with previous models
In this paper we perform a large exploration of the parameter space. We find that our results include most other individual models presented in other works.
Dilutecore models were first modelled by Wahl et al. (2017), who, using the MH13 equation of state, found a value of M_{core} + M_{dilute} between 10 and 24 M_{Earth} and of 2427 M_{Earth}. These results are consistent with our results, even though most of them have subsolar atmospheric metallicities.
Ni (2019) performed fourlayer models using the MH13 and CMS19 equations of state. They found a range of M_{core}+M_{dilcore} of 6.527 M_{Earth} and of 2428 M_{Earth} with the MH13 equation of state and M_{core} + M_{core} of 312 M_{Earth} and of 812 M_{Earth} with the CMS19 equation of state. For both equations of state, the results are in good qualitative agreement with ours, even though, in this case again, their optimisation led to subsolar atmospheric metallicities as opposed to our results, where Z_{1} = 1 Z_{solar}.
Debras & Chabrier (2019) imposed the constraint of a minimum atmospheric metallicity Z = 0.02. They found that in order to fit this and all other constraints using the CMS19 equation of state, an inward decrease of the abundance of heavy elements between the H_{2} and H_{metallic}dominated regions was required. We also find these types of solutions when using the CMS19 equation of state and threelayer models (Fig. 5), but, since they represent only 2.4% of our sample, we find it to be unlikely. The found by Debras & Chabrier (2019) are 2530 and 4045 M_{Earth} for models without and with inner cores, respectively. Our results with CMS19 always led to considerably smaller values. A similar discrepancy with the Debras & Chabrier (2019) study is found by Nettelmann et al. (2021).
Nettelmann et al. (2021) used CMS19 and dilutecore models. They found M_{core} of up to 3.8 M_{Earth} and up to 13 M_{Earth}, comparable to our results. We note that in Nettelmann et al. (2021) the authors calculate the gravitational harmonics using an expansion of the ToF of seventh order, different from the method used here (Sect. 2.4). Similar to our conclusions, they also find that higher temperatures at 1 bar help in reaching higher metallicities in the atmosphere. Their models have separation pressures between the H_{2} and H_{metallic}dominated regions close to 6 Mbar, while our models use separation pressures of around 3 Mbar.
Fig. 7 Mass of heavy elements in the different layers as a function of the total mass of heavy elements in Jupiter for a random sample of 1000 models extracted from our threelayer models (a, b, and c) and from models with a dilute core (d, e, f, and g). Different colours show models calculated with different equations of state. Panels a and d show the mass in the H_{2}dominated region in the yaxis, panels b and e show the mass in the H_{metallic}dominated region, panel f shows the mass in the dilute core, and panels c and g show the mass in the inner core. 
4 Conclusions
This study comprehensively reproduces observational constraints from the Juno measurements (the even and odd gravity harmonics and water abundance in the atmosphere), along with helium measurements from the Galileo probe, exploring different models for Jupiter’s interior and considering all recent equations of state. We show that the gravity constraints point to a deep entropy of Jupiter that corresponds to a 1 bar temperature that is higher than traditionally assumed (i.e. 170180 K rather than 166K). We robustly demonstrate that the heavyelement abundance is not homogeneous in Jupiter’s envelope. Our results imply that Jupiter continued to accrete heavy elements in large amounts while its hydrogenhelium envelope was growing, contrary to predictions based on the pebbleisolation mass in its simplest incarnation (Lambrechts et al. 2014), favouring instead planetesimalbased or more complex hybrid models (Alibert et al. 2018; Guilera et al. 2020). Furthermore, the envelope did not mix completely during the planet’s subsequent evolution, not even when Jupiter was young and hot (Vazan et al. 2018). Our result clearly shows the need for further exploration of nonadiabatic interior models for the giant planets, and it provides a base example for exoplanets: a nonhomogeneous envelope implies that the metallicity observed is a lower limit to the planet bulk metallicity. Therefore, metallicities inferred from remote atmospheric observations in exoplanets might not represent the bulk metallicity of the planet. Moreover, we demonstrate that knowledge of the equation of state is crucial in determining the mass of heavy elements in the interior of Jupiter, and we put important constraints on Jupiter’s inner core, which is found to be up to 7 M_{Earth}, a result that is independent of the interior model and equation of state adopted in the calculations.
Acknowledgments
The authors would like to thank the Juno Interior Working Group for useful discussions. We also thank S. Mazevet for providing the equation of state used in his paper (Mazevet et al. 2022) and for useful discussions about this topic. JIL was supported by the Juno Project through a subcontract with the Southwest Research Institute.
Appendix A Bayesian statistical model
Taking into account the uncertainties in the equation of state, the temperature at 1 bar, the pressure where heliumrain happens, and the error bars on the measurements themselves, we employed a Bayesian statistical approach to explore a large ensemble of Jupiter models and find which are the possible scenarios for Jupiter’s internal structure that fit all observational constraints. The parameters of CEPAM we allowed to vary are (A.1)
and the data we allowed to vary are (A.2)
The quantity of interest is the posterior density function of the parameters that are conditional on the data, p(θX). Bayes’ theorem tells us that the posterior density function is completely specified by the likelihood, p(Χθ), and the prior density function, p(θ), through the relation (A.3)
The likelihood is the probability density of the data, which is conditional on the parameters. In general, it is not a probability density for the parameter. The prior density encodes the information we have on the parameters before the inference processes (e.g. theoretical limits and previous measurements).
It is possible, under some assumptions, to write analytical forms for the likelihood and the prior. To define the former, we first assumed an additive statistical model, (A.4)
Here J is a mapping from the space of parameters to the space of observables. In practice, it represents our code solving the equations of planetary structure. The vector є represent a random noise. The distribution of this noise is thus the distribution of X _{J}(θ). We assume now that the distribution of є is normal and that its covariance matrix is diagonal. The likelihood function is therefore (A.5)
The represent the diagonal elements of the covariance matrix and are given by the observational uncertainties corresponding to the components of X. The adopted values are given in Tables A.1 and A.2.
In order to simplify the prior density, we assumed that all parameters are independent. It can thus be written as the product of univariate densities, (A.6)
The chosen priors are given in Tables A.1 and A.2. The boundaries of all priors were chosen using a testandtrial stage, ensuring that they were wide enough to avoid numerical issues during the sampling stage and that they were narrow enough so that the parameter space to sample is not too vast.
We used an affine invariant MCMC algorithm (Christen & Fox 2010; Goodman & Weare 2010) to sample the space of parameters and to approximate the posterior probability functions. It uses the relative positions in the parameter space of several Markov chains, run in parallel, to efficiently adapt the proposition law. Using 512 walkers, we found that good convergence for the MCMC simulations can be reached with approximately 10000 iterations of the algorithm.
Figures A.1, A.2, A.3, A.4, A.5, and A.6 show the corner plots with all the parameters that result from the MCMC calculations for the six cases considered (see the main text for an extended discussion on the results).
Parameters explored in our MCMC calculations for threelayer models. 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. Y_{proto}=0.277, Y_{1}=0.238, and Z_{1} =0.0153.
Parameters explored in our MCMC calculations for dilutecore models. 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. Y_{proto} = 0.277, , Z_{1} = 0.0153, and Z_{2} = 0.0153.
Fig. A.1 MCMC corner plot showing the posterior distribution of variables obtained with the threelayer model and the MH13 equation of state. Red points with error bars show the observed parameters as a reference. 
Fig. A.2 Posterior distribution results of the MCMC runs with three layers and the CMS 19 equation of state. Observed parameters are indicated with the red points. 
Fig. A.3 Distribution results of the MCMC runs with three layers and the MLS20 equation of state. 
Fig. A.4 Corner plot resulting from the runs with the dilutecore model and the MH13 equation of state. Red dots show the observed parameters. 
Fig. A.5 Distribution of different parameters resulting from the runs with the dilutecore model and the CMS 19 equation of state. 
Fig. A.6 Resulting distributions obtained with the dilutecore model and the MLS20 equation of state. 
References
 Alibert, Y., Venturini, J., Helled, R., et al. 2018, Nat. Astron., 2, 873 [NASA ADS] [CrossRef] [Google Scholar]
 Bazot, M., Bourguignon, S., & ChristensenDalsgaard, J. 2012, MNRAS, 427, 1847 [NASA ADS] [CrossRef] [Google Scholar]
 Brygoo, S., Loubeyre, P., Millot, M., et al. 2021, Nature, 593, 517 [NASA ADS] [CrossRef] [Google Scholar]
 Chabrier, G., Mazevet, S., & Soubiran, F. 2019, ApJ, 872, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Christen, J. A., & Fox, C. 2010, Bayesian Anal., 5, 263 [CrossRef] [Google Scholar]
 Debras, F., & Chabrier, G. 2019, ApJ, 872, 100 [Google Scholar]
 Duer, K., Galanti, E., & Kaspi, Y. 2020, J. Geophys. Res. (Planets), 125, 8 [Google Scholar]
 Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e86572 [NASA ADS] [CrossRef] [Google Scholar]
 Fletcher, L. N., Orton, G. S., Greathouse, T. K., et al. 2020, J. Geophys. Res., 125, 8 [Google Scholar]
 Galanti, E., & Kaspi, Y. 2016, ApJ, 820, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Galanti, E., & Kaspi, Y. 2017, Icarus, 286, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Galanti, E., Kaspi, Y., Miguel, Y., et al. 2019, Geophys. Res. Lett., 46, 2 [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
 Guilera, O. M., Sándor, Z., Ronco, M. P., et al. 2020, A&A, 642, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guillot, T., & Morel, P. 1995, A&AS, 109, 109 [NASA ADS] [Google Scholar]
 Guillot, T., Miguel, Y., Militzer, B., et al. 2018, Nature, 555, 227 [Google Scholar]
 Guillot, T., Li, C., Bolton, S., et al. 2020, J. Geophys. Res., 125, 8 [Google Scholar]
 Gupta, P., Atreya, S., Steffes, P. G. et al. 2022, PSJ, submitted [Google Scholar]
 Hubbard, W. B. 2013, ApJ, 768, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Hubbard, W. B., & Militzer, B. 2016, ApJ, 820, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Iess, L., Folkner, W. M., Durante, D., et al. 2018, Nature, 555, 220 [NASA ADS] [CrossRef] [Google Scholar]
 Kaspi, Y. 2013, Geophys. Res. Lett., 40, 676 [NASA ADS] [CrossRef] [Google Scholar]
 Kaspi, Y., Davighi, J. E., Galanti, E., & Hubbard, W. B. 2016, Icarus, 276, 170 [NASA ADS] [CrossRef] [Google Scholar]
 Kaspi, Y., Guillot, G., Galanti, E., et al. 2017, Geophys. Res. Lett., 44, 12 [Google Scholar]
 Kaspi, Y., Galanti, E., Hubbard, W. B., et al. 2018, Nature, 555, 223 [Google Scholar]
 Kaspi, Y., Galanti, E., Showman, A. P., et al. 2020, Space Sci. Rev., 216, 84 [Google Scholar]
 Kulowski, L., Cao, H., & Bloxham, J. 2020, J. Geophys. Res.: Planets, 125 [Google Scholar]
 Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Leconte, J., Selsis, F., Hersant, F., & Guillot, T. 2017, A&A, 598, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Li, C., Ingersoll, A., Bolton, S., et al. 2020, Nat. Astron., 4, 609 [Google Scholar]
 Lindal, G. F., et al. 1981, J. Geophys. Res., 86, 8721 [NASA ADS] [CrossRef] [Google Scholar]
 Lyon, S., & Johnson, J. 1992, LANL Report LAUR923407, Los Alamos [Google Scholar]
 Mazevet, S., Licari, A., & Soubiran, F. 2022, A&A, in press, https://doi.org/10.1051/00046361/201935764 [Google Scholar]
 Miguel, Y., Guillot, T., Fayon, L. 2016, A&A, 596, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Militzer, B., & Hubbard, W. B. 2013, ApJ, 774, 148 [NASA ADS] [CrossRef] [Google Scholar]
 Morales, M. A., Hamel, S., Caspersen, K., & Schwegler, E. 2013, Phys. Rev. B, 87, 174105 [NASA ADS] [CrossRef] [Google Scholar]
 Nettelmann, N. 2017, A&A, 606, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Nettelmann, N., Movshovitz, N., Ni, D., et al. 2021, Planet. Sci. J., 2, 241 [NASA ADS] [CrossRef] [Google Scholar]
 Ni, D. 2019, A&A, 632, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ormel, C. W., Vazan, A., & Brouwers, M. G. 2021, A&A, 647, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
 Seiff, A., Kirk, D. B., Knight, T. C. D., et al. 1998, J. Geophys. Res., 103, 22857 [NASA ADS] [CrossRef] [Google Scholar]
 Serenelli, A. M., & Basu, S. 2010, ApJ, 719, 865 [NASA ADS] [CrossRef] [Google Scholar]
 Schöttler, M., Redmer, R. 2018, Phys. Rev. Lett., 120 [Google Scholar]
 Tollefson, J., Wong, M. H., Pater, J., et al. 2017, Icarus, 296, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Valletta, C., & Helled, R. 2019, ApJ, 871, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Vazan, A., Helled, R., & Guillot, 2018, A&A, 610, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Venturini, J., & Helled, R. 2020, A&A, 634, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 von Zahn, U., Hunten, D. M., & Lehmacher, G. 1998, J. Geophys. Res., 103, 22815 [NASA ADS] [CrossRef] [Google Scholar]
 Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649 [CrossRef] [Google Scholar]
 Zharkov, V. N., & Trubitsyn, V. P. 1978, Phys. Planet. Inter. Astron. Astrophys. Ser., Tucson: Pachart [Google Scholar]
All Tables
Parameters explored in our MCMC calculations for threelayer models. 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. Y_{proto}=0.277, Y_{1}=0.238, and Z_{1} =0.0153.
Parameters explored in our MCMC calculations for dilutecore models. 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. Y_{proto} = 0.277, , Z_{1} = 0.0153, and Z_{2} = 0.0153.
All Figures
Fig. 1 Schematic view of the Jupiter structure models used in this paper. Panel a: threelayer model. Panel b: dilutecore model. The different regions are indicated with arrows. 

In the text 
Fig. 2 Offsets between ToF and CMS for two subsamples (MCMC models and optimised models). Left panels: offsets on the gravitational moments. The dashed black line shows the linear fit of our models. The green dot shows the origin. The yellow star corresponds to the linear fit of the optimised model median of ΔJ_{2}. The pink error bars show the uncertainty of the Juno measurements that account for the differential rotation. The purple error bars show the uncertainty of the Juno measurements. The green cross shows the former offsets from Guillot et al. (2018). Right panels: residuals. 

In the text 
Fig. 3 Results for four sets of 100 models: random models from our preferred MCMC runs, optimised models with new offsets, optimised models with former offsets from Guillot et al. (2018), and optimised models with CMS calculation. Top left panel: equatorial radius vs. J_{2}. Top right panel: core mass vs. fraction of ices in the molecular hydrogen layer. Grey lines show the pairing between each model and its optimised version. Middle and bottom panels: gravitational moments. The black error bars show the uncertainty of the Juno measurements that account for the differential rotation. 

In the text 
Fig. 4 Contribution to the gravity harmonics from the deep atmospheric flow. Shown are the solutions for the cloudlevel wind (a), the radial decay profile (b), and the windinduced even gravity harmonics (cf). The grey lines and dots show all of the 3000 cases, and the black lines and dots show the plausible solutions. The red line in (a) shows the observed cloudlevel wind, and the green line in (b) shows the solution in Kaspi et al. (2018). The red crosses in (ef) show the expected mean value of the differential contribution. 

In the text 
Fig. 5 Details of models with three layers. Panels a, b, and c: J_{4} vs. J_{6} for threelayer models calculated using different equations of state for H, indicated in the figure. Cyan squares show the Juno measurements , and gravitational harmonics corrected by differential rotation  the ones to match with our interior models  are shown in red. Panel d: histogram showing the difference in heavy elements (ΔZ = Z_{2}  Z_{1}) between the H_{metallic} and H_{2}dominated layers. 

In the text 
Fig. 6 Specifics of models with dilute cores. Panels a, b, and c: J_{4} vs. J_{6} obtained considering different equations of state for H. The Juno measurement is indicated with the cyan square. The effective gravity harmonics to be matched by our interior models are shown in red. Panel d: histogram that shows the difference in heavy elements (ΔZ = Z_{dilute}  Z_{1}) between the dilutecore region and the H_{2}dominated layer. 

In the text 
Fig. 7 Mass of heavy elements in the different layers as a function of the total mass of heavy elements in Jupiter for a random sample of 1000 models extracted from our threelayer models (a, b, and c) and from models with a dilute core (d, e, f, and g). Different colours show models calculated with different equations of state. Panels a and d show the mass in the H_{2}dominated region in the yaxis, panels b and e show the mass in the H_{metallic}dominated region, panel f shows the mass in the dilute core, and panels c and g show the mass in the inner core. 

In the text 
Fig. A.1 MCMC corner plot showing the posterior distribution of variables obtained with the threelayer model and the MH13 equation of state. Red points with error bars show the observed parameters as a reference. 

In the text 
Fig. A.2 Posterior distribution results of the MCMC runs with three layers and the CMS 19 equation of state. Observed parameters are indicated with the red points. 

In the text 
Fig. A.3 Distribution results of the MCMC runs with three layers and the MLS20 equation of state. 

In the text 
Fig. A.4 Corner plot resulting from the runs with the dilutecore model and the MH13 equation of state. Red dots show the observed parameters. 

In the text 
Fig. A.5 Distribution of different parameters resulting from the runs with the dilutecore model and the CMS 19 equation of state. 

In the text 
Fig. A.6 Resulting distributions obtained with the dilutecore model and the MLS20 equation of state. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.