New Models of Jupiter in the Context of Juno and Galileo

Our Jupiter models are calculated with the concentric Maclaurin spheroid method (Hubbard 2012, 2013). As demonstrated in Debras & Chabrier (2018), in order to yield valid models, the method must fulfill several mathematical and numerical constraints in terms of numbers and spacing of the spheroids and of the treatment of the outermost spheroids. Accordingly, the spheroids implemented in our calculations are spaced exponentially, their equatorial radius is ${\lambda }_{i}=1\,-({e}^{i\beta }-1)/({e}^{N\beta }-1)$ with N the number of spheroids, i ranging from 0 to $N-1$, β = 6/N, and the upper atmosphere is neglected.³ In this paper, we examine which type of model is compatible with the Juno observations, provided the difference can be explained by the maximum allowed amount of differential rotation, i.e., differential rotation penetrating down to 10,000 km (see Guillot et al. 2018 and Kaspi et al. 2017). Said differently, we wish that the uncertainties on the J values obtained for an acceptable model are smaller than the uncertainty that is due to this maximum possible level of differential rotation. At this level, we checked that 512 spheroids yield a sufficient precision and that using 1000 spheroids or changing the β parameter does not significantly affect the conclusions. Deriving more precise models that precisely meet all Juno's and Galileo's constraints with lower levels of differential rotation, however, requires at least 1000 spheroids to ensure that the discretization error is negligible compared with the other sources of error on the evaluation of the gravitational moments. The various parameters used for Jupiter throughout this work are reported in Table 1.

Table 1.  Values of the Planetary Parameters of Jupiter

Parameter	Value
G (a) (global parameter)	$6.672598\times {10}^{-11}\pm 2\times {10}^{-17}$ m3 kg−1 s−2
$G\times {M}_{J}$(b)	(126686533 ± 2) ×109 m3 s−2
MJ	1.89861 × 1027 kg
Req (c)	71492 ± 4 km
Rpolar (c)	66854 ± 10 km
ω (d)	1.7585324 × 10−4 ± 6 × 10−10 s−1
$\bar{\rho }$	1326.5 kg m−3
$m=3{\omega }^{2}/4\pi G\bar{\rho }$	0.083408
$q={\omega }^{2}{R}_{\mathrm{eq}}^{3}/{{GM}}_{J}$	0.0891954
${J}_{2}\times {10}^{6}$ (e)	14696.572 ± 0.014
$-{J}_{4}\times {10}^{6}$ (e)	586.609 ± 0.004
J6 × 106 (e)	34.198 ± 0.009
−J8 × 106 (e)	2.426 ± 0.025
J10 × 106 (e)	0.172 ± 0.069

Note. R_eq and R_polar are observed at 1 bar. The pulsation value is chosen following Archinal et al. (2011). (a) Cohen & Taylor (1987), (b) Folkner et al. (2017), (c) Archinal et al. (2011), (d) Riddle & Warwick (1976), and (e) Iess et al. (2018).

Download table as:  ASCIITypeset image

Throughout this work, we use for the H/He mixture a combination of the new EOS recently derived by Chabrier et al. (2019), based on semi-analytical models in the low (molecular/atomic) and high (fully ionized) temperature-density domains and quantum molecular dynamic (QMD) calculations in the intermediate pressure dissociation/ionization regime, and the Militzer & Hubbard (2013) EOS that takes into account non-ideal correlation effects. As in Miguel et al. (2016), we first calculated a pure H table by calculating at each P–T point

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{{\rho }_{\mathrm{MH}13}} & = & \displaystyle \frac{{X}_{\mathrm{MH}13}}{{\rho }_{{\rm{H}}}}+\displaystyle \frac{{Y}_{\mathrm{MH}13}}{{\rho }_{\mathrm{He},\mathrm{New}}}\Rightarrow \displaystyle \frac{1}{{\rho }_{{\rm{H}}}}\\ & = & \displaystyle \frac{1}{{X}_{\mathrm{MH}13}}\left(\displaystyle \frac{1}{{\rho }_{\mathrm{MH}13}}-\displaystyle \frac{{Y}_{\mathrm{MH}13}}{{\rho }_{\mathrm{He},\mathrm{New}}}\right),\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}\begin{array}{rcl}{S}_{\mathrm{MH}13} & = & {X}_{\mathrm{MH}13}{S}_{{\rm{H}}}+{Y}_{\mathrm{MH}13}{S}_{\mathrm{He},\mathrm{New}}+{S}_{\mathrm{mix}}\Rightarrow {S}_{{\rm{H}}}\\ & = & \displaystyle \frac{1}{{X}_{\mathrm{MH}13}}\left({S}_{\mathrm{MH}13}-{Y}_{\mathrm{MH}13}{S}_{\mathrm{He},\mathrm{New}}-{S}_{\mathrm{mix}}\right),\end{array}\end{eqnarray} \tag{ 2 }$$

where ρ_H is the mass density for pure hydrogen, ρ_MH13 is the density derived from MH13 by spline procedures, ρ_He,New is the helium density in the Chabrier et al. (2019) EOS, S_H is the sought pure hydrogen specific entropy, S_MH13 is the splined specific entropy from MH13, ${S}_{\mathrm{He},\mathrm{New}}$ is the helium-specific entropy in the new EOS, all at the same (P,T), and X_MH13 = 0.7534, Y_MH13 = 0.2466 are the mass fractions of hydrogen and helium in the MH13 simulations. Finally, S_mix is the mixing-specific entropy defined as

$$\begin{eqnarray}&&\displaystyle \frac{{S}_{\mathrm{mix}}}{{k}_{b}}=\displaystyle \frac{1}{{M}_{{\rm{H}},\mathrm{He}}}\left[{N}_{{\rm{H}}}\mathrm{ln}\left(1+\displaystyle \frac{{N}_{\mathrm{He}}}{{N}_{{\rm{H}}}}\right)+{N}_{\mathrm{He}}\mathrm{ln}\left(1+\displaystyle \frac{{N}_{{\rm{H}}}}{{N}_{\mathrm{He}}}\right)\right],\end{eqnarray} \tag{ 3 }$$

with N_H and N_He the numbers of H and He particles, respectively, of number fractions ${x}_{i}={N}_{i}/({N}_{{\rm{H}}}+{H}_{\mathrm{He}})$ ($i\equiv {\rm{H}}$ or He), ${M}_{{\rm{H}},\mathrm{He}}=\bar{A}{m}_{{\rm{H}}}$ the total mass of H+He, $\bar{A}={\sum }_{i}{x}_{i}{A}_{i}$ the mean atomic number, and ${m}_{{\rm{H}}}=1.660\times {10}^{-27}$ kg the atomic mass unit. This "mixed" (Chabrier et al./MH13) pure hydrogen EOS is then combined with the new pure helium EOS (Chabrier et al. 2019; F. Soubiran et al. 2019, in preparation) to obtain a complete EOS for the H/He mixture at any given helium mass fraction Y.

Figure 1 displays the relative error on the density between the EOS used by us or Miguel et al. (2016) and the EOS of MH13, $(\rho -{\rho }_{\mathrm{MH}13})/{\rho }_{\mathrm{MH}13}$, for Y = 0.2466, which is the helium fraction used in MH13. For Miguel et al. (2016), we have combined their published pure H table with a He table from SCvH with a cubic order spline. The comparisons are made for 32 (T, P) points from MH13 corresponding to an entropy characteristic of the Jupiter interior, 7–8 k_B/proton. These points are used as inputs in the mixed EOS used by us or Miguel et al. (2016) to calculate the corresponding density and entropy, which is to be compared with the EOS of MH13. As seen in the figure, above $500\,\mathrm{kg}\,{{\rm{m}}}^{-3}$, which is the pressure ionization domain in Jupiter, the difference between our results and those of MH13 is always <0.5%, which is smaller than the numerical error in MH13, whereas for the EOS of Miguel et al. (2016), the differences are significant. This is a major issue in the present context, where a very accurate density profile is required to derive reliable gravitational moments.

Zoom In Zoom Out Reset image size

For the heavy elements in the H/He-rich envelope, which is essentially composed of volatiles (H₂O, CH₄, and NH₃), we use a recent EOS for water that is based on QMD simulations at high density (Licari 2016; Mazevet et al. 2018). Here also, this EOS has been shown to adequately reproduce available Hugoniot experiments (Knudson et al. 2012). Given the small number fraction of CH₄ and NH₃ compared with H, He and even H₂O, we do not expect the assumption to generically treat all the volatiles with the water EOS to be consequential on the results. When considering a diluted core, we combine the above water EOS with the Sesame "drysand" EOS (Lyon & Johnson 1992) to take into account additional heavy elements such as silicates and iron. Finally, when including a central compact core, we use a 100% "drysand" EOS. We verified that using an EOS for pure iron instead of the drysand EOS for this region, for instance, does not make a noticeable difference. Note that ab initio simulations have shown that under the characteristic temperature and pressure conditions typical of the deep interior of Jupiter, water is in a liquid state and is fully soluble in metallic hydrogen (Wilson & Militzer 2012b). Similarly, solid SiO₂ and MgO, which are representative examples of planetary silica and rocky material, are also found to be soluble in H⁺ under similar conditions (Wilson & Militzer 2012a; Wahl et al. 2013; González-Cataldo et al. 2014). These thermodynamic considerations support a core erosion for typical central conditions on Jupiter and thus a mixed ${\rm{H}}/\mathrm{He}/Z$ eos in such a region.

As shown by Soubiran & Militzer (2016), the inclusion of heavy elements in a H/He/Z mixture under Jupiter-like internal temperature and density conditions can be performed with the so-called additive-volume-law (AVL) provided we use an effective volume (density) for the heavy species. We verified that our EOS is consistent with the work of Soubiran & Militzer (2016), and hence that our water EOS, which is representative of volatiles in Jupiter, can be used throughout the entire T–P domain from Jupiter's atmosphere to the center.

For H/He/Z mixtures, our EOS are then combined at each given (P,T) point throughout the AVL:

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{1}{\rho } & = & \displaystyle \frac{1-{Z}_{\mathrm{water}}-{Z}_{\mathrm{drysand}}}{{\rho }_{{\rm{H}},\mathrm{He}}}+\displaystyle \frac{{Z}_{\mathrm{water}}}{{\rho }_{\mathrm{water}}}\\ & & +\,\displaystyle \frac{{Z}_{\mathrm{drysand}}}{{\rho }_{\mathrm{drysand}}}\quad \mathrm{at}\,(P,T)=\mathrm{constant},\end{array}\end{eqnarray} \tag{ 4 }$$

where ${\rho }_{{\rm{H}},\mathrm{He}}$, ${\rho }_{\mathrm{water}}$ and ${\rho }_{\mathrm{drysand}}$ are the densities of the H/He mixture, water and drysand, respectively, and ${Z}_{\mathrm{water}}\,={M}_{\mathrm{water}}/M$, ${Z}_{\mathrm{drysand}}={M}_{\mathrm{drysand}}/M$ the mass fractions of water and drysand, respectively, with M the mass of the planet. Note that the accuracy of the AVL for the hydrogen/water mixture under the relevant conditions for Jupiter interior has been verified with QMD simulations (Soubiran & Militzer 2015).

Given the small number fraction of heavy elements compared with H and He, the P and T used to calculate the densities in the H/He/Z mixture are those obtained with the H/He mixture only. Similarly, the entropy of heavy elements can be neglected (see Soubiran & Militzer 2016), and even when their mass fraction becomes Z ≳ 0.2, they affect the total mixing entropy by at most 2%, which represents a few per thousand of the total entropy. Moreover, this occurs only in the deepest part of the planet, with little impact on the gravitational moments. Hence, the total entropy is evaluated as the entropy of a pure H–He mixture with effective hydrogen and helium mass fractions ${X}^{\mathrm{eff}}=X/(1-Z)$ and ${Y}^{\mathrm{eff}}=Y/(1-Z)$, respectively, with ${X}^{\mathrm{eff}}+{Y}^{\mathrm{eff}}=1$, and $X={M}_{{\rm{H}}}/M$, $Y\,={M}_{\mathrm{He}}/M$, and $Z={M}_{Z}/M={Z}_{\mathrm{water}}+{Z}_{\mathrm{drysand}}$.

For the outer element abundances, the observations of Galileo give

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{{Y}_{\mathrm{Gal}}}{({X}_{\mathrm{Gal}}+{Y}_{\mathrm{Gal}})} & = & 0.238\pm 0.005,\\ {Z}_{\mathrm{Gal}} & = & 0.0167\pm 0.006,\end{array}\end{eqnarray*}$$

where ${X}_{\mathrm{Gal}}$ and ${Y}_{\mathrm{Gal}}$ are the observed mass abundances of hydrogen and helium, respectively. ${Z}_{\mathrm{Gal}}$ is the abundance of heavy elements in the high envelope measured by Galileo, but the real abundance of heavy elements should be larger. This implies that ${Y}_{\mathrm{Gal}}$ and ${X}_{\mathrm{Gal}}$ are only defined relative to each other and that ${X}_{\mathrm{Gal}}+{Y}_{\mathrm{Gal}}+{Z}_{\mathrm{Gal}}$ can be larger than 1. In all the following models of this paper, except if stated otherwise, we impose the external atmosphere to have helium and heavy element mass fractions of

$$\begin{eqnarray*}\begin{array}{rcl}{Y}_{\mathrm{ext}} & = & 0.23\\ {Z}_{\mathrm{ext}} & = & 0.02,\end{array}\end{eqnarray*}$$

which corresponds to ${Y}_{\mathrm{ext}}/({X}_{\mathrm{ext}}+{Y}_{\mathrm{ext}})=0.2347$. As just mentioned, this Z value is most likely a lower limit for the heavy element content in the external envelope of Jupiter (see Section 5.4 for a detailed exploration of this issue). Neglecting Galileo's constraints, i.e., reducing the observed amount of heavy elements, drastically reduces the constraints on the models and allows the derivation of a large range of models that are compatible with Juno data. Relaxing these constraints thus drastically simplifies the calculations of models consistent with only Juno observations. Such simplifications, as done in all recent studies (Wahl et al. 2017; Guillot et al. 2018, except when using an ad hoc modification of their EOS), however, can hardly be justified (see Section 6). In all our calculations, the mean helium fraction of the planet is fixed to the protosolar value: $\bar{Y}/(\bar{X}+\bar{Y})={Y}_{\odot }=0.275$ (see, e.g., Anders & Grevesse 1989).

We first calculate an isentropic model, composed of one homogeneous convective isentropic gaseous envelope, with $\bar{Y}=0.275$, and a spherical compact core of constant density. The total heavy material content is determined to obtain the correct mass of the planet. This model reproduces the J₂ observed value within 10⁻⁷, the maximum intrinsic precision of the CMS method (Debras & Chabrier 2018). The J₄ and J₆ values are compared to those of Juno in Figure 2 under the labels "Isentrope." The differences between the observed and calculated values are about 3% and 6%, respectively, which is well above any numerical source of error.

Zoom In Zoom Out Reset image size

The only possibility to reconcile the observed and theoretical values would be strong differential rotation that affects the calculation of the gravitational moments, which is not implemented in the CMS calculations. However, the results of Kaspi et al. (2017) show that with various flow profiles extending more than 10,000 km within the planet, the change in J₄ is at most 0.7% (see their Figure 4). Furthermore, the study of Cao & Stevenson (2017) excludes such a deep differential rotation.

• 1.  
  This yields the first robust conclusion: Jupiter's interior is not isentropic.

The next step in increasing complexity is to change the composition, then the entropy, somewhere in the planet. There are two physically plausible domains: a diluted core extending throughout a substantial fraction of the interior, and/or a region of either layered convection, hydrogen metallization, or H/He phase separation somewhere within the envelope.

3.2.1. Two Possible Locations

3.2.2. Results

Following up on the above analysis, we have explored two types of models with entropy and compositional changes either in the central region (the "diluted core") or in the gaseous envelope, as schematically illustrated in Figure 3. In case of a diluted core that is unstable to double-diffusive behavior, Moll et al. (2017) showed that a central seed could survive erosion longer than the lifetime of Jupiter. The fact that under Jupiter central T and P conditions, iron and silicates are in solid form (see, e.g., Musella et al. 2018) tends to favor the presence of such a central seed. We thus consider the presence of a compact core at the very center of the planet. For the entropy variation within the gaseous envelope, we have considered either an abrupt (ΔS) or a gradual (δS) change in entropy and composition at the metallization boundary, between at most 0.1 and 3.0 Mbar. The obtained J₄ and J₆ values are given in Figure 2. For the models with a change of composition in the gaseous envelope, the values for a gradual or a sharp entropy change are plotted under the labels "Gradual" and "Sharp," respectively.

Zoom In Zoom Out Reset image size

In case of a gradual (continuous) change, which implies a continuous molecular H₂ to metallic H⁺ transition and no H/He phase separation, the smooth change in entropy is simply due to a composition change (see Equations (2)–(4)). In case of a first-order molecular-metallic transition, the abrupt change in entropy ΔS is used as a free parameter, discretized over a certain number of spheroids to obtain the proper J values.

In all cases, none of these two types of models, regardless of the type of change in composition and entropy, sharp or gradual, was found to be able to yield J₄ and J₆ values sufficiently close to the observed values (labeled "Juno") to be explained by differential rotation or deep winds.

Another "simple" possible interior structure model is the model suggested by Leconte & Chabrier (2012): the entire planet would be made of alternating convective and diffusive layers. These authors, however, pointed out that this entire double-diffusive interior model could be replaced by a model with a localized double-diffusive buffer in the envelope, surrounded by large-scale convective envelopes (see their Section 4.3), similar to the type of model explored above, which is found to be excluded. We return to this point in Sections 5 and 6.

This yields the following conclusions:

• 1.  
  Conclusion 1: The first conclusion of this section is that the models of Jupiter displayed in Figures 3(a) and (b) cannot fulfill the observational constraints of both Juno and Galileo. A mix of these two types of models is required: Jupiter is at least composed of an envelope split in two parts (an outer molecular/atomic envelope with Galileo element composition and an inner ionized envelope) separated by a region of compositional change, and a diluted core extending throughout a significant fraction of the planet. A compact core can also be present at the center of the planet.
• 2.  
  Conclusion 2: In order to decrease the values of $| {J}_{4}|$ and $| {J}_{6}|$ in the models of Figure 3, we realized that either ΔS had to be substantial (an issue explored in Section 5), or the heavy element content must decrease with depth in the outer part of the planet. This local decrease in Z is balanced by an increase in Y so that the density (and the molecular weight, see Section 6), of course, increases with depth. As mentioned in the introduction, one of the most stringent constraints on the models are the Y and Z values observed by Galileo, which are surprisingly high for the observed values of the high-order gravitational moments. A local inward decrease in metal content in some region of the planet's gaseous envelope appears to be the favored solution to resolve this discrepancy. This is examined in detail in the next section.

It is essential at this stage to stress the crucial role played by the H/He EOS. Using the SCvH EOS, Chabrier et al. (1992) and subsequent similar models, which fulfilled the observational constraints of Galileo and Voyager, could relatively easily fulfill those of Juno as well. This is entirely due to the SCvH EOS (or, similarly, R-EOS; Nettelmann et al. 2012), see Figure 11 of Militzer & Hubbard (2013) or Figure 27 of Chabrier et al. (2019): for a given entropy, such an EOS has a lower density (pressure) in the ∼Mbar region than our new EOS, enabling a larger amount of heavy element repartition in the interior of Jupiter, which appreciably relaxes the constraints on possible models. The constraints become much more stringent with stiffer, more accurate EOSs.

Two physical processes can lead to a locally decreasing abundance of heavy elements with depth in Jupiter's outer envelope, i.e., a locally positive gradient ∇Z > 0: a "static" process, based on thermodynamic stability criteria, and a "dynamic" process, which involves non-equilibrium processes. Both these processes can be enhanced by an external accretion event, which we also discuss in the following (Section 4.1.3).

4.1.1. Thermodynamic Stability

4.1.2. Upward Atomic Motions

4.1.3. Accretion

To be considered as plausible, our models with a region of locally inward-decreasing abundance of heavy elements must be consistent with what is known of Jupiter's long-term evolution. If, as expected, Jupiter formed through core accretion (Pollack et al. 1996), the primordial abundance of heavy elements in the planet should be increasing with depth (see Figure 4). As explored by, e.g., Leconte & Chabrier (2013) and Vazan et al. (2018), the differential core-envelope cooling of the planet leads to a redistribution of heavy elements with time, yielding an Increasing heavy element content in the gas-rich envelope. To explain our ∇Z > 0, a physical process must have inhibited convection within the envelope and prevented a homogeneous redistribution of elements. Three possibilities have been discussed in the previous subsections: first-order metallization, immiscibility, and accretion. We examine whether they are compatible with the evolution of the planet.

Zoom In Zoom Out Reset image size

If a first-order metallization (PPT) stopped the convective motions, our discussion on the bubbles in Section 4.1.2 shows that it is possible to deplete the inner envelope and enrich the outer envelope. This would be in adequation with any evolutionary scenario.

In the case of hydrogen/helium immiscibility, the outer envelope will be depleted in helium because of helium sedimentation (see Section 4.1.1, Stevenson & Salpeter 1977a), yielding an enrichment in $Z={M}_{Z}/M$ in this region. Note that in order to obtain a locally steep enough gradient of heavy elements ∇Z > 0 between the outer Z-enriched and the inner Z-depleted envelopes, the He-rich falling dropplets must be as Z-poor as possible. This implies limited miscibility between some Z-components and He, in addition to H and He, and typical Z/He phase diagrams yielding very low concentrations of heavy elements in He-rich dropplets, as discussed in Section 4.1.1. Under such conditions, it is possible to preserve (and even increase) a positive Z-gradient with time.

Finally, in the absence of PPT or immiscibility, if inhibition of convection, leading to layered convection (Leconte & Chabrier 2012) or even partly radiative interiors during the planet growth, is the only reason for the difference in composition between the outer and inner envelopes, it must have persisted since the accretion event(s). Although, as examined in Section 6.2, the present models fulfill the constraints required for the onset and persistence of layered convection (Leconte & Chabrier 2012, 2013), whether such structures can persist during Jupiter's or in fact any gaseous planet's cooling history (Chabrier & Baraffe 2007) needs to be explored with extreme care (e.g., Rosenblum et al. 2011; Mirouh et al. 2012; Wood et al. 2013; Kurokawa & Inutsuka 2015) and requires that the key physical processes at play are handled with great accuracy. A fantastic challenge for numerical simulations.

To conclude this section, we should mention that recent evolutionary calculations (e.g., Vazan et al. 2018) converge to a structure profile for the current Jupiter with a monotonically outward-decreasing compositional gradient, ∇Z < 0, which is at odds with our suggestion of a local ∇Z > 0. These models, however, use the SCvH EOS. As discussed in Section 3.2.2, this latter yields a too hot thermal structure along an adiabat (see Militzer & Hubbard 2013 or Figure 27 of Chabrier et al. 2019), favoring convection and allowing larger metal fractions. These models of the internal structure of Jupiter are excluded by our Juno+Galileo analysis and thus cannot be used as reliable evolutionary constraints. Further evolutionary calculations with the proper physics, including a proper treatment of double-diffusive convection and of H/He phase separation, are needed to verify the consistency of our models with the thermal history of Juipter. As just mentioned, however, properly handling such complex physical processes is a task of major difficulty (see Section 6.5).

As briefly examined in Section 3.2, a brutal inward increase of S can have several physical foundations. Assuming that Jupiter's outermost thermal profile is isentropic (because of adiabatic convection in this region), the observed condition, T = 165 K, P = 1 bar, according to our EOS corresponds to T ≃ 5000 K at 1 Mbar. As mentioned above, recent first-principle simulations (Mazzola et al. 2018) predict a critical temperature for the metallization of hydrogen in the range T_c ≃ 2000–5000 K at P ≃ 1 Mbar. Simulations (Soubiran et al. 2013; Mazzola et al. 2018) and experiments (Loubeyre et al. 1985) both suggest, however, that even for a low helium concentration as is found in Jupiter (x_He < 0.1), the critical pressure increases while the critical temperature decreases with increasing helium concentration, which probably excludes a PPT between molecular and metallic hydrogen in Jupiter. However, given the current uncertainties in these determinations, we must still explore this possibility.

In case of a first-order transition, ΔS is given by the Clausius-Clapeyron relation along the critical line P(T)

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}S}{{\rm{\Delta }}(1/\rho )}=\displaystyle \frac{{dP}}{{dT}}.\end{eqnarray} \tag{ 5 }$$

Analytical calculations (Saumon & Chabrier 1992) suggest

$$\begin{eqnarray}&&{\rm{\Delta }}{S}_{\mathrm{metallization}}\sim 0.5\,{k}_{{\rm{B}}}/\mathrm{proton}.\end{eqnarray} \tag{ 6 }$$

Because, as mentioned above, the temperature in this region of Jupiter's interior should be close to T_c, we expect ΔS to be lower than this value.

If hydrogen pressure ionization does not occur through a first-order transition inside Jupiter, a sharp entropy change can be due to H/He phase separation (also a first-order transition). As shown by Stevenson & Salpeter (1977a), drowning nucleated helium dropplets lead to a release of gravitational energy, and even though their analysis suggests that most of this energy is radiated away, part of it contributes to heating up the inner part of the planet, raising the entropy (see the detailed discussion in Section 6.5).

The shape of the H/He phase diagram is a major uncertainty in this context. The rather limited helium depletion with respect to the solar value in Jupiter's external envelope, x_He ≃ 0.1, suggests that the variation ΔY in the immiscible region should be modest (about ∼10%). In that case, according to Figure 2 of Morales et al. (2013b), the mixing entropy should depart only slightly from the ideal mixing entropy, by ∼0.03 k_B/at at 5000 K for x_He = 0.1. Because the maximum value of the ideal mixing entropy for a concentration x_He = 0.5 is ${S}_{\mathrm{mix}}^{\mathrm{id}}/N=-[x\mathrm{ln}x+(1-x)\mathrm{ln}(1-x)]=0.7\,{k}_{B}$/at (about 0.3 k_B/at for x_He = 0.1), we see that the entropy jump due to H/He immiscibility should be ≲0.5 k_B/proton. The entropy change due to helium dropplet sedimentation is more difficult to evaluate and requires numerical explorations. Guidance is provided by the calculations of Fortney & Hubbard (2003) for the case of Saturn. In the case of a maximum temperature gradient in the inhomogeneous region and no formation of a helium layer on top of the core (both the most likely current situation), these authors find that a change in composition $Y=0.21\to 0.36$ corresponds to a global increase in entropy ΔS ∼ 0.3 k_B/proton. For Jupiter, we expect helium sedimentation (i) to have occurred, if ever, more recently than for Saturn, (ii) to encompass a much smaller fraction of the planet (see Section 6.1) and thus to induce a much smaller entropy variation. Adding up these two contributions, it seems difficult to justify an entropy jump arising from H/He phase separation much larger than

$$\begin{eqnarray}&&{\rm{\Delta }}{S}_{{\rm{H}}/\mathrm{He}}\lesssim 1.0\,{k}_{{\rm{B}}/\mathrm{proton}}.\end{eqnarray} \tag{ 7 }$$

Clearly, more experimental and numerical exploration of hydrogen pressure metallization, the H/He phase diagram, and the He sedimentation process are strongly needed to help constrain these processes.

Finally, if the mean molecular weight gradient that is due to the change in composition is large enough to hamper adiabatic convection, a regime of layered convection can develop and lead to a superadiabatic temperature structure that is similar, at least in some part of the planet, to the structure derived in Leconte & Chabrier (2012). The detailed treatment of layered convection in our calculations is presented in Section 6.2. By varying the location and extent of layered convection between 0.1 Mbar and 2 Mbar, i.e., in the vicinity of hydrogen pressure ionization, we numerically obtain a maximum entropy increase from layered convection,

$$\begin{eqnarray}&&{\rm{\Delta }}{S}_{\mathrm{layered}}\lesssim 0.6\,{k}_{B}/\mathrm{proton},\end{eqnarray} \tag{ 8 }$$

with a decreasing metal abundance, i.e., ∇Z > 0, in this region. An increasing metal abundance in this region yields higher values of ΔS, but in this case, α ≲ 10⁻⁷, where α denotes the ratio of the size of the convective layer to the pressure scale height, α = l/H_P (see Section 6). This implies the presence of a diffusive buffer or a regime of turbulent diffusion within the ionization boundary layer. Although detailed calculations are lacking, it seems difficult to reconcile such a structure with Jupiter's thermal history.

In summary, if either a first-order transition and/or layered convection is present within some part (most likely around the Mbar) of the gas-rich planetary envelope, an inhomogeneous zone where the total increase in entropy can reach $\sim 0.5\mbox{--}1\,{k}_{B}/\mathrm{proton}$ is expected.

Figure 5 shows the typical structure of our final Jupiter models. They all share the following features:

• 1.  
  An outer homogeneous convective envelope characterized by the Galileo helium and heavy element abundances and the adiabat condition at T = 165 K, P = 1 bar.
• 2.  
  An inhomogeneous region between ∼0.1 and ∼2 Mbar associated with (i) a change in composition, most likely characterized by an inward-decreasing metal abundance (∇Z > 0), and (ii) a non-negligible entropy (≳5%) and temperature increase. These gradients stem from layered convection and/or H–He immiscibility, even though a PPT cannot be entirely excluded for now.
• 3.  
  An inner homogeneous convective envelope lying on a warmer isentrope than the outer region, with an average larger helium fraction and most likely a lower metal fraction than in the outermost region. Indeed, even though a larger Z fraction in this region than in the outer region is not entirely excluded, it requires an uncomfortably large entropy increase, ΔS ≳ 1.4 k_B/proton (case (c)), according to the estimates derived above. As shown below, there is a degeneracy between the entropy jump in the inhomogeneous region and the helium and metal fractions in the inner envelope. The larger Y and Z in the inner envelope, the larger ΔS needs to be.
• 4.  
  A diluted (eroded) core extending throughout a significant fraction of the planet. A small entropy jump in the inhomogeneous region, ΔS ≲ 0.5 k_B/proton (case (a)), yields an inward-increasing helium abundance in the core, while a higher value (case (b)) implies an inward-decreasing helium abundance in the core.
• 5.  
  Probably, but not necessarily, a central compact, solid core.

Zoom In Zoom Out Reset image size

A quantitative analysis of these models is given in the next subsection.

The results of our optimized models with an entropy discontinuity ${\rm{\Delta }}S\in [0\mbox{--}2]{k}_{{\rm{B}}}/\mathrm{proton}$ in the ionization boundary region, projected in different J_k–${J}_{k+2}$ plans, are displayed in Figure 6 for two values of the external heavy element abundance, namely ${Z}_{\mathrm{ext}}={Z}_{\odot }$ and ${Z}_{\mathrm{ext}}=2\times {Z}_{\odot }$. The first obvious conclusion from this figure is that our range of models consistent with the gravitational moments observed by Juno differs from the range derived by Guillot et al. (2018) with 200,000 models (see their Figure 1 of the Extended Data, reported as a hashed area in Figure 6). We have verified that this is not a discretization issue: with interior structures calculated with 1000 spheroids, this conclusion is unaltered. Even though the difference between the two analyses should partly stem from the different EOS used by these authors, it arises essentially from our different representation of the planet interior. Indeed, we recall that these authors did not take the constraint from Galileo on the heavy element abundances into account. Therefore, if Galileo's observations are correct, Figure 6 shows that the qualitative conclusions that these authors draw about differential rotation could be altered.

Zoom In Zoom Out Reset image size

1. Impact on the low-order moments (≤J₄). We found out that the entropy change, ΔS, is strongly affected by the composition in the inner convective envelope, i.e., the region between the inhomogeneous envelope and the diluted core. Therefore the size and composition of the diluted core, the composition of the inner envelope, and the entropy change ΔS in the region of compositional variation are intrinsically linked. To better understand this result, we must recall that the main problem of the models is to decrease J₄ at constant J₂. Figure 7(a) shows the value of the contribution function ${J}_{2}-{J}_{4}$ in the planet as a function of pressure.

Zoom In Zoom Out Reset image size

This figure shows that in order to decrease J₄ with respect to J₂, one needs to enrich the planet deeper than ∼2 Mbar, and the region around ∼10 Mbar is where it is most efficient. Therefore, an enriched inner envelope decreases J₄ at constant J₂ (see the pressure range in Figure 7(a)); but enriching the inner envelope implies a steeper compositional gradient in the boundary region between the outer and inner envelopes, which has the opposite effect on J₄ compared to J₂. Furthermore, this boundary region between ∼0.1 and 2 Mbar has a much stronger contribution on J₂ and J₄ than the deeper region. This stems from the fact that this region has a large mean radius, hence the mass of a spherical shell is much higher than in the 5 Mbar region, and the impact on J₂ and J₄ is enhanced. In consequence, a small change in the ∼0.1–1 Mbar region must be compensated for by a strong change in the diluted core.

2. Impact on the high-order moments (>J₄). The high-order gravitational moments strongly depend on the value of J₄, as shown by Guillot et al. (2018). For a given J₄, the other parameters affecting these moments are the external abundance of metals (as expected), and the mass of the central compact core. Changing the helium content within the inner convective envelope has almost no impact, as there is a trade-off between the inner abundances of helium and heavy elements and the entropy increase, without affecting the high-order gravitational moments. Similarly, the position and extent of the boundary region of compositional change is a second-order correction to the J₆ to J₁₀ values. As a whole, we found out that the J₆ to J₁₀ values are not much affected by the composition in the inner part of Jupiter, deeper than where the compositional change occurs.

As mentioned above, we found out that some models with an inward increase of heavy elements in the inhomogeneous region (∇Z < 0) of the envelope can fulfill all observational constraints (case (c)) if the entropy change around ∼1 Mbar reaches values ${\rm{\Delta }}S\gt 1{k}_{{\rm{B}}}/\mathrm{proton}$. This requires a strong entropy discontinuity induced either by a PPT for T ≃ 4500 K (thus a critical point T_c ≫ 4500 K) or by H/He differentiation and sedimentation, as layered convection alone cannot yield such an entropy jump. Therefore, although not entirely excluded, models with ∇Z < 0 throughout the entire envelope are rather uncomfortable, as discussed in Section 5.1. In contrast, models with an inward-decreasing abundance of heavy elements (∇Z > 0) in this region require a more modest entropy change.

The fact that, surprisingly, the mass of the compact core affects the high-order moments can be explained as follows. Because we consider the central compact core as spherical, it has no direct influence on the gravitational moments. However, in that case, a smaller fraction of the planetary mass is available to satisfy the J values. Since the outer envelope composition is constrained by Galileo, one can only enrich the inner envelope or the diluted core to compensate for this. Figure 7(a) shows, however, that if the increase in density deeper than 2 Mbar is too strong, the increase in J₂ is larger than the increase in J₄ (and even larger than the increase in higher-order moments, not shown). This leads to

conclusion 1: for given J₂, a central compact core tends to decrease the moments of order ≥J₄ compared to a model with only a diluted core.

Increasing the mass of the compact core thus implies adding more heavy elements in the inner regions of the planets, diluted core, or inner envelope to reproduce the J values. The larger the amount of heavy elements in the deep layers, however, the larger the required entropy jump ΔS between the outer and inner envelopes. This in turn has consequences for the high-order moments: a model with a high ΔS in the metallization region implies a higher temperature and thus a lower density for a given composition at a given pressure in the Mbar region than a model with a smaller ΔS. The lower density means that this region has a smaller contribution to the gravitational moments than in the case of a small ΔS. Although for J₂ and J₄, this can be balanced by a higher contribution of the internal layers (deeper than a few Mbars), this is not the case for the higher-order moments.

To illustrate this result, we show in Figure 7(b) the differences ΔJ_n in the contribution functions J₂ to J₈ between a model with a small entropy change, $\delta S\lt 0.1\,{k}_{{\rm{B}}}/\mathrm{proton}$, in the ionization boundary region that is only due to a composition change, and a model with a total entropy discontinuity ${\rm{\Delta }}S=0.9\,{k}_{{\rm{B}}}/\mathrm{proton}$. We see that the region of the outer (molecular) envelope always has a stronger contribution to the Js when the entropy change is small, as expected, regardless of the order of the moment. On the other hand, the inner region of the diluted core has a far greater impact on J₂ than on the other moments. A strong ΔS thus requires more heavy elements in the diluted core to preserve J₂, while the high-order moments are almost insensitive to the enhanced composition in the diluted core. This leads to

conclusion 2: an entropy jump in the envelope tends to decrease the value of the high-order moments at a given J₂. A large enough entropy change is necessary to preserve the correct balance between the moments.

As discussed in Section 5.1, however, the possible entropy increase in the metallization boundary region is limited by physics principles. For central compact cores larger than $\gtrsim 5\,{M}_{\oplus }$, one needs ${\rm{\Delta }}S\gt 1\,{k}_{{\rm{B}}}/\mathrm{proton}$, which, as discussed above, seems to be hardly possible at these temperatures. Figure 8 displays the values of the high-order gravitational moments as a function of the entropy jump ΔS. Low (absolute) values of J₆, J₈, or J₁₀ allways require a significant ΔS, except when we decrease the atmospheric ${Z}_{\mathrm{ext}}$ and in that case violate the Galileo constraint, as done in all recent studies. Models with no entropy jump in the gaseous envelope thus seem to be excluded, as mentioned previously.

Zoom In Zoom Out Reset image size

As seen in Figure 8, none of our models can match the 3σ error bars on J₆ for ${\rm{\Delta }}S\lt 1\,{k}_{{\rm{B}}}/\mathrm{proton}$ when the contribution from the winds derived in Kaspi et al. (2018) is considered. This is particularly true when the external abundance of heavy elements Z_ext is supersolar (see Section 5.4). When the dynamical correction from Kaspi et al. (2017) is considered, however, flows extending down to 3000 km are sufficient to explain the discrepancy with the observed gravitational moments. Therefore, either the ΔJ₆ correction due to the wind contribution in Kaspi et al. (2018) is underestimated because of an erroneous estimation of the winds or the presence of north-south symmetric zonal flows that will affect the even gravitational moments, or the entropy increase must reach at least $\sim 1.5\,{k}_{{\rm{B}}}/\mathrm{proton}$. In any case, a continuously increasing heavy element mass fraction with depth, i.e., ∇Z < 0, in the Mbar region is hard to justify (in Figure 8, these models all have ${\rm{\Delta }}S\gtrsim 1\,{k}_{{\rm{B}}/\mathrm{proton}}$).

As shown in Figure 5, the valid models predict a size for the metallization boundary region, l_b ≈ 15% of Jupiter's radius. Clearly, this is orders of magnitude larger than any possible interface due to a PPT. It can, however, be consistent with the size of the inhomogeneous H/He region, as this latter continues to expand during the planetary cooling. Finally, as shown in Section 6.2, this region, which is characterized by a compositional gradient, is prone to layered convection, which by itself is characterized by a superadiabaticity and thus by its own entropy variation, to be added to the convection issued from a phase transition, and thus contributing to the total ΔS.

Figure 9 portrays the thermodynamic and composition profiles of our models that are consistent with all Galileo and Juno constraints, taking into account for this latter the correction due to differential rotation from Kaspi et al. (2017, 2018), respectively. Profiles for an isentropic interior structure, which of course is inconsistent with the observed gravitational moments, are shown for comparison. The blue curves represent our favored models, with J₆ compatible with Kaspi et al. (2017) but not with Kaspi et al. (2018), while the red curves are the profiles obtained from a model with a lower J₆, which is at the limit of what can be reached according to Kaspi et al. (2018). Globally, the pressure and density profiles differ by a few percent at most from those of the isentropic model, which are barely visible in the figure. However, it is worth stressing that the density of the optimized model is lower in the Mbar region than the density of the isentropic model, whereas the opposite is true in the central regions (diluted and compact core). This is a direct consequence of the constraints arising from the gravitational moments and the Galileo observations, as it allows decreasing the J₄ to J₁₀ values for the correct J₂. In contrast, the temperature departs from the isentropic profile for R ≲ 0.9 × R_J, i.e., within most of the interior, by a difference ΔT ≃ +1000–2000 K. Interestingly enough, this temperature increase agrees very well with the value obtained by Fortney & Hubbard (2003) in the H/He inhomogeneous layer for a helium enrichment in the interior from Y = Y_⊙ to Y = 0.35, and a temperature gradient leading to overstable convection. As a consequence, the specific entropy increases from the outer to the inner envelope. This increase is steeper for the model with a lowered J₆ (consistent with Kaspi et al. 2018). For this latter, the inner isentropic envelope occupies a very limited fraction of the planet, 0.89 × R_J ≲ R ≲ 0.92 × R_J, and the diluted core extends up to 85% of the planet. Within the diluted core, the specific entropy decreases drastically because of the strong increase in heavy elements (Figure 9(b)). We do not show the specific entropy and compositional profiles in the diluted core because of the degeneracy between helium and heavy element distributions in this region, which yield similar results for the gravitational moments, in addition to the fact that the entropy profile is of no real interest in this region. The important parameter is the steepness of the gradient of composition between the inner envelope and the diluted core. The steeper the gradient, the smaller the diluted core needs to be to obtain the correct J₂. The mean heavy element mass fraction is displayed in Figure 9(b). As discussed earlier, Z decreases between the outer and inner envelopes (∇Z > 0) because models with a continuously increasing Z (∇Z < 0) in the envelope, although not strictly excluded, require a very large entropy jump (ΔS > 1k_B/proton), which is difficult to reconcile with the examined physical processes (Section 5.1). Future work in this direction will certainly help to clarify this issue.

Zoom In Zoom Out Reset image size

Figure 10(a) shows the corresponding mass profile of our typical optimized interior structure of Jupiter that fulfills the constraints of Juno and Galileo with the wind correction of Kaspi et al. (2017). An isentropic profile is shown for comparison. The black circle indicates the inner limit of the outer convective zone, the two crosses bracket the inner convective zone, and the diamond corresponds to the limit of the compact core, if present. The zoom on the right-hand side shows the inner convective zone, encompassing about 30% of the mass of the planet. The heavy element distribution is displayed in Figure 10(b). For models with no central compact core, the total amount of heavy elements in the planet is M_Z = 25–30 M_⊕. Adding up a compact core yields up to M_Z = 40–45 M_⊕.

Zoom In Zoom Out Reset image size

As discussed in Section 3, we have taken a very conservative lower limit for the true average metal content in the atmosphere, ${Z}_{\mathrm{ext}}$, in our calculations. We have taken a solar value, while the abundances measured by Galileo of individual heavy elements, excluding oxygen and neon, instead yield ${Z}_{\mathrm{Gal}}^{\mathrm{true}}\approx (2\mbox{--}3)\times {Z}_{\odot }\approx 0.04\mbox{--}0.06$. Interestingly, such a highly oversolar value seems to be supported even for oxygen by the latest observations of the great red spot (Bjoraker et al. 2018). Increasing ${Z}_{\mathrm{ext}}$ for a given structure increases $| {J}_{6}|$ to $| {J}_{10}|$ and thus implies either a very strong differential rotation or a very large ΔS to preserve the moments. If ${Z}_{\mathrm{ext}}\simeq 0.05$, models with a constant inward increase of Z in the metallization boundary region lead to a ΔS much larger than the 1 ${k}_{B}/\mathrm{proton}$ maximum value derived before that was consistent with physical estimates. This reinforces our previous conclusion:

• 1.  
  The internal structures of Jupiter with an inward increase of heavy elements within the Mbar boundary region imply uncomfortable physical constraints: the entropy jump or amount of differential rotation required to be compatible with the high-order gravitational moments needs to be very large. In contrast, models with a locally decreasing abundance of heavy elements within Jupiter's metallization boundary region fulfill all constraints with acceptable levels of entropy variation and differential rotation.

The J values for five models with ${Z}_{\mathrm{ext}}\gt {Z}_{\odot }$ are shown in Figure 8 (orange circles). We see that for a given ΔS, these models have higher $| {J}_{6}|$ to $| {J}_{10}|$ values than models with ${Z}_{\mathrm{ext}}={Z}_{\odot }$. Although some of these models are compatible with the correction that is due to the differential rotation estimated in Kaspi et al. (2017), they are hardly compatible with the observations when the correction to the even gravitational moments estimated in Kaspi et al. (2018) is considered. We recall, however, that all the models of Figure 8 have J₂ and J₄ values that are consistent with Kaspi et al. (2018). Because of the strong correlation between J₄ and J₆, further decreasing $| {J}_{4}|$, which would then be consistent with Kaspi et al. (2017) but not with Kaspi et al. (2018), would allow us to decrease the $| {J}_{6}|$ and higher-order moment values and expand the range of plausible models. The derivation of precise constraints on the depth penetration of differential rotation and its effect on the Js as a function of ${Z}_{\mathrm{ext}}$ will be examined in a subsequent paper.

First, following the nomenclature of Stevenson & Salpeter (1977a), we have assumed that Jupiter had a "hot start," meaning that the initial inner temperature of the planet was higher than the critical temperature of both hydrogen metallization through a PPT, ${T}_{c}({\rm{H}}\mbox{--}{{\rm{H}}}_{2})$, and H/He demixion, ${T}_{c}({\rm{H}}\mbox{--}\mathrm{He})$ (for ${x}_{\mathrm{He}}=0.08$). According to all existing numerical simulations aimed at exploring these issues, this is quite a safe assumption. Further work on the metallization of hydrogen and the H/He phase diagram will enable us to discriminate between the sectors I, II, and III of Figure 1 of these authors, namely:

• 1.  
  Sector I: if ${T}_{c}({\rm{H}}\mbox{--}{{\rm{H}}}_{2})\lt T(P)\lt {T}_{c}({\rm{H}}\mbox{--}\mathrm{He})$, where T(P) is the local temperature at pressure P, hydrogen metallization is occurring smoothly but probably triggers H/He or Z_i/He immiscibility. The only possibility of depleting (enriching) the inner (outer) envelope in metals is to invoke two- or three-body immiscibility diagrams between partially pressure-ionized heavy elements ${Z}^{n+}$ and neutral He, similar to what is occurring for H⁺–He phase separation, as explored in Section 4.1.1, and/or external accretion events, depending on the exact value of ${Z}_{\mathrm{ext}}$.
• 2.  
  Sector II: ${T}_{c}({\rm{H}}\mbox{--}{{\rm{H}}}_{2})\sim {T}_{c}({\rm{H}}\mbox{--}\mathrm{He})$. As for the Sector I case, an inhomogeneous region forms that is depleted in He and some Z-components, but in that case, there is a possibly of H⁺-rich bubble nucleation and thus uplifting He-poor, Z-rich bubbles and dropping He-rich dropplets (Section 4.1.2).
• 3.  
  Sector III: if ${T}_{c}({\rm{H}}\mbox{--}{{\rm{H}}}_{2})\gt {T}_{c}({\rm{H}}\mbox{--}\mathrm{He})$, H/He demixion has not yet started, the redistribution is due to the H⁺-rich, He-poor bubbles described before. This is probably the most unlikely situation.

These situations are imposed by the necessity of globally increasing the metal content of the upper envelope (and conversely, deplete the lower envelope) to fulfill the Galileo constraints, ${Z}_{\mathrm{ext}}\geqslant {Z}_{\mathrm{Gal}}$, but also to enrich the inner helium content, Y, to balance the Z decrease. According to current work on metallization and immiscibility of hydrogen and helium, even though substantial uncertainty remains, and if, as found in numerical simulations, hydrogen (or any heavy component) ionization triggers immiscibility with He atoms (or He-like atoms), the Sector I case is the most likely one. This shows the need for numerical explorations of this type of diagram and, more generally, of the stability of H/He/Z mixtures for the internal conditions of Jupiter.

Noticeable differences still exist between modern ab initio calculations aimed at characterizing the H/He phase diagram. Figure 11 shows the immiscibility region predicted by some of these calculations with the T–P profiles obtained with our favored models that fulfill all constraints of Galileo and Juno, taking into account either the Kaspi et al. (2018) (low J₆) or Kaspi et al. (2017) (high J₆) correction due to differential rotation. Figure 11(a) corresponds to interior structures with strongly superadiabatic layered convection occurring at P ≥ 0.1 Mbar. Figure 11(b) displays two models (labeled "Morales" and "Lorenzen," respectively) for which the change of entropy is only due to the H/He phase separation, i.e., where it occurs at the corresponding critical pressures, without any layered convection above this layer. A Jupiter isentropic profile is shown for comparison. As seen in the figure, while according to the Lorenzen et al. (2009) calculations, H/He phase separation could take place in some fraction of our favored Jupiter interior models, this is not the case with the results of Morales et al. (2013b; or Schöttler & Redmer 2018, not shown), which predict no H/He immiscibility in the current Jupiter. For the models with no layered convection above the phase separation (dash–dotted lines in Figure 11(b)), the temperature gradient is probably too high for overstable modes to persist, and convection will prevail (see, e.g., Figure 3 of Stevenson & Salpeter 1977a). Although the lack of excess (non-ideal) mixing entropy in the calculations of Lorenzen et al. (2009) casts doubt on the reliability of their phase diagram, it is worth noting that a ∼200–800 K underestimation of the critical temperature in the 1–2 Mbar domain by Morales et al. (2013b; no temperature error bar is shown in these calculations) would be consistent with immisciblity for our model with Y₂ = 0.25. Therefore, the current internal structure of Jupiter could entail a region of layered convection starting around ∼0.1 Mbar, associated with some change in composition, and a (probably small) region of H/He (most likely H/He/Z) immiscibility at deeper levels. Although more numerical exploration of this major issue is certainly needed, key diagnostics on H/He phase separation under the relevant conditions might be provided by existing experiments (Soubiran et al. 2013).

Zoom In Zoom Out Reset image size

As seen in Figure 11, it seems difficult to reconcile a H/He phase separation, according to the most recent calculations, with a model that reproduces the Kaspi et al. (2018) J₆ value. Furthermore, the required entropy increase for this model leads to such a steep temperature gradient that unstable convection will prevail. It is thus very difficult to justify the very large entropy change ${\rm{\Delta }}S$ in the gaseous envelope of Jupiter that is required in this model on physical grounds. This suggests either a revision of the Kaspi et al. (2018) analysis, or the presence of north-south symmetric winds, which are inconsequential for the odd gravitational moments, but would increase the correction to the even gravitational moments, rejoining the corrections obtained in Kaspi et al. (2017).

As found out in the previous sections, fulfilling the constraints of Galileo and Juno while preserving a global mean helium protosolar value $\bar{Y}={Y}_{\odot }$ and a physically acceptable entropy increase ΔS in the hydrogen metallization region requires an inward decrease of the heavy element abundance in this region, i.e., a locally positive gradient, ∇Z > 0. We verified that because of the ∼1/40 ratio of heavy elements to helium number, this region still exhibits a positive molecular weight gradient, ${{\rm{\nabla }}}_{\mu }=(d\mathrm{log}\mu /d\mathrm{log}P)\gt 0$. In that case, large-scale adiabatic convection can be inhibited and lead to a regime of small-scale superadiabatic double-diffusive convection (also called semi-convection) to transport heat. As mentioned previously, although a first-order transition is not required to trigger such a process, it strongly favors it, as suggested for instance at the Earth's mantle boundary (e.g., Christensen & Yuen 1985).

The condition for the onset of double-diffusive convection reads (e.g., Stern 1960)

$$\begin{eqnarray}&&0\lt {{\rm{\nabla }}}_{T}-{{\rm{\nabla }}}_{\mathrm{ad}}\lt \displaystyle \frac{{\alpha }_{\mu }}{{\alpha }_{T}}{{\rm{\nabla }}}_{\mu },\end{eqnarray} \tag{ 9 }$$

where ${({\alpha }_{\mu }=(\partial \mathrm{ln}\rho /\partial \mathrm{ln}\mu )}_{P,T}$ and ${\alpha }_{T}={(\partial \mathrm{ln}\rho /\partial \mathrm{ln}T)}_{P,\mu }$. In geophysics, it is well known that double-diffusive convection generally takes the form of oscillatory convection or layered convection, i.e., a stack of small-scale convective layers of size l separated by diffusive interfaces (e.g., Rosenblum et al. 2011). In astrophysical objects, because they are essentially of the lower Prandtl number, this layering is more blurred, and according to simulations, double-diffusive convection instead takes the form of ill-defined turbulent interfaces, even though finite-amplitude layering remains a possibility (Moll et al. 2016). In the absence of simulations in this context, we use the analytical formalism derived by Leconte & Chabrier (2012) to verify the presence of layered convection in our models. As shown by these authors, this is controlled by the parameter α, which is the ratio of the size of the convective layer to the pressure scale height, α = l/H_P. From their Equation (21), we can relate this parameter to the superadiabatic gradient, $(\langle {{\rm{\nabla }}}_{{\rm{T}}}\rangle -{{\rm{\nabla }}}_{\mathrm{ad}})$, by

$$\begin{eqnarray}&&\langle {{\rm{\nabla }}}_{{\rm{T}}}\rangle -{{\rm{\nabla }}}_{\mathrm{ad}}={\epsilon }_{d}\times \left[{\left({{\rm{\Phi }}}_{0}{\alpha }^{4}{\epsilon }_{d}\right)}^{\tfrac{-1}{4(1+a)}}+{\left({{\rm{\Phi }}}_{0}{\alpha }^{4}{\epsilon }_{d}\right)}^{\tfrac{-a}{(1+a)}}\right],\end{eqnarray} \tag{ 10 }$$

with layered convection occurring when

$$\begin{eqnarray}&&{10}^{-9}\mbox{--}{10}^{-6}\lesssim \alpha \lesssim {10}^{-4}\mbox{--}{10}^{-2},\end{eqnarray} \tag{ 11 }$$

(see Table 1 of Leconte & Chabrier 2012).

As MacLaurin spheroids by definition have a constant density, layered convection cannot be prescribed very accurately with the CMS method. As for the case of a first-order phase transition/separation, we have implemented a sharp entropy and composition change at constant T and P between consecutive layers. We can then verify for the appropriate models whether conditions (11) and (9) are fulfilled. Figure 12(a) displays the values of $({\alpha }_{\mu }/{\alpha }_{T})$, calculated with our EOS, for an isentropic profile and for two profiles with an entropy increase in the Mbar region of 0.5 and $1\,{k}_{{\rm{B}}}/\mathrm{proton}$, respectively. We see that this quantity increases with depth by an order of magnitude, between ${\alpha }_{\mu }/{\alpha }_{T}=1$ in the external layers and the values prevailing at depth in Jupiter, due essentially to the onset of of H₂ dissociation (see Chabrier et al. 2019). This favors the onset of layered convection deeper than ∼0.1 Mbar (see Equation (9)). Figure 12(b) displays the values of α and of the parameter ${R}_{\rho }^{-1}=({\alpha }_{\mu }/{\alpha }_{T}){{\rm{\nabla }}}_{\mu }/({{\rm{\nabla }}}_{T}-{{\rm{\nabla }}}_{\mathrm{ad}})$ (overstable convection occurs for ${R}_{\rho }^{-1}\gt 1$, see Rosenblum et al. 2011; Leconte & Chabrier 2012; Mirouh et al. 2012), for a composition change from (${Y}_{\mathrm{ext}}=0.23,{Z}_{\mathrm{ext}}=0.03$) to (${Y}_{2}\,=0.3,{Z}_{2}=0.01$), with an entropy increase ${\rm{\Delta }}S=0.45\,{k}_{{\rm{B}}/\mathrm{proton}}$ between 0.1 and 1 Mbar. We see that α and ${R}_{\rho }^{-1}$ fulfill the conditions for the presence of layered convection in this domain. Models with higher ΔS (of $\sim 0.6\,{k}_{{\rm{B}}}/\mathrm{proton}$ at most, see Section 5.1) require larger ΔY. Globally, we verified that all our favored models do fulfill the conditions for the occurrence of layered convection derived in Leconte & Chabrier (2012).

Zoom In Zoom Out Reset image size

One word of caution should be noted: when H₂ dissociates into atomic H⁺, the mean molecular weight μ decreases very strongly. According to Nellis et al. (1995), for instance, the fraction of dissociation is about 10% at 1.4 Mbar, however. The molecular weight thus remains barely affected up to this pressure and the decrease in μ due to H₂ dissociation should occur over a rather limited region between ∼1.4 and 2 Mbar. Whether layered convection is still present in this domain is less clear (although overshoot probably occurs), but we consider it to be localized enough to not significantly modify the temperature profile calculated above.

We examine here whether the high heavy element abundances reported by Galileo in the outermost part of the envelope might be due to recent impact or recently accreted material that would not yet have had time to be redistributed within the planet and thus would not affect its gravitational potential. In that case, the Juno contraints could be examined without taking into account those from Galileo. Using standard equations of the mixing length theory (Kippenhahn & Weigert 1990; Hansen & Kawaler 1994), the typical convective velocity in Jupiter scales as

$$\begin{eqnarray}&&{v}_{\mathrm{conv}}\approx 10{\left(\displaystyle \frac{{\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}}}{{10}^{-8}}\right)}^{1/2}\left(\displaystyle \frac{{H}_{P}}{{10}^{6}\,{\rm{m}}}\right)\,{\rm{m}}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 12 }$$

with (∇ − ∇_ad) ranging from ∼10⁻¹⁰ at the center to ∼10⁻⁶ near the surface, yielding v_conv ∼ a few to about ∼100 cm s⁻¹ from the center to the surface.

This means that within at most a few years, the external extra material will be mixed throughout most of the planet. It is quite clear that Jupiter has not accreted a few Earth masses of heavy elements in the past 20 years. The only source of uncertainty is Shoemaker–Levy 9. Although its mass is ridiculously low compared to the mass of the envelope and even to the mass of heavy elements in the envelope of Jupiter, the energy deposited when it crashed onto the planet (at the end of July 1994) triggered an uplifting of deep material (e.g., Bézard et al. 2002 or Moreno et al. 2003). As Galileo entered Jupiter's atmosphere 1.5 years later (7 December 1995), the material had probably be mixed again throughout the upper envelope (we recall that v_conv ≈ 1 ${\rm{m}}\,{{\rm{s}}}^{-1}$ in this region). Therefore, it seems to be mandatory to take the Galileo observations of Jupiter's heavy element external abundance into account when we try to recover the gravitational moments reported by Juno.

Atmospheric dynamical effects (e.g., a vortex), on the other hand, could have produced a localized maximum concentration of heavy material that is not representative of the (lower) average value. Additional information is provided by the Juno microwave instrument (Li et al. 2017), which suggests a rather unexpected vertical amonia profile. However, when we compare Figure 4 of these authors with the location at which the Galileo entry probe dived (around 5°–10° in planetocentric latitude), we see that compared to the region deeper than 100 bar, Galileo can only have measured a lower limit of the amonia content. Indeed, there seems to be a maximum amonia concentration in the equatorial regions, followed by a minimum in the mid-latitudes. The Galileo entry probe dived at the limit between these two regions, and most importantly, the concentration anywhere above the 25 bar level is lower than the concentration deeper than 100 bars, where convective mixing definitely takes place. Therefore it seems unlikely that the measurements of Galileo, consistent with an increase in ammonia with depth of up to 25 bar, are an upper bound of the heavy element composition.

One can also wonder whether the composition varies between the P = 1 bar and P = 100 bar levels, affecting the gravitational moments, and then should be parameterized instead of being assumed to be constant. Between two models with (Y, Z) = (0.23, 0.02) and (Y, Z) = (0.23, 0.04) in the external envelope, we obtain changes in J₂ and J₄ due to the first 100 bar variation of at most 2 × 10⁻⁶ and 2 × 10⁻⁷, respectively, which is about an order of magnitude smaller than the change that is due to differential rotation (Figure 3 of Kaspi et al. 2017). Variations in composition in the external envelope are thus a second-order correction compared to differential rotation.

Although the magnetic field at the surface of Jupiter has been shown to vary with latitude and longitude within an order of magnitude (Connerney et al. 2018), the leading feature is a dipolar field with moment M = 4.170 × 10⁻⁴ T. According to numerical simulations, self-consistent dynamo action is generally found to start when the convective magnetic Reynolds number, i.e., the ratio of magnetic field production to Ohmic dissipation, ${Rm}={v}_{\mathrm{rms}}d/\eta$, where ${v}_{\mathrm{rms}},d$ and η denote the rms flow velocity, the thickness of the shell and the magnetic diffusivity, respectively, exceeds a critical value Rm ≳ 50 (Christensen & Aubert 2006). With simple scaling arguments (e.g., Chabrier et al. 2007 and Section 6.3), it is easily verified that this condition is well fulfilled at the ionization boundary, located around P ≈ 1 Mbar, where the density is $\rho \,\approx 800\,\mathrm{kg}\,{{\rm{m}}}^{-3}$, the radius R ≈ 0.85 × R_J, and that deeper in Jupiter convective, metallic zone we obtain Rm ≫ 10⁵. This suggests that the primary dipole-dominated magnetic field is created at a depth where Rm is significant, the electrical conductivity is high, and the density contrast is relatively mild. This is indeed what is found in state-of-the-art numerical simulations that reveal that the combination of a deep-seated dipolar dynamo and a magnetic banding associated with the equatorial jet reproduce the field geometry of Jupiter with a realistic relative axial dipole, an equatorial dipole, and quadrupole and octupole field contributions (Gastine et al. 2014). These simulations are also consistent with the suggestion that the mean internal field strength as well as the mean convective velocity scale with the available convective power (Christensen & Aubert 2006). Gastine et al. (2014) find that Jupiter's surface magnetic field strength, B_J ≃ 7 G, is consistent with a typical rms flow velocity ∼3 cm s⁻¹ for a shell thickness extending from 0.2 to 0.99 R_J. Such a velocity is largely consistent with the maximum value derived in Section 6.3 around the metallization boundary. Although a dedicated study is necessary to explore this issue in the presence of an outer layered convection region, the rms velocity and the average conductivity should remain high enough for Rm to still exceed the critical value R_m ≈ 50, and thus for the reservoir of convective power to still contribute appreciably to the dynamo action.

Defining R₅₀ as the radius in the planet above which Rm ≲ 50, Duarte et al. (2018) find that R₅₀ ≃ 0.9 R_J, due essentially to the great change in conductivity when molecular H₂ fully recombines, while values below R₅₀ ≃ 0.85 R_J seem to be excluded with some confidence. This is consistent with our favorite models (Figure 6) without inclusion of the region of compositional change.

Interestingly enough, recent observations of the hemispheric difference in Jupiter's field geometry (Moore et al. 2018) are consistent with the superposition of two types of dynamo action, namely a thick-shell dynamo that is reponsible for the strong axial dipole and occurs at depth in the metallic region, and a thin-shell dynamo that yields the observed hemispheric asymmetry and occurs farther up in the envelope. A region of density gradient/layering between these two regions provides a very plausible explanation for such a separation.

Note that the large-scale field generation also constrains the size and/or the magnetic/electric properties of the diluted core. When this latter is unable to sustain dynamo action, this implies that the inner convective envelope must be sufficiently large to generate the magnetic field; this limits the maximum extension of the diluted core. This in turn limits the maximum mass of the central dense core. Indeed, as shown in the previous sections, the larger the compact core, the smaller (in absolute value) the high-order gravitational moments, but the larger the diluted core. If, however, the conductivity in the diluted core is high enough to generate electric currents, flow motions generated by density contrasts (due, e.g., to ill-mixed elements) and the Coriolis force resulting from Jupiter's spin might be able to produce magnetic fields and sustain a geodynamo process similar to the one taking place near the terrestrial iron core. In that case, the diluted core might contribute substantially to the field generation, extending the available domain for global dynamo action. We realize that at this stage, such a discussion is purely speculative, but we hope to motivate dedicated explorations of these issues as the magnetic field generation on Jupiter can certainly help to assess the reliability of our structure models.

Constraints due to Jupiter's evolution have been briefly examined in Section 4.2. Our favored models suggest an entropy jump ΔS ≈ 0.3–1 k_B/proton between the outer molecular envelope and the inner metallic envelope (see Figure 8), which yields a warmer inner temperature for the planet than in the absence of ΔS. This temperature difference, due to the entropy gradient in the inhomogeneous region, will continue to increase with time as either layered convection and/or H/He phase separation and helium sedimentation will continue to progress. This yields a heat release from inside out during the planetary evolution. The luminosity observed in Jupiter today, however, suggests that if H/He phase separation does occur in the planet, it must contribute only modestly to its cooling history. This condition can be fulfilled for several reasons. (1) If a significant fraction of this energy is devoted to heating up the interior (keeping the inner convective envelope on a high isentrope), the energy release, regardless of its source, remains modest throughout cooling, not to mention the fact that H/He might become miscible again. (2) H/He separation in Jupiter may have started only recently, contributing negligibly to Jupiter's luminosity (power), regardless of the He sedimentation energy release. (3) The H/He immiscible region, yielding a temperature gradient, might encompass a relatively modest fraction of Jupiter's interior. (4) More interestingly in the context of our favored models, core erosion, if it occurs, implies that a fraction of the internal energy of the planet must be consumed in the redistribution of heavy elements against gravity, and thus be transformed in potential energy (see, e.g., Soubiran & Militzer 1985; Guillot et al. 2004). This consumption of Jupiter's available internal energy will speed up the cooling of the planet. Even in the presence of layered convection, the final energetic balance might eventually decrease or increase the planetary cooling rate (Leconte & Chabrier 2013).

In summary, if our final models are representative of Jupiter's current internal structure and composition, its cooling history should include (i) layered convection, (ii) H/He (or other elements) phase separation plus helium dropplet sedimentation, and (iii) core erosion. Determining the resulting impact of these three processess upon the global cooling history of the planet is a highly non-trivial task that can hardly be intuited or inferred with simplistic models.

In this section, we raise another issue regarding the Galileo constraints. The Galileo measurements are taken from 1 to about 25 bar, and in all existing models, including those derived in the previous sections, the temperature profile is assumed to follow an isentrope, starting from the observed value 165 K at 1 bar (for reversible processes such as convection, an adiabat is equivalent to an isentrope, ${dQ}={TdS}=0$). It is not obvious, however, that the deeper profile (e.g., between 1 and 100 bar) does follow an isentrope. The measurements of Galileo show an increase in heavy element abundance with depth, indicating that at P ≃ 25 bar, the probe has not yet reached a well-mixed region. In case of departure from adiabaticity, the outermost temperature gradient could then be larger than the isentropic gradient, implying that the real inner entropy profile lies on a warmer isentrope than the profile that is obtained when one assumes that it is given by the P = 1 bar, T = 165 K condition.

We consider, notably, the impact of rotation. There is currently no well-defined theory for turbulent convection in the presence of rotation, so that we can only rely on estimates. At the pressure level P ≃ 10 bar in the external envelope of Jupiter, the optical depth is τ ≫ 100, so that except for the possible impact of rotation, one can safely assume that the profile is isentropic at this level. This pressure corresponds to ρ ∼ 1 kg m⁻³, $g={{GM}}_{J}/{R}^{2}\simeq 20$ m s⁻², T ∼ 200 K, and thus a typical convective length $l\sim {H}_{P}=P/(\rho g)\approx {10}^{5}$ m, about 1% of the planetary radius. Assuming that the entire internal flux of Jupiter, ${ \mathcal F }\simeq 5.4$ W m⁻², is transported by convection and using the usual equations of the mixing length theory (Kippenhahn & Weigert 1990), this yields a typical superadiabatic gradient ${\rm{\Delta }}{\rm{\nabla }}{T}_{0}={({\rm{\nabla }}-{{\rm{\nabla }}}_{\mathrm{ad}})}_{0}\approx {10}^{-6}$ at this pressure level, i.e., a characteristic convective velocity v_conv ≲ 1 m s⁻¹. Because the angular velocity of Jupiter is ${{\rm{\Omega }}}_{J}\,={v}_{\mathrm{rot}}/{R}_{J}=1.76\times {10}^{-4}$ rad s⁻¹, the ratio of inertial to Coriolis forces, known as the Rossby number, at 10 bar is thus ${Ro}={v}_{\mathrm{conv}}/(l{\rm{\Omega }})\gtrsim 0.1$. Convection at the top of the upper envelope, where the Galileo measurements have been made, should thus be only moderately affected by rotation. It is also easily verified that the Coriolis acceleration is much smaller than the gravity, R_JΩ² ≪ g, which allows performing a linear stability analysis of the MLT equations in the presence of a Coriolis force, $2{\rm{\Omega }}\times {v}_{\mathrm{conv}}$ (Chandrasekhar 1961). Additionally, given the value of the Rossby number, this linear analysis can be performed in the strong rotation limit (Stevenson 1979). This yields for the suradiabtic gradient in the presence of rotation in the region probed by Galileo (Stevenson 1979)

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\nabla }}{T}_{\simeq }{\left({\rm{\Delta }}{\rm{\nabla }}{T}_{0}\right)}^{3/5}{\left(\displaystyle \frac{{{\rm{\Omega }}}^{2}l}{g}\right)}^{2/5}\simeq 6\times {10}^{-6}.\end{eqnarray} \tag{ 13 }$$

This estimate shows that rotation cannot yield a strong departure from adiabaticity in the outermost envelope layers of Jupiter, as expected from the inferred Rossby number value, in contrast to deep convective regions (see, e.g., Chabrier et al. 2007).

Moreover, at a pressure of 1 bar, the atmosphere is composed of alternative superrotating and underrotating stripes in latitude. Both at the equator, where advection dominates, and in the mid-to-high latitudes, where geostrophy applies, one can show that the horizontal variation of temperature due to the winds is on the order of a few percents in latitude and longitude, with a maximum at the equator (as confirmed by the observations of Fisher et al. 2016 and the temperature profiles from GCM calculations of Schneider & Liu 2009). In that regard, the measurements of Galileo are more an upper bound than a lower bound on the temperature, and deviations from these measurements are small.

In conclusion, it seems quite safe to assume that the external condition defined by the T = 165 K, P = 1 bar condition lies on an adiabatic profile.