EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 662 (June 2022) A&A, 662 (2022) A18 Full HTML
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/0004-6361/202243207
Published online 08 June 2022

© 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 kilometre-sized plan-etesimals (Pollack et al. 1996) or centimetre-sized 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 planetesimal-driven 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 non-homogenous distribution of heavy elements in the envelope (Alibert et al. 2018). In contrast, in the pebble-driven 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 so-called pebble-isolation 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 non-homogeneous 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: Three-layer and dilute-core 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 three-layer Jupiter model, where Jupiter’s interior consists of an outer H2-dominated region, an intermediate Hmetallic-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 (Zdilute), a parameter that controls the extent of the dilute core in terms of mass (mdilute), and a parameter controlling how steep the heavy-element gradient (δmdil) is. In the end, the ratio of the heavy-element mass mixing between the inner core outer boundary and the helium transition is given by Z=Z2+Zdilute−Z22[1−erf(m−mdiluteδmdil)],$Z = {Z_2} + {{{Z_{{\rm{dilute}}}} - {Z_2}} \over 2}\left[{1 - {\rm{erf}}\left({{{m - {m_{{\rm{dilute}}}}} \over {\delta {m_{{\rm{dil}}}}}}} \right)} \right],$(1)

with δmdil = 0.075. We tested possible variations of this parameter and found that the results of the calculation are not sensitive to the choice of δmdil. In addition, for all the models, we allowed for an increase in the interior entropy in the helium-demixing region by including a temperature jump (ΔTHe), 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 MH13-H (Militzer & Hubbard 2013; Miguel et al. 2016), CMS19-H (Chabrier et al. 2019), and MLS21-H (Mazevet et al. 2022) for the equation of state of hydrogen, and for helium we used SCVH95-He (Saumon et al. 1995) and CMS19-He (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 (T1bar). Voyager radio occultations (Lindal 1981) initially indicated T1bar = 165 ± 5 K. In situ measurements from the Galileo probe measured, in a so-called 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 3-5 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 T1bar 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 H2-dominated layer. We took these constraints into account and used a He abundance of Y1 = 0.238 (von Zahn et al. 1998) and metallicity Ζ1 = 1 Zsolar in our calculations (Li et al. 2020). For the Hmetalic -dominated region, we required that the overall abundance of He be consistent with the protosolar abundance (Yproto = 0.277, Serenelli et al. 2010) and allowed either different metallicity from the H2-dominated region or the same (the latter option for the dilute-core models).

thumbnail Fig. 1

Schematic view of the Jupiter structure models used in this paper. Panel a: three-layer model. Panel b: dilute-core model. The different regions are indicated with arrows.
thumbnail Fig. 2

Offsets between ToF and CMS for two subsam-ples (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 ΔJ2. 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 (J2n,n=1,…,5), which were ultimately compared with the Juno-derived 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 (J2nstatic)$\left({J_{2n}^{{\rm{static}}}} \right)$ by subtracting the computed dynamical contributions from the Juno measurements. Therefore, the effective gravity harmonics that our interior model calculations must match are J2nstatic=J2nJuno−ΔJ2ndifferential$J_{2n}^{{\rm{static}}} = J_{2n}^{{\rm{Juno}}} - \Delta J_{2n}^{{\rm{differential}}}$ (n = 1,…; 5), where J2nJuno$J_{2n}^{{\rm{Juno}}}$ are the values provided by Juno mission measurements (Iess et al. 2018; Durante et al. 2020) and ΔJ2ndifferential$\Delta J_{2n}^{{\rm{differential}}}$ 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 Gauss-Legendre 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 (δJ2n = J2n(ToF) - J2n(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 Mcore and Z1 to perfectly fit the observational measurements of the equatorial radius Req and J2 - 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 sub-samples (models with and without the optimisation procedure). Our results show a correlation between the offsets, where δJ4 and δJ6 depend strongly on δJ2. Higher order gravitational moments (δJ8, δJ10) are less dependent on δJ2. Thanks to these linear relationships, we can quantify the impact of an offset on J2 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 δJ2 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 (ΔJndifferential)$\left({\Delta J_{\rm{n}}^{{\rm{differential}}}} \right)$, we considered the widest possible range of flow profiles that still matched the gravity harmonics. The Juno gravity measurements revealed significant north-south asymmetries in Jupiter’s gravity field (Iess et al. 2018). Such hemispherical asymmetries (represented through the measurement of J3, J5, J7, and J9) can only be due to interior flows. The resemblance of the measured gravity signature of the flow to that obtained by extending the east-west (zonal) observed cloud-level 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 cloud-level wind structure) and its vertical profile (Galanti et al. 2019). This was achieved through randomly varying ΔJndifferential$\Delta J_{\rm{n}}^{{\rm{differential}}}$, 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 J3Juno$J_3^{{\rm{Juno}}}$ to J9Juno$J_9^{{\rm{Juno}}}$ and even differential contributions ΔJ2differential$\Delta J_{\rm{2}}^{{\rm{differential}}}$# to ΔJ10differential$\Delta J_{{\rm{10}}}^{{\rm{differential}}}$). 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 cloud-level 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: ΔJ2differential=1.039±0.354$\Delta J_{\rm{2}}^{{\rm{differential}}} = 1.039 \pm 0.354$, ΔJ4differential=−0.076±0.083$\Delta J_{\rm{4}}^{{\rm{differential}}} = - 0.076 \pm 0.083$, ΔJ6differential=0.016±0.076$\Delta J_{\rm{6}}^{{\rm{differential}}} = 0.016 \pm 0.076$, ΔJ8differential=0.053±0.062$\Delta J_{\rm{8}}^{{\rm{differential}}} = 0.053 \pm 0.062$, and ΔJ10differential=−0.080±0.042$\Delta J_{{\rm{10}}}^{{\rm{differential}}} = - 0.080 \pm 0.042$) (where ΔJ2ndifferential/10−6$\Delta J_{{\rm{2}}n}^{{\rm{differential}}}/{10^{ - 6}}$).

thumbnail 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. J2. 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.
Table 1

Offsets between ToF and CMS.

thumbnail Fig. 4

Contribution to the gravity harmonics from the deep atmospheric flow. Shown are the solutions for the cloud-level wind (a), the radial decay profile (b), and the wind-induced even gravity harmonics (c-f). 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 cloud-level wind, and the green line in (b) shows the solution in Kaspi et al. (2018). The red crosses in (e-f) show the expected mean value of the differential contribution.

3 Results

3.1 Envelope inhomogeneity

All of our models reproduce J2 and Jupiter’s radius, and Figs. 5 (panels a-c) and 6 (panels a-c) show the values of J4 and J6 for the MCMC solutions in the case of a three-layer model and with a dilute-core model, respectively (the other J2n 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 H2- and Hmetallic -dominated layers for the three-layer models and to consist of the H2 - and Hmetallic-dominated layers and the dilute core for the dilute-core models. In Fig. 5d we show results of the three-layer models, and we see that the difference in heavy elements between the H2- and Hmetallic-dominated layers (ΔΖ = Z2 - Z1) 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 H2-dominated and dilute-core regions (ΔΖ = Zdilute - Z1). 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 (MZtatal)$\left({{M_{{Z_{{\rm{tatal}}}}}}} \right)$ varies between 11 and 30 MEarth, with differences resulting from the choice of the equation of state. Calculations done with the MH13 equation of state have 18<MZtatal<30MEarth$18 < {M_{{Z_{{\rm{tatal}}}}}} < 30{M_{{\rm{Earth}}}}$, models with MLS21 have 14<MZtatal<24MEarth$14 < {M_{{Z_{{\rm{tatal}}}}}} < 24{M_{{\rm{Earth}}}}$, and models with CMS19 have 11<MZtatal<18MEarth$11 < {M_{{Z_{{\rm{tatal}}}}}} < 18{M_{{\rm{Earth}}}}$, independently of the model of Jupiter adopted (three layers or dilute core). For the three-layer models, the differences mostly arise from discrepancies in the mass of heavy elements in the Hmetallic-dominated region (MZ2)$\left({{M_Z}_{_2}} \right)$, which varies between 2 and 23 MEarth 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 dilute-core region (MZdilute)$\left({{M_{{Z_{{\rm{dilute}}}}}}} \right)$, which is found to vary between 1 and 25 MEarth depending on the equation of state adopted in the calculation. The mass of heavy elements in the H2-dominated region (MZ1)$\left({{M_{{Z_1}}}} \right)$ 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 MEarth for all models, with its exact value depending on the model adopted for Jupiter’s interior. The three-layer models have Mcore - 6 MEarth, 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 MEarth and another group with small masses of up to 2 MEarth.

Interestingly, models tend to have temperatures at 1 bar going from ~169 to ~188 K for the three-layer models to values of between ~165 and ~182K for the dilute-core 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 super-adiabaticity 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, dilute-core 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).

thumbnail Fig. 5

Details of models with three layers. Panels a, b, and c: J4 vs. J6 for three-layer models calculated using different equations of state for H, indicated in the figure. Cyan squares show the Juno measurements (J2nJuno)$\left({J_{2n}^{{\rm{Juno}}}} \right)$, and gravitational harmonics corrected by differential rotation - the ones to match with our interior models (J2nstatic)$\left({J_{2n}^{{\rm{static}}}} \right)$ - are shown in red. Panel d: histogram showing the difference in heavy elements (ΔZ = Z2 - Z1) between the Hmetallic- and H2-dominated layers.
thumbnail Fig. 6

Specifics of models with dilute cores. Panels a, b, and c: J4 vs. J6 obtained considering different equations of state for H. The Juno measurement (J2nJuno)$\left({J_{2n}^{{\rm{Juno}}}} \right)$ is indicated with the cyan square. The effective gravity harmonics to be matched by our interior models (J2nstatic)$\left({J_{2n}^{{\rm{static}}}} \right)$ are shown in red. Panel d: histogram that shows the difference in heavy elements (ΔZ = Zdilute - Z1) between the dilute-core region and the H2-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.

Dilute-core models were first modelled by Wahl et al. (2017), who, using the MH13 equation of state, found a value of Mcore + Mdilute between 10 and 24 MEarth and MZtotal${M_{{Z_{{\rm{total}}}}}}$ of 24-27 MEarth. These results are consistent with our results, even though most of them have sub-solar atmospheric metallicities.

Ni (2019) performed four-layer models using the MH13 and CMS19 equations of state. They found a range of Mcore+Mdilcore of 6.5-27 MEarth and of 24-28 MEarth with the MH13 equation of state and Mcore + Mcore of 3-12 MEarth and MZtotal${M_{{Z_{{\rm{total}}}}}}$ of 8-12 MEarth 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 sub-solar atmospheric metallicities as opposed to our results, where Z1 = 1 Zsolar.

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 H2- and Hmetallic-dominated regions was required. We also find these types of solutions when using the CMS19 equation of state and three-layer models (Fig. 5), but, since they represent only 2.4% of our sample, we find it to be unlikely. The MZtotal${M_{{Z_{{\rm{total}}}}}}$ found by Debras & Chabrier (2019) are 25-30 and 40-45 MEarth 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 dilute-core models. They found Mcore of up to 3.8 MEarth and MZtotal${M_{{Z_{{\rm{total}}}}}}$ up to 13 MEarth, 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 H2- and Hmetallic-dominated regions close to 6 Mbar, while our models use separation pressures of around 3 Mbar.

thumbnail 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 three-layer 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 H2-dominated region in the y-axis, panels b and e show the mass in the Hmetallic-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. 170-180 K rather than 166K). We robustly demonstrate that the heavy-element abundance is not homogeneous in Jupiter’s envelope. Our results imply that Jupiter continued to accrete heavy elements in large amounts while its hydrogen-helium envelope was growing, contrary to predictions based on the pebble-isolation mass in its simplest incarnation (Lambrechts et al. 2014), favouring instead planetesimal-based 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 non-adiabatic interior models for the giant planets, and it provides a base example for exoplanets: a non-homogeneous 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 MEarth, 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 helium-rain 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 θ={[Mcore,Z2rock,Z2ice,T1bar,PHe,ΔTHe]three-layer[Mcore,Zdiluterock,Zdiluteice,T1bar,PHe,ΔTHe]dilute-core,$\theta = \left\{ {\matrix{ {\left[{{M_{core}},Z_2^{{\rm{rock}}},Z_2^{{\rm{ice}}},{T_{1{\rm{bar}}}},{P_{{\rm{He}}}},\Delta {T_{{\rm{He}}}}} \right]}{{\rm{three - layer}}} \cr {\left[{{M_{core}},Z_{{\rm{dilute}}}^{{\rm{rock}}},Z_{{\rm{dilute}}}^{{\rm{ice}}},{T_{1{\rm{bar}}}},{P_{{\rm{He}}}},\Delta {T_{{\rm{He}}}}} \right]}{{\rm{dilute - core,}}} \cr } } \right.$(A.1)

and the data we allowed to vary are X=[J2static,J4static,J6static,J8static,J10static,Reqobs]$X = \left[{J_2^{static},J_4^{static},J_6^{static},J_8^{static},J_{10}^{static},R_{eq}^{obs}} \right]$(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 p(θ|X)∝p(θ)p(X|θ).$p\left({\theta \left| X \right.} \right) \propto p\left(\theta \right)p\left({X\left| \theta \right.} \right).$(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, X=J(θ)+ϵ.$X{ = _J}\left(\theta \right) + .$(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 p(X|θ)∝∏iexp[−(Xi−Ji(θ))22σi2].$p\left({X\left| \theta \right.} \right) \propto \prod\limits_i {\exp \left[{ - {{{{\left({{X_i}{ - _{Ji}}\left(\theta \right)} \right)}^2}} \over {2\sigma _i^2}}} \right]} .$(A.5)

The σi2$\sigma _i^2$ 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, p(θ)∝∏ipi(θi).$p\left(\theta \right) \propto \prod\limits_i {{p_i}\left({{\theta _i}} \right)} .$(A.6)

The chosen priors are given in Tables A.1 and A.2. The boundaries of all priors were chosen using a test-and-trial 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).

Table A.1

Parameters explored in our MCMC calculations for three-layer 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. Yproto=0.277, Y1=0.238, and Z1 =0.0153.

Table A.2

Parameters explored in our MCMC calculations for dilute-core 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. Yproto = 0.277, YH2=0.238${Y_{{H_2}}} = 0.238$, Z1 = 0.0153, and Z2 = 0.0153.

thumbnail Fig. A.1

MCMC corner plot showing the posterior distribution of variables obtained with the three-layer model and the MH13 equation of state. Red points with error bars show the observed parameters as a reference.
thumbnail 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.
thumbnail Fig. A.3

Distribution results of the MCMC runs with three layers and the MLS20 equation of state.
thumbnail Fig. A.4

Corner plot resulting from the runs with the dilute-core model and the MH13 equation of state. Red dots show the observed parameters.
thumbnail Fig. A.5

Distribution of different parameters resulting from the runs with the dilute-core model and the CMS 19 equation of state.
thumbnail Fig. A.6

Resulting distributions obtained with the dilute-core model and the MLS20 equation of state.

References

  1. Alibert, Y., Venturini, J., Helled, R., et al. 2018, Nat. Astron., 2, 873 [NASA ADS] [CrossRef] [Google Scholar]
  2. Bazot, M., Bourguignon, S., & Christensen-Dalsgaard, J. 2012, MNRAS, 427, 1847 [NASA ADS] [CrossRef] [Google Scholar]
  3. Brygoo, S., Loubeyre, P., Millot, M., et al. 2021, Nature, 593, 517 [NASA ADS] [CrossRef] [Google Scholar]
  4. Chabrier, G., Mazevet, S., & Soubiran, F. 2019, ApJ, 872, 1 [NASA ADS] [CrossRef] [Google Scholar]
  5. Christen, J. A., & Fox, C. 2010, Bayesian Anal., 5, 263 [CrossRef] [Google Scholar]
  6. Debras, F., & Chabrier, G. 2019, ApJ, 872, 100 [Google Scholar]
  7. Duer, K., Galanti, E., & Kaspi, Y. 2020, J. Geophys. Res. (Planets), 125, 8 [Google Scholar]
  8. Durante, D., Parisi, M., Serra, D., et al. 2020, Geophys. Res. Lett., 47, e86572 [NASA ADS] [CrossRef] [Google Scholar]
  9. Fletcher, L. N., Orton, G. S., Greathouse, T. K., et al. 2020, J. Geophys. Res., 125, 8 [Google Scholar]
  10. Galanti, E., & Kaspi, Y. 2016, ApJ, 820, 91 [NASA ADS] [CrossRef] [Google Scholar]
  11. Galanti, E., & Kaspi, Y. 2017, Icarus, 286, 46 [NASA ADS] [CrossRef] [Google Scholar]
  12. Galanti, E., Kaspi, Y., Miguel, Y., et al. 2019, Geophys. Res. Lett., 46, 2 [Google Scholar]
  13. Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
  14. Guilera, O. M., Sándor, Z., Ronco, M. P., et al. 2020, A&A, 642, A140 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Guillot, T., & Morel, P. 1995, A&AS, 109, 109 [NASA ADS] [Google Scholar]
  16. Guillot, T., Miguel, Y., Militzer, B., et al. 2018, Nature, 555, 227 [Google Scholar]
  17. Guillot, T., Li, C., Bolton, S., et al. 2020, J. Geophys. Res., 125, 8 [Google Scholar]
  18. Gupta, P., Atreya, S., Steffes, P. G. et al. 2022, PSJ, submitted [Google Scholar]
  19. Hubbard, W. B. 2013, ApJ, 768, 43 [NASA ADS] [CrossRef] [Google Scholar]
  20. Hubbard, W. B., & Militzer, B. 2016, ApJ, 820, 80 [NASA ADS] [CrossRef] [Google Scholar]
  21. Iess, L., Folkner, W. M., Durante, D., et al. 2018, Nature, 555, 220 [NASA ADS] [CrossRef] [Google Scholar]
  22. Kaspi, Y. 2013, Geophys. Res. Lett., 40, 676 [NASA ADS] [CrossRef] [Google Scholar]
  23. Kaspi, Y., Davighi, J. E., Galanti, E., & Hubbard, W. B. 2016, Icarus, 276, 170 [NASA ADS] [CrossRef] [Google Scholar]
  24. Kaspi, Y., Guillot, G., Galanti, E., et al. 2017, Geophys. Res. Lett., 44, 12 [Google Scholar]
  25. Kaspi, Y., Galanti, E., Hubbard, W. B., et al. 2018, Nature, 555, 223 [Google Scholar]
  26. Kaspi, Y., Galanti, E., Showman, A. P., et al. 2020, Space Sci. Rev., 216, 84 [Google Scholar]
  27. Kulowski, L., Cao, H., & Bloxham, J. 2020, J. Geophys. Res.: Planets, 125 [Google Scholar]
  28. Lambrechts, M., Johansen, A., & Morbidelli, A. 2014, A&A, 572, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  29. Leconte, J., Selsis, F., Hersant, F., & Guillot, T. 2017, A&A, 598, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Li, C., Ingersoll, A., Bolton, S., et al. 2020, Nat. Astron., 4, 609 [Google Scholar]
  31. Lindal, G. F., et al. 1981, J. Geophys. Res., 86, 8721 [NASA ADS] [CrossRef] [Google Scholar]
  32. Lyon, S., & Johnson, J. 1992, LANL Report LA-UR-92-3407, Los Alamos [Google Scholar]
  33. Mazevet, S., Licari, A., & Soubiran, F. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/201935764 [Google Scholar]
  34. Miguel, Y., Guillot, T., Fayon, L. 2016, A&A, 596, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  35. Militzer, B., & Hubbard, W. B. 2013, ApJ, 774, 148 [NASA ADS] [CrossRef] [Google Scholar]
  36. Morales, M. A., Hamel, S., Caspersen, K., & Schwegler, E. 2013, Phys. Rev. B, 87, 174105 [NASA ADS] [CrossRef] [Google Scholar]
  37. Nettelmann, N. 2017, A&A, 606, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Nettelmann, N., Movshovitz, N., Ni, D., et al. 2021, Planet. Sci. J., 2, 241 [NASA ADS] [CrossRef] [Google Scholar]
  39. Ni, D. 2019, A&A, 632, A76 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Ormel, C. W., Vazan, A., & Brouwers, M. G. 2021, A&A, 647, A175 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
  42. Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
  43. Seiff, A., Kirk, D. B., Knight, T. C. D., et al. 1998, J. Geophys. Res., 103, 22857 [NASA ADS] [CrossRef] [Google Scholar]
  44. Serenelli, A. M., & Basu, S. 2010, ApJ, 719, 865 [NASA ADS] [CrossRef] [Google Scholar]
  45. Schöttler, M., Redmer, R. 2018, Phys. Rev. Lett., 120 [Google Scholar]
  46. Tollefson, J., Wong, M. H., Pater, J., et al. 2017, Icarus, 296, 163 [NASA ADS] [CrossRef] [Google Scholar]
  47. Valletta, C., & Helled, R. 2019, ApJ, 871, 127 [NASA ADS] [CrossRef] [Google Scholar]
  48. Vazan, A., Helled, R., & Guillot, 2018, A&A, 610, L14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Venturini, J., & Helled, R. 2020, A&A, 634, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. von Zahn, U., Hunten, D. M., & Lehmacher, G. 1998, J. Geophys. Res., 103, 22815 [NASA ADS] [CrossRef] [Google Scholar]
  51. Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649 [CrossRef] [Google Scholar]
  52. Zharkov, V. N., & Trubitsyn, V. P. 1978, Phys. Planet. Inter. Astron. Astrophys. Ser., Tucson: Pachart [Google Scholar]

All Tables

Table 1

Offsets between ToF and CMS.

In the text
Table A.1

Parameters explored in our MCMC calculations for three-layer 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. Yproto=0.277, Y1=0.238, and Z1 =0.0153.

In the text
Table A.2

Parameters explored in our MCMC calculations for dilute-core 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. Yproto = 0.277, YH2=0.238${Y_{{H_2}}} = 0.238$, Z1 = 0.0153, and Z2 = 0.0153.

In the text

All Figures

thumbnail Fig. 1

Schematic view of the Jupiter structure models used in this paper. Panel a: three-layer model. Panel b: dilute-core model. The different regions are indicated with arrows.
In the text
thumbnail Fig. 2

Offsets between ToF and CMS for two subsam-ples (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 ΔJ2. 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
thumbnail 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. J2. 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
thumbnail Fig. 4

Contribution to the gravity harmonics from the deep atmospheric flow. Shown are the solutions for the cloud-level wind (a), the radial decay profile (b), and the wind-induced even gravity harmonics (c-f). 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 cloud-level wind, and the green line in (b) shows the solution in Kaspi et al. (2018). The red crosses in (e-f) show the expected mean value of the differential contribution.
In the text
thumbnail Fig. 5

Details of models with three layers. Panels a, b, and c: J4 vs. J6 for three-layer models calculated using different equations of state for H, indicated in the figure. Cyan squares show the Juno measurements (J2nJuno)$\left({J_{2n}^{{\rm{Juno}}}} \right)$, and gravitational harmonics corrected by differential rotation - the ones to match with our interior models (J2nstatic)$\left({J_{2n}^{{\rm{static}}}} \right)$ - are shown in red. Panel d: histogram showing the difference in heavy elements (ΔZ = Z2 - Z1) between the Hmetallic- and H2-dominated layers.
In the text
thumbnail Fig. 6

Specifics of models with dilute cores. Panels a, b, and c: J4 vs. J6 obtained considering different equations of state for H. The Juno measurement (J2nJuno)$\left({J_{2n}^{{\rm{Juno}}}} \right)$ is indicated with the cyan square. The effective gravity harmonics to be matched by our interior models (J2nstatic)$\left({J_{2n}^{{\rm{static}}}} \right)$ are shown in red. Panel d: histogram that shows the difference in heavy elements (ΔZ = Zdilute - Z1) between the dilute-core region and the H2-dominated layer.
In the text
thumbnail 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 three-layer 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 H2-dominated region in the y-axis, panels b and e show the mass in the Hmetallic-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
thumbnail Fig. A.1

MCMC corner plot showing the posterior distribution of variables obtained with the three-layer model and the MH13 equation of state. Red points with error bars show the observed parameters as a reference.
In the text
thumbnail 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
thumbnail Fig. A.3

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

Corner plot resulting from the runs with the dilute-core model and the MH13 equation of state. Red dots show the observed parameters.
In the text
thumbnail Fig. A.5

Distribution of different parameters resulting from the runs with the dilute-core model and the CMS 19 equation of state.
In the text
thumbnail Fig. A.6

Resulting distributions obtained with the dilute-core model and the MLS20 equation of state.
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.