<figure><img src="image_1.png" /></figure>

# New Models of Jupiter in the Context of *Juno* and *Galileo*

Florian Debras<fnref n="1" />,<fnref n="2" /> ID and Gilles Chabrier<fnref n="1" />,<fnref n="2" /> ID

Received 2018 October 5; revised 2018 December 27; accepted 2018 December 28; published 2019 February 14

## Abstract

Observations of Jupiter's gravity field by *Juno* have revealed surprisingly low values for the high-order gravitational moments, considering the abundances of heavy elements measured by *Galileo* 20 years ago. The derivation of recent equations of state for hydrogen and helium, which are much denser in the megabar region, exacerbates the conflict between these two observations. In order to circumvent this puzzle, current Jupiter model studies either ignore the constraint from *Galileo* or invoke an ad hoc modification of the equations of state. In this paper, we derive Jupiter models that satisfy constraints of both *Juno* and *Galileo*. We confirm that Jupiter's structure must encompass at least four different regions: an outer convective envelope, a region of compositional and thus entropy change, an inner convective envelope, an extended diluted core enriched in heavy elements, and potentially a central compact core. We show that in order to reproduce *Juno* and *Galileo* observations, one needs a significant entropy increase between the outer and inner envelopes and a lower density than for an isentropic profile, which is associated with some external differential rotation. The best way to fulfill this latter condition is an inward-decreasing abundance of heavy elements in this region. We examine in detail the three physical mechanisms that can yield such a change of entropy and composition: a first-order molecular-metallic hydrogen transition, immiscibility between hydrogen and helium, or a region of layered convection. Given our present knowledge of hydrogen pressure ionization, a combination of the two latter mechanisms seems to be the most favored solution.

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

## 1. Introduction

For more than 30 years, guided by the observations of *Voyager* and *Pioneer* (Campbell & Synott 1985), all traditional models of Jupiter have been described as two- or three-layer models, namely a homogeneous, convective gas-rich envelope, generally split into a molecular/atomic outer part and an ionized/metallic inner part, and a core that is assumed to be solid (e.g., Chabrier et al. 1992; Saumon & Guillot 2004), as first intuited by Stevenson & Salpeter (1977a, 1977b). Later on, *Galileo* provided new constraints on Jupiter's outer layer composition (von Zahn et al. 1998 for helium and Wong et al. (2004) for the reanalyzed results of the heavy elements). Finally, in 2017 and 2018, the observations of *Juno* reported in Bolton et al. (2017) and Iess et al. (2018) stressed the need to resolve a real puzzle: how can an internal structure of Jupiter be modeled that matches the observations for the gravitational moments of *Juno* and for composition of *Galileo*, which latter reveal a highly supersolar outer element abundance?

The trouble indeed is to reconcile the low value of *Juno*'s high-order even gravitational moments, $J_4$ to $J_{10}$, and the high value of helium and heavy elements observed by *Galileo*, $Y_{Gal}$ and $Z_{Gal}$. The higher the order of a gravitational moment, the more sensitive it is to the outermost part of the planet. Hence, the most important physical parameters for determining the values of $J_4$ to $J_{10}$, for a given mass and $J_2$, are the abundances of helium and heavy elements in the external envelope of the planet.

In order to resolve this puzzle, Wahl et al. (2017) either invoked an ad hoc modification of their H/He equation of state (EOS) or reduced the outer heavy element content compared with *Galileo*'s observations. Guillot et al. (2018) also allowed

the outer heavy element content to vary from 0 to $Z_{Gal}$, but their model that matches all *Juno* $J_n$ values has an amount of heavy elements in the atmosphere that is not compatible with the *Galileo* constraints (T. Guillot 2018, private communication).

In this paper, we present models of Jupiter that do fulfill the observational constraints of both *Juno* and *Galileo*. We describe the method and the different physics inputs in Section 2. In Section 3 we demonstrate the necessity of having several different regions in Jupiter's interior and show that traditional two- or three-layer models fail to reproduce the observations. In Section 4 we show that a locally inward-decreasing abundance of heavy elements in the megabar (Mbar) region is the favored solution for resolving this puzzle. We explore the possibility of such an element distribution in detail.

Our final models are presented in Section 5. We first show that without a sharp entropy increase somewhere within the gaseous envelope, the values of $J_6$ to $J_{10}$ are too high compared with those of *Juno*, which implies that an implausibly large amount of differential rotation needs to be invoked. Indeed, a strong entropy increase in the region of hydrogen metallization (around 1 Mbar) yields higher internal temperatures, allowing a larger amount of heavy elements in the central region (Section 5). This in turn affects the high-order gravitational moments and enables us to derive Jupiter models that satisfy the observational constraints of both *Juno* and *Galileo*. We examine the possible physical mechanisms in detail that lead to this type of internal structure and discuss their implications for the physics of hydrogen pressure ionization. We also examine the possible amount of differential rotation in Jupiter. In Section 6 we summarize and examine the validity of the main physical assumptions that are made throughout this study. Section 7 is devoted to the conclusion.

<footnote n="1">¹ Ecole normale supérieure de Lyon, CRAL, UMR CNRS 5574, F-69364 Lyon Cedex 07, France; florian_debras@hotmail.com</footnote>

<footnote n="2">² School of Physics, University of Exeter, Exeter, EX4 4QL, UK</footnote>
<table><caption>**Table 1**
Values of the Planetary Parameters of Jupiter</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>G (a) (global parameter)</td><td>6.672598 × 10<sup>-11</sup> ± 2 × 10<sup>-17</sup> m<sup>3</sup> kg<sup>-1</sup> s<sup>-2</sup></td></tr><tr><td>G × M<sub>J</sub>(b)</td><td>(126686533 ± 2) × 10<sup>9</sup> m<sup>3</sup> s<sup>-2</sup></td></tr><tr><td>M<sub>J</sub></td><td>1.89861 × 10<sup>27</sup> kg</td></tr><tr><td>R<sub>eq</sub> (c)</td><td>71492 ± 4 km</td></tr><tr><td>R<sub>polar</sub> (c)</td><td>66854 ± 10 km</td></tr><tr><td>ω (d)</td><td>1.7585324 × 10<sup>-4</sup> ± 6 × 10<sup>-10</sup> s<sup>-1</sup></td></tr><tr><td>ρ&#773;</td><td>1326.5 kg m<sup>-3</sup></td></tr><tr><td>m = 3ω<sup>2</sup>/4πGρ&#773;</td><td>0.083408</td></tr><tr><td>q = ω<sup>2</sup>R<sub>eq</sub><sup>3</sup>/GM<sub>J</sub></td><td>0.0891954</td></tr><tr><td>J<sub>2</sub> × 10<sup>6</sup> (e)</td><td>14696.572 ± 0.014</td></tr><tr><td>-J<sub>4</sub> × 10<sup>6</sup> (e)</td><td>586.609 ± 0.004</td></tr><tr><td>J<sub>6</sub> × 10<sup>6</sup> (e)</td><td>34.198 ± 0.009</td></tr><tr><td>-J<sub>8</sub> × 10<sup>6</sup> (e)</td><td>2.426 ± 0.025</td></tr><tr><td>J<sub>10</sub> × 10<sup>6</sup> (e)</td><td>0.172 ± 0.069</td></tr></tbody></table>

**Note.** R<sub>eq</sub> and R<sub>polar</sub> 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).

## 2. Method

### 2.1. Concentric MacLaurin Spheroids

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, $\beta = 6/N$, and the upper atmosphere is neglected.<fnref n="3" /> 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 $\beta$ 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.

### 2.2. Equations of State

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{aligned} \frac{1}{\rho_{\text{MH13}}} &= \frac{X_{\text{MH13}}}{\rho_{\text{H}}} + \frac{Y_{\text{MH13}}}{\rho_{\text{He,New}}} &\Rightarrow \frac{1}{\rho_{\text{H}}} \\ &= \frac{1}{X_{\text{MH13}}} \left( \frac{1}{\rho_{\text{MH13}}} - \frac{Y_{\text{MH13}}}{\rho_{\text{He,New}}} \right), \end{aligned} \quad (1) $$

$$ \begin{aligned} S_{\text{MH13}} &= X_{\text{MH13}}S_{\text{H}} + Y_{\text{MH13}}S_{\text{He,New}} + S_{\text{mix}} \Rightarrow S_{\text{H}} \\ &= \frac{1}{X_{\text{MH13}}}(S_{\text{MH13}} - Y_{\text{MH13}}S_{\text{He,New}} - S_{\text{mix}}), \end{aligned} \quad (2) $$

where $\rho_H$ is the mass density for pure hydrogen, $\rho_{\text{MH13}}$ is the density derived from MH13 by spline procedures, $\rho_{\text{He,New}}$ is the helium density in the Chabrier et al. (2019) EOS, $S_H$ is the sought pure hydrogen specific entropy, $S_{\text{MH13}}$ is the splined specific entropy from MH13, $S_{\text{He,New}}$ is the helium-specific entropy in the new EOS, all at the same (P,T), and $X_{\text{MH13}} = 0.7534$, $Y_{\text{MH13}} = 0.2466$ are the mass fractions of hydrogen and helium in the MH13 simulations. Finally, $S_{\text{mix}}$ is the mixing-specific entropy defined as

$$ \frac{S_{\text{mix}}}{k_b} = \frac{1}{M_{\text{H,He}}} \left[ N_{\text{H}} \ln \left( 1 + \frac{N_{\text{He}}}{N_{\text{H}}} \right) + N_{\text{He}} \ln \left( 1 + \frac{N_{\text{H}}}{N_{\text{He}}} \right) \right], \quad (3) $$

with $N_H$ and $N_{\text{He}}$ the numbers of H and He particles, respectively, of number fractions $x_i = N_i/(N_H + H_{\text{He}})$ ($i \equiv H$ or He), $M_{\text{H,He}} = \bar{A}m_H$ the total mass of H+He, $\bar{A} = \sum_i x_i A_i$ the mean atomic number, and $m_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_{\text{MH13}})/\rho_{\text{MH13}}$, 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 kg m<sup>-3</sup>, 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.

<footnote n="3">³ This implies an irreducible error on the order of $10^{-7}$ on $J_2$ and a few $10^{-8}$ on higher-order moments, which is negligible compared to the possible impact of differential rotation (Kaspi et al. 2017; Debras & Chabrier 2018).</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. Relative error on the density between MH13 and the mixed EOS used by us or Miguel et al. (2016) for Y = 0.2466.</figcaption></figure>

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 H/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:

$$ \frac{1}{\rho} = \frac{1 - Z_{\text{water}} - Z_{\text{drysand}}}{\rho_{\text{H,He}}} + \frac{Z_{\text{water}}}{\rho_{\text{water}}} + \frac{Z_{\text{drysand}}}{\rho_{\text{drysand}}} \quad \text{at } (P, T) = \text{constant,} \quad (4) $$

where $\rho_{H,He}$, $\rho_{water}$ and $\rho_{drysand}$ are the densities of the H/He mixture, water and drysand, respectively, and $Z_{water} = M_{water}/M$, $Z_{drysand} = M_{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 \gtrsim 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^{\text{eff}} = X/(1-Z)$ and $Y^{\text{eff}} = Y/(1-Z)$, respectively, with $X^{\text{eff}} + Y^{\text{eff}} = 1$, and $X = M_H/M$, $Y = M_{He}/M$, and $Z = M_Z/M = Z_{water} + Z_{drysand}$.

## 2.3. Galileo Constraints on the Composition

For the outer element abundances, the observations of Galileo give

$$ \frac{Y_{\text{Gal}}}{(X_{\text{Gal}} + Y_{\text{Gal}})} = 0.238 \pm 0.005, \\ Z_{\text{Gal}} = 0.0167 \pm 0.006, $$

where $X_{Gal}$ and $Y_{Gal}$ are the observed mass abundances of hydrogen and helium, respectively. $Z_{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_{Gal}$ and $X_{Gal}$ are only defined relative to each other and that $X_{Gal} + Y_{Gal} + Z_{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

$$ Y_{\text{ext}} = 0.23 \\ Z_{\text{ext}} = 0.02, $$

which corresponds to $Y_{ext}/(X_{ext} + Y_{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
<figure><img src="image_3.png" /><figcaption>**Figure 2.** $J_4 \times 10^6$ vs. $J_6 \times 10^6$ for different models (see text) and *Juno* values (Iess et al. 2018); the error bars too small to be visible in this figure. The shaded area correspond to the uncertainty on the gravitational moments arising from differential rotation shallower than 10,000 km, as evaluated by Kaspi et al. (2017).</figcaption></figure>

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

# 3. Simple Benchmark Models

In this section, we show that traditional homogeneous, adiabatic two- or three-layer interior models for Jupiter are excluded by the new observations and that the planet must consist of several different regions.

## 3.1. Homogeneous Adiabatic Gaseous Envelope

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_2$ observed value within $10^{-7}$, the maximum intrinsic precision of the CMS method (Debras & Chabrier 2018). The $J_4$ and $J_6$ 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.

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_4$ 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.

## 3.2. Region of Compositional and Entropy Variation within the Planet

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

A diluted core or a region of layered convection farther up in the planet can emerge during the evolution of the planet or can be inherited from the formation process (see, e.g., Soubiran & Militzer 1985; Chabrier & Baraffe 2007 and references therein, Helled & Stevenson 2017) and is characterized by a compositional gradient $\nabla\mu = (\partial\mu/\partial r) \propto \nabla Z = (\partial Z/\partial r)$, with $\mu$ the mean molecular weight, thus an entropy gradient $\nabla S = (\partial S/\partial r)$. For hydrogen metallization, first-principle numerical simulations suggest that it might occur through a first-order phase transition (usually denominated plasma phase transition, PPT, or liquid-liquid transition, LLT) in a domain $P_c \sim 1-2$ Mbar for temperatures below the critical temperature $T \le T_c \approx 2000-5000$ K (Morales et al. 2010, 2013b; Lorenzen et al. 2011; Knudson et al. 2015; Mazzola et al. 2018). Experiments on liquid deuterium, $D_2$, seem to be consistent with these figures, even though significant differences still persist between various experiments (see, e.g., Knudson et al. 2015; Celliers et al. 2018). Until this issue is resolved definitively, an entropy discontinuity, $\Delta S$ due to hydrogen pressure ionization in Jupiter's envelope, although unlikely, can thus not be definitely excluded. In a similar vein, H/He phase separation, also a first-order transition, will also yield an entropy discontinuity if the local temperature is lower than the critical temperature for the appropriate He concentration (see below). Last but not least, a regime of double-diffusive layered convection could develop somewhere within the planet interior, triggered either by one of these two transitions (or by any phase separation involving some heavy component insoluble in metallic hydrogen) and/or simply by a local compositional gradient (Leconte & Chabrier 2012). Phase transitions, indeed, notably endothermic ones, are suspected of enforcing layered convection, for instance in the Earth's mantle (Christensen & Yuen 1985). The physical reason is the release of latent heat at the transition, which leads to thermal expansion and temperature advection, which tends to hamper convection. It is interesting to note that due essentially to the higher entropy in the plasma phase than in the molecular phase, a PPT is an endothermic transition, i.e., $dP/dT < 0$ along the critical line of the transition, according to the Clausius-Clapeyron equation, $dP/dT \propto -\Delta S/\Delta\rho$. As for the H/He immiscibility, ab initio simulations, while still differing substantially, seem to suggest a critical temperature for $x_{He} = 0.08$ in the range $\sim 2000-8000$ K for $P \gtrsim 1$ Mbar (100 GPa), with a weak dependence upon pressure in the $T, P$ domain relevant for Jupiter, suggesting $dP/dT \sim 0$ for the protosolar helium value (Lorenzen et al. 2009; Morales et al. 2009, 2013a; see Figure 13 below).

Therefore, the entropy variation in the gaseous envelope could occur either within a region of layered convection due to compositional gradients or because of either a PPT or a H/He phase separation. Needless to say, not only are these three physical processes not exclusive, but they are likely to be
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Simple structures of Jupiter with an internal entropy variation either in the envelope (left) or in the core (right) (see Section 3.2.2). None of these models can match the observations of both *Juno* and *Galileo*.</figcaption></figure>

tightly linked and thus to take place in the same more or less extended region that we hereafter denominate “metallization boundary region” near the Mbar.

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

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*₄| and |*J*₆| 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 observa-
tional 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.

4. Locally Inward-decreasing Z-abundance in the Gaseous Envelope

4.1. Inward-decreasing Abundance of Heavy Elements in Some Part of the Outer Envelope

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 $\nabla 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

Salpeter (1973) and Stevenson & Salpeter (1977a; see also
Stevenson 1979) suggested the occurrence of helium differ-
entiation in giant planet interiors, in either the same or a
different region than the H<sub>2</sub>-H<sup>+</sup> pressure metallization of
hydrogen. These authors suggested that minor constituents,
namely the heavy elements, could suffer differentiation in a
way similar to or even larger than helium. Unfortunately, phase
diagram calculations of two or more components under the
typical relevant conditions for Jupiter (about 1 Mbar and
5000 K) are scarce or even nonexistent, which means that
determining which element under which molar concentration
prefers the H<sup>+</sup>-rich or He-rich phase is an open question. The
only existing study has been published by Wilson & Militzer
(2010). Ab initio simulations by these authors suggest that Ne
association with He is thermodynamically favored, while the
opposite is true for Ar, which is found to be more soluble with
H<sup>+</sup>. The underlying physical reason is the additional electron
shell of the argon atom, which increases its effective volume
with respect to He due to the Pauli exclusion principle. If this
explanation is correct, Kr and Xe should likewise be soluble in
metallic hydrogen, which is consistent with their observed
nondepletion in Jupiter's atmosphere.

It is indeed intuitively appealing to think that in case some
species, $Z_i$, is pressure ionized, it might become immiscible
with neutral helium, as for H<sup>+</sup>/He, due to the strongly
repulsive pseudo-potential, as in the case, e.g., of alkali metals
(Stevenson 1979). For some element (atom or molecule) to
differentiate in the midst of some mixture, one needs its
interaction energy in the mixture, typically the molecule or
electron binding energy, to be higher (in absolute value) than
the ideal mixing entropy, -k<sub>B</sub> T ln x<sub>Z</sub>. Since the most abundant

heavy elements have a number fraction $x_Z \approx 0.1\%$, this yields
near the metallization boundary, ~5000 K, $|E_b| \gtrsim 0.5$ eV, a
condition that is rather easy to fulfill. As mentioned above, all
heavy elements, however, do not necessarily behave similarly.
Heavy noble gases, indeed, are more likely to form compounds
(Hyman 1964; Blackburn 1966; Wilson & Militzer 2010),
suggesting that species such as neon, which act like helium,
and argon, behave differently in the H<sup>+</sup>/He mixture.

In case of element differentiation, according to the Gibbs
phase rule, $x_{II} = x_I \exp\{-\Delta G(P, T)/k_B T\}$, where $x_i$ denotes
the number abundance of a given species in phase I or II
(H<sup>+</sup>-rich/poor, conversely, He-poor/rich in the present con-
text) and ΔG is the excess mixing enthalpy in the mixture, the
differentiation of a given heavy element can be similar or
opposite to that of helium, yielding an increasing (decreasing)
abundance with depth in the former (latter) case. In all cases,
this yields a gradient of abundance $\nabla Z$ within some part of the
planet envelope, with $\nabla Z > 0$ for some heavy elements.

If H-He immiscibility, leading to a depletion of helium in
the outer envelope, is triggered by the metallization of
hydrogen, the fact that hydrogen metallization in a H/He
mixture is found to occur at lower pressures with decreasing
helium fraction (e.g., Mazzola et al. 2018) implies that the
pressure range of immiscibility will extend with time because
of the planet's decreasing internal temperature and the
decreasing abundance (depletion) of helium in the upper
layers. The region of immiscibility can thus be relatively broad
in Jupiter's interior, depending on when it started. We note in
passing that if H/He differentiation occurs and is at least partly
responsible for the redistribution of heavy elements in Jupiter's
envelope, this excludes the H/He diagrams suggested by
Morales et al. (2009) and Schöttler & Redmer (2018), which
predict no immiscibility within the current Jupiter. This point is
discussed in detail in Section 6.1.

It should be stressed that given the small number abundances
of helium and heavy elements in the Jupiter interior, whereas
the aforedescribed demixing processes could lead to some
Z-enrichment in the outer envelope of the planet, consistent with
or even higher than the observed *Galileo* value, this enrichment
will remain modest and could not explain values significantly
higher than $Z_{\text{Gal}}$. In this latter case, external accretion seems
necessary, as examined in Section 4.1.3 below.

4.1.2. Upward Atomic Motions

As explored thoroughly by Stevenson & Salpeter (1977a) in
the case of a first-order hydrogen pressure ionization (PPT) and
a H/He phase separation occurring in the same region, the
following process might occur. At the onset of hydrogen
metallization, characterized by a pressure P+, latent heat
release will lead to the superposition of an overheated
(supercooled) H+-rich (H₂-rich) layer, thus less (more) dense
than the surrounding medium, underneath (above) the metalli-
zation boundary (see Figure 2 of Stevenson & Salpeter 1977a).
Under such conditions, nucleation of bubbles might occur.
Concomitantly, He atoms will differentiate from H+. If such
H+-rich/He-poor bubbles form, they will absorb heat by
thermal diffusivity and will be lighter than the surrounding gas.
The bubbles will then rise by buoyancy, up to a pressure P
lower than the metallization pressure, P < P+. They will then
break, and H+ will recombine to form H₂, depleting the upper
envelope in He little by little by mixing this convective region
with H+-rich bubbles while enriching the lower envelope in
He. Consequently, the heavy elements, which for chemical and/or thermodynamic reasons have a preference for these H⁺-rich/He-poor bubbles rather than for the He-rich/H⁺-poor surrounding medium, will be transported upward and be depleted little by little in the deep envelope, whereas the opposite will be true for species favoring association with helium atoms. This is similar to an ongoing distillation process in the sense that the redistribution of elements arises from a physical separation rather than a chemical reaction and mass is not locally conserved beneath the uppermost convective envelope. This occurs only if the heavy elements do not significantly affect the density of the bubbles, which must remain lighter than the surrounding gas. In the typical conditions of Jupiter's outer envelope, there are about 500 times more atoms of hydrogen than of heavy elements. For a typical atomic weight ratio $\dot{A}_Z/A_H \sim 15/1$ (average between C, N, and O atoms), such a process is therefore possible. Heavier molecules (such as iron) being even rarer than H, gravitational considerations are still consistent with this scenario.

Such a scenario has further theoretical support. First, noble gases have been known to be almost insoluble in metals since the end of the nineteenth century (Ramsay & Travers 1897 or Blackburn 1966 for a review, and Wilson & Militzer 2010 for the case of neon and argon). On the other hand, at high pressure, hydrogen can very efficiently form complex polyhydrides molecules, with many different atoms (sulfur, lithium, sodium, iron, etc.; see, e.g., Ashcroft 2004 or Pépin et al. 2017). Therefore, at metallization, non-inert heavy elements tend to form polyhydrides within metallic H⁺-rich bubbles. If, as discussed above, the density of these bubbles is lower than the density of the surrounding gas, these heavy elements will be transported upward, enriching Jupiter's outer envelope while depleting the inner envelope. The formation of polyhydrides, however, has been probed experimentally so far up to ~1500 K (Pépin et al. 2017) and remains to be explored up to T ~ 4500 K, the onset of H metallization in Jupiter. Further numerical or experimental work on the formation of polyhydrides at high pressures and temperatures would help to assess the validity of this process.

Concomitantly with hydrogen metallization and the formation of polyhydrides, we also expect reduction-oxidation (redox) reactions to occur. The loss of its 1 s electron at hydrogen metallization makes H⁺ prone to react with other heavy elements through electron transfer. The H⁺ bubbles could then trap N, O, or other elements, for instance, participating also in an enrichment (depletion) of these elements in the upper (lower) envelope. We recall, however, that in the absence of dynamical variations such as gravity waves or upward plume penetrations, the amount of overheating due to a PPT would be insufficient to yield homogeneous nucleation (Stevenson & Salpeter 1977a).

In all cases, the “rising bubble” process described above requires hydrogen molecular-metallic transition to occur through a first-order transition, leading to a local release of latent heat. As mentioned in Section 3.2.1, although a PPT is indeed found in some modern ab initio numerical simulations, its critical temperature remains to be determined precisely, being predicted within the range 2000 $\lesssim T_c \lesssim 5000$ K for a critical pressure $1 \lesssim P_c \lesssim 2$ Mbar. Interestingly enough, a critical temperature $T_c \approx 5000$ K around ~1 Mbar would be consistent with a PPT in the outer part of Jupiter. In the absence of a PPT in the envelope, H⁺/He phase separation might still

occur, but the outer observed oversolar abundance of heavy elements is ruled out due to upward bubble motions. Phase separation of these elements with He, as discussed in Section 4.1.1, will thus be the favored explanation.

### 4.1.3. Accretion

Finally, the overabundance of heavy elements in Jupiter's upper envelope can have a third explanation, namely one or several giant impacts (Iaroslavitz & Podolak 2007) or, similarly, ongoing accretion of planetesimals (e.g., Bézard et al. 2002). This scenario, however, implies that global internal convective motions must be inhibited somewhere in Jupiter, preventing the extra accreted material from being redistributed homogeneously throughout the planet. Indeed, a global Z-abundance throughout the planet equal to the Galileo value would yield low-order gravitational moments that are inconsistent with observations. If convection inhibition is due to H/He immiscibility, this latter must already have started when the external event took place. This in turn places an important constraint on the H/He phase diagram, notably on the critical P, T values for $x_{He} = 0.08$, the helium protosolar concentration in Jupiter. If inhibition is due to hydrogen metallization, it implies that this latter must be a first-order phase transition (yielding an entropy jump). Convection can also be inhibited by the onset of double-diffusive convection, either as an enhanced diffusive process (oscillatory convection), or as layered convection, a process that might be triggered by extensive planetesimal accretion (e.g., Soubiran & Militzer 1985; Chabrier & Baraffe 2007 and references therein) and/or by deposition of high-entropy material onto the growing planet (Berardo & Cumming 2017), which prevents homogeneization of the envelope composition. In order to explain a genuine abundance significantly higher than the observed Galileo value, $Z_{ext} \approx 2.5 \times Z_{\odot}$, the total accreted mass must be $M_{excess} \lesssim 1.5 M_{\oplus}$, a significant but not implausible value.

## 4.2. Constraints from the Evolution

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 $\nabla 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.

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),
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Typical evolution of the heavy element content Z as a function of the radial distance r with time, assuming that a first-order phase transition (PPT) or immiscibility will occur during the cooling history. At $t = t_0$, Jupiter just formed. A small convective external envelope is connected to a gradually diluted planet, structure inherited from the core accretion. At $t = t_c$, immiscibility or first-order metallization is about to occur. The convective zone has somewhat expanded, redistributing the metal content in the planet. Later on ("metallization" curve), the rising bubbles have enriched the outer envelope in heavy elements and depleted the inner envelope. If immiscibility takes place ("immiscibility" curve), the outer envelope will eventually lose part of its mass because of drowning helium droplets, increasing the heavy element mass fraction in the outer envelope, whereas the dilution of these elements in the inner envelope is almost negligible.</figcaption></figure>

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 $\nabla Z > 0$ between the outer Z-enriched and the inner Z-depleted envelopes, the He-rich falling droplets 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 droplets, 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, $\nabla Z < 0$, which is at odds with our suggestion of a local $\nabla 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 Jupiter. As just mentioned, however, properly handling such complex physical processes is a task of major difficulty (see Section 6.5).

## 5. Models with at Least Four Layers and an Entropy Discontinuity in the Gaseous Envelope

As shown in Section 3, the interior structure of Jupiter must entail at least four different regions, namely two outer and inner homogeneous adiabatic envelopes, separated by a region of compositional and thus entropy variation, a diluted core that also harbors a more or less extended domain of compositional/entropy change, and potentially a solid rocky seed. One of the unknowns in these models is the amount of entropy change in the envelope. With an entropy change that is only due to a change in composition (see Section 2), the lowest value of $J_4$ we could obtain lies within the limit of what can be explained by a differential rotation shallower than 10,000 km (Kaspi et al. 2017). The higher-order moments, however, remain far too high. This yields another conclusion:

1. Models with a small entropy change in the Mbar region seem to be excluded as a possible internal structure of Jupiter.

In this section the inward increase of entropy (due to H-He immiscibility or the onset of superadiabatic layered convection) is therefore now used as a free parameter in the calculations, and is discretized throughout a certain number of spheroids across the ionization boundary region. That is, we assume an entropy gradient $\nabla S = \Delta S/\Delta R < 0$ within the relevant pressure range.

### 5.1. Physical Expectations

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 \approx 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 \approx 2000-5000$ K at $P \approx 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, $\Delta S$ is given by the Clausius-Clapeyron relation along the critical line $P(T)$

$$ \frac{\Delta S}{\Delta(1/\rho)} = \frac{dP}{dT}. \qquad (5) $$
Analytical calculations (Saumon & Chabrier 1992) suggest

$$ \Delta S_{\text{metallization}} \sim 0.5 k_{\text{B}}/\text{proton}. \qquad (6) $$

Because, as mentioned above, the temperature in this region of Jupiter's interior should be close to $T_c$, we expect $\Delta 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 droplets 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_{\text{He}} \simeq 0.1$, suggests that the variation $\Delta 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_{\text{He}} = 0.1$. Because the maximum value of the ideal mixing entropy for a concentration $x_{\text{He}} = 0.5$ is $S_{\text{mix}}^{\text{id}}/N = -[x \ln x + (1-x)\ln(1-x)] = 0.7 k_B$/at (about 0.3 $k_B$/at for $x_{\text{He}} = 0.1$), we see that the entropy jump due to H/He immiscibility should be $\lesssim 0.5 k_B$/proton. The entropy change due to helium droplet 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 \rightarrow 0.36$ corresponds to a global increase in entropy $\Delta S \sim 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

$$ \Delta S_{\text{H/He}} \lesssim 1.0 k_{\text{B}}/\text{proton}. \qquad (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,

$$ \Delta S_{\text{layered}} \lesssim 0.6 k_{\text{B}}/\text{proton}, \qquad (8) $$

with a decreasing metal abundance, i.e., $\nabla Z > 0$, in this region. An increasing metal abundance in this region yields higher values of $\Delta S$, but in this case, $\alpha \lesssim 10^{-7}$, where $\alpha$ denotes the ratio of the size of the convective layer to the pressure scale height, $\alpha = 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 ~0.5–1 $k_B$/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 ($\nabla Z > 0$), and (ii) a non-negligible entropy ($\gtrsim 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, $\Delta S \gtrsim 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 $\Delta 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, $\Delta S \lesssim 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.

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

## 5.2. Quantitative Results on the Gravitational Moments

The results of our optimized models with an entropy discontinuity $\Delta S \in [0-2]k_B$/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_{\text{ext}} = Z_{\odot}$ and $Z_{\text{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
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Schematic internal structure of our final Jupiter models. $Y_{\text{ext}} = 0.23$, $Z_{\text{ext}} = 0.02$, as stated in the text. (a) The modest entropy increase between the outer and inner envelopes yields a moderate helium increase in this latter, and an inward helium enrichment in the diluted core (see text). (b) The larger entropy increase in the inhomogeneous boundary region yields a supersolar helium fraction in the inner envelope, but then the helium abundance decreases in the diluted core. (c) Our least favored model. An increase in both helium and heavy element abundances in the inner envelope requires a strong entropy increase, at the limit of what can physically be achieved. A mixture of structures (a) and (c) is also possible, with a small increase in both helium and heavy elements. The required $\Delta S$ would be comparable to (b).</figcaption></figure>

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.

1. *Impact on the low-order moments* ($\le J_4$). We found out that the entropy change, $\Delta 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 $\Delta 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_4$ at constant $J_2$. Figure 7(a) shows the value of the contribution function $J_2 - J_4$ in the planet as a function of pressure.

This figure shows that in order to decrease $J_4$ with respect to $J_2$, one needs to enrich the planet deeper than $\sim 2$ Mbar, and the region around $\sim 10$ Mbar is where it is most efficient. Therefore, an enriched inner envelope decreases $J_4$ at constant $J_2$ (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_4$ compared to $J_2$. Furthermore, this boundary region between $\sim 0.1$ and $2$ Mbar has a much stronger contribution on $J_2$ and $J_4$ 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_2$ and $J_4$ is enhanced. In
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Gravitational moments that we obtain with various models in the $J_k - J_{k+2}$ plans for two values of $Z_{ext}$. All the values are multiplied by $1 \times 10^6$, as was done in, e.g., Guillot et al. (2018). The green dot is the *Juno* value with the observed error bars (Iess et al. 2018), without any dynamical correction. The hatched area corresponds to the models of Guillot et al. (2018).</figcaption></figure>

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 ($\nabla Z < 0$) of the envelope can fulfill all observational constraints (case (c)) if the entropy change around ~1 Mbar reaches values $\Delta S > 1k_B$/proton. This requires a strong

entropy discontinuity induced either by a PPT for $T \simeq 4500$ K (thus a critical point $T_c \gg 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 $\nabla 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 ($\nabla 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
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Subtraction of the contribution functions for $J_2$ minus $J_4$, both normalized to their maxima, with respect to pressure or radius for an isentropic model. (b) Normalized difference of contribution functions for $J_2$ to $J_8$ between a model with a small entropy change and a model with a strong entropy change.</figcaption></figure>

**conclusion 1:** for given $J_2$, a central compact core tends to decrease the moments of order $\ge J_4$ 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 $\Delta S$ between the outer and inner envelopes. This in turn has consequences for the high-order moments: a model with a high $\Delta 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 $\Delta S$. The lower density means that this region has a smaller contribution to the gravitational moments than in the case of a small $\Delta S$. Although for $J_2$ and $J_4$, 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 $\Delta J_n$ in the contribution functions $J_2$ to $J_8$ between a model with a small entropy change, $\delta S < 0.1 \ k_B$/proton, in the ionization boundary region that is only due to a composition change, and a model with a total entropy discontinuity $\Delta S = 0.9 \ k_B$/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_2$ than on the other moments. A strong $\Delta S$ thus requires more heavy elements in the diluted core to preserve $J_2$, 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_2$. 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 $\ge 5 \ M_\oplus$, one needs $\Delta S > 1 \ k_B$/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 $\Delta S$. Low (absolute) values of $J_6$, $J_8$, or $J_{10}$ allways require a significant $\Delta S$, except when we decrease the atmospheric $Z_{\text{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.

As seen in Figure 8, none of our models can match the $3\sigma$ error bars on $J_6$ for $\Delta S < 1 \ k_B$/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_{\text{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 $\Delta J_6$ 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_B$/proton. In any case, a continuously increasing heavy element mass fraction with depth, i.e., $\nabla Z < 0$, in the Mbar region is hard to justify (in Figure 8, these models all have $\Delta S \gtrsim 1 \ k_B$/proton).

As shown in Figure 5, the valid models predict a size for the metallization boundary region, $l_b \approx 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 $\Delta S$.
<figure><img src="image_9.png" /></figure>

(a)

<figure><img src="image_10.png" /></figure>

(b)

<figure><img src="image_11.png" /></figure>

(c)

**Figure 8.** High-order gravitational moments as a function of the entropy jump $\Delta S$ in the envelope. All these models have $J_2$ and $J_4$ values within the allowed range from *Juno* observations with the Kaspi et al. (2018) correction, except for those with $\Delta S = 0$, for which we could not even match $J_4$. The $Z_{\text{ext}}$ value is the atmospheric Z, “Juno” corresponds to the observations of *Juno* with the $3\sigma$ error bars, “Juno + Kaspi2018” are the observations corrected by the winds estimated in Kaspi et al. (2018), and “Juno + Kaspi2017” includes the correction to the observed gravitational moments due to a differential rotation slower than 3000 km, from Kaspi et al. (2017). While several models are compatible with the *Juno* observations when the corrections due to differential rotation estimated by Kaspi et al. (2017) are taken into account, this becomes much more difficult when the correction derived from the odd gravitational moments by Kaspi et al. (2018) are considered. Note that in this latter case, none of the models with $\Delta S = 0$ or $\nabla Z < 0$ can reproduce the gravitational moments of Jupiter.

## 5.3. Optimized Jupiter Models

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_6$ 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_6$, 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_4$ to $J_{10}$ values for the correct $J_2$. In contrast, the temperature departs from the isentropic profile for $R \lesssim 0.9 \times R_J$, i.e., within most of the interior, by a difference $\Delta T \simeq +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_\odot$ 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_6$ (consistent with Kaspi et al. 2018). For this latter, the inner isentropic envelope occupies a very limited fraction of the planet, $0.89 \times R_J \lesssim R \lesssim 0.92 \times R_J$, and the diluted core
<figure><img src="image_12.png" /><figcaption>**Figure 9.** (a) Pressure, temperature, density, and specific entropy as a function of radius for an isentropic (dashed gray) structure of Jupiter, a model with a high *J*<sub>6</sub> value (blue), compatible with Kaspi et al. (2017), and a model with a lowered *J*<sub>6</sub> compatible with Kaspi et al. (2018) (red). (b) Hydrogen (X), helium (Y), and heavy element (Z) mass abundances for the same models. The black crosses correspond to *P* = 1 Mbar, about the region of hydrogen pressure dissociation/ionization. The outer and inner convective zones correspond to the regions of constant (homogeneous) composition and entropy, and the gradients are representative of the inhomogeneous region of the envelope and the outer part of the diluted core, respectively.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 10.** (a) Normalized mass as a function of radius for a model that matches the observational constraints and an isentropic model. (b) Mass of heavy elements (in Earth masses) as a function of radius for the same models.</figcaption></figure>

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*<sub>2</sub>. The mean heavy element mass fraction is

displayed in Figure 9(b). As discussed earlier, *Z* decreases between the outer and inner envelopes ($\nabla Z > 0$) because models with a continuously increasing *Z* ($\nabla Z < 0$) in the envelope, although not strictly excluded, require a very large entropy jump ($\Delta 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.

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_⊕$.

## 5.4. Supersolar Atmospheric Abundance of Heavy Elements

As discussed in Section 3, we have taken a very conservative lower limit for the true average metal content in the atmosphere, $Z_{\text{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_{\text{Gal}}^{\text{true}} \approx (2–3) \times Z_{\odot} \approx 0.04–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_{\text{ext}}$ for a given structure increases $|J_6|$ to $|J_{10}|$ and thus implies either a very strong differential rotation or a very large $\Delta S$ to preserve the moments. If $Z_{\text{ext}} \simeq 0.05$, models with a constant inward increase of $Z$ in the metallization boundary region lead to a $\Delta S$ much larger than the $1 k_B$/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_{\text{ext}} > Z_{\odot}$ are shown in Figure 8 (orange circles). We see that for a given $\Delta S$, these models have higher $|J_6|$ to $|J_{10}|$ values than models with $Z_{\text{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_2$ and $J_4$ values that are consistent with Kaspi et al. (2018). Because of the strong correlation between $J_4$ and $J_6$, 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_{\text{ext}}$ will be examined in a subsequent paper.

# 6. Discussion

In this section, we examine the reliability of the various assumptions we used in the models in detail.

## 6.1. Hydrogen Pressure Metallization and H/He Phase Separation

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(\text{H-H}_2)$, and H/He demixion, $T_c(\text{H-He})$ (for $x_{\text{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(\text{H-H}_2) < T(P) < T_c(\text{H-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^+–\text{He}$ phase separation, as explored in Section 4.1.1, and/or external accretion events, depending on the exact value of $Z_{\text{ext}}$.

2. Sector II: $T_c(\text{H-H}_2) \sim T_c(\text{H-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 droplets (Section 4.1.2).

3. Sector III: if $T_c(\text{H-H}_2) > T_c(\text{H-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_{\text{ext}} \ge Z_{\text{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_6$) or Kaspi et al. (2017) (high $J_6$) correction due to differential rotation. Figure 11(a) corresponds to interior structures with strongly superadiabatic layered convection occurring at $P \ge 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
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Temperature-pressure profiles of different Jupiter models, overplotted on the immiscibility diagrams of Lorenzen et al. (2009; pale yellow area) and Morales et al. (2013b; orange area). In all models (except for the isentrope), $J_2$ and $J_4$ match the *Juno* values corrected by the estimated differential rotation of Kaspi et al. (2017), and a mass abundance of heavy elements in the inner envelope $Z_2 = 0.005$. The $Y_2$ value is the mass abundance of helium in the inner envelope, as in Figure 5. “Isentrope”: fully isentropic structure, with $Y = 0.275$. “High $J_6$”: models with $J_6$ compatible with Kaspi et al. (2017) but not Kaspi et al. (2018). “Low $J_6$”: model with $J_6$ compatible with Kaspi et al. (2018). This latter model seems difficult to justify physically (see text). (a) Models with layered convection of $P \ge 0.1$ Mbar. (b) The “Morales” and “Lorenzen” profiles are isentropic up to the onset of immiscibility according to the two related phase diagrams (see text).</figcaption></figure>

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 immiscibility for our model with $Y_2 = 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).

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_6$ 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 $\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).

## 6.2. Layered Convection

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 $\Delta S$ in the hydrogen metallization region requires an inward decrease of the heavy element abundance in this region, i.e., a locally positive gradient, $\nabla 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, $\nabla_\mu = (d \log \mu/d \log P) > 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)

$$0 < \nabla_T - \nabla_{ad} < \frac{\alpha_\mu}{\alpha_T} \nabla_\mu, \qquad (9)$$

where $(\alpha_\mu = (\partial \ln \rho/\partial \ln \mu)_{P,T}$ and $\alpha_T = (\partial \ln \rho/\partial \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).
<figure><img src="image_15.png" /><figcaption>**Figure 12.** (a) $\alpha_\mu/\alpha_T$ as a function of pressure for an isentropic model and two models with $\Delta S = 0.5$ and $1\,k_B/\text{proton}$, respectively. (b) $\alpha$ (Equation (11)) and $R_\rho^{-1}$ (see Leconte & Chabrier 2012) for the model described in the text.</figcaption></figure>

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 $\alpha$, which is the ratio of the size of the convective layer to the pressure scale height, $\alpha = l/H_p$. From their Equation (21), we can relate this parameter to the superadiabatic gradient, $(\langle\nabla_T\rangle - \nabla_{ad})$, by

$$ \langle\nabla_T\rangle - \nabla_{ad} = \epsilon_d \times \left[ (\Phi_0 \alpha^4 \epsilon_d)^{\frac{-1}{4(1+\alpha)}} + (\Phi_0 \alpha^4 \epsilon_d)^{\frac{-\alpha}{(1+\alpha)}} \right], \quad (10) $$

with layered convection occurring when

$$ 10^{-9} - 10^{-6} \lesssim \alpha \lesssim 10^{-4} - 10^{-2}, \quad (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_B/\text{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 $H_2$ 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 $\alpha$ and of the parameter $R_\rho^{-1} = (\alpha_\mu/\alpha_T)\nabla_\mu/(\nabla_T - \nabla_{ad})$ (overstable convection occurs for $R_\rho^{-1} > 1$, see Rosenblum et al. 2011; Leconte & Chabrier 2012; Mirouh et al. 2012), for a composition change from ($Y_{\text{ext}} = 0.23$, $Z_{\text{ext}} = 0.03$) to ($Y_2 = 0.3$, $Z_2 = 0.01$), with an entropy increase $\Delta S = 0.45\,k_B/\text{proton}$ between 0.1 and 1 Mbar. We see that $\alpha$ and $R_\rho^{-1}$ fulfill the conditions for the presence of layered convection in this

domain. Models with higher $\Delta S$ (of ~0.6 $k_B/\text{proton}$ at most, see Section 5.1) require larger $\Delta Y$. Globally, we verified that all our favored models do fulfill the conditions for the occurrence of layered convection derived in Leconte & Chabrier (2012).

One word of caution should be noted: when $H_2$ dissociates into atomic $H^+$, the mean molecular weight $\mu$ 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 $\mu$ due to $H_2$ 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.

### 6.3. External Impacts and Atmospheric Dynamical Effects

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 constraints 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

$$ v_{\text{conv}} \approx 10 \left( \frac{\nabla - \nabla_{\text{ad}}}{10^{-8}} \right)^{1/2} \left( \frac{H_p}{10^6 \text{ m}} \right) \text{m s}^{-1}, \quad (12) $$

with $(\nabla - \nabla_{ad})$ ranging from ~$10^{-10}$ at the center to ~$10^{-6}$ near the surface, yielding $v_{\text{conv}} \sim$ a few to about ~100 cm s$^{-1}$ 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_{\text{conv}} \approx 1 \text{ m 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 ammonia profile. However, when we compare Figure 4 of these authors with the location at which the *Galileo* entry probe dived (around $5^\circ-10^\circ$ in planetocentric latitude), we see that compared to the region deeper than 100 bar, *Galileo* can only have measured a lower limit of the ammonia content. Indeed, there seems to be a maximum ammonia 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 \times 10^{-6}$ and $2 \times 10^{-7}$, 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.

## 6.4. Magnetic Field

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*<sub>rms</sub>*d*/*η*, where *v*<sub>rms</sub>, *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 *ρ* ≈ 800 kg m⁻³, the radius *R* ≈ 0.85 × *R*<sub>J</sub>, 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*<sub>J</sub> ≈ 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*<sub>J</sub>. 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 *Rm* ≈ 50, and thus for the reservoir of convective power to still contribute appreciably to the dynamo action.

Defining *R*<sub>50</sub> as the radius in the planet above which *Rm* ≲ 50, Duarte et al. (2018) find that *R*<sub>50</sub> ≈ 0.9 *R*<sub>J</sub>, due essentially to the great change in conductivity when molecular H₂ fully recombines, while values below *R*<sub>50</sub> ≈ 0.85 *R*<sub>J</sub> 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 responsible 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.
## 6.5. Evolution

Constraints due to Jupiter's evolution have been briefly examined in Section 4.2. Our favored models suggest an entropy jump $\Delta S \approx 0.3-1$ kB/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 $\Delta 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 droplet sedimentation, and (iii) core erosion. Determining the resulting impact of these three processes upon the global cooling history of the planet is a highly non-trivial task that can hardly be intuited or inferred with simplistic models.

## 6.6. Does the Observed Outer Condition Lie on an Adiabat?

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 \simeq 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 \simeq 10$ bar in the external envelope of Jupiter, the optical depth is $\tau \gg 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 $\rho \sim 1$ kg m$^{-3}$, $g = GM_J/R^2 \simeq 20$ m s$^{-2}$, $T \sim 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$^{-2}$, is transported by convection and using the usual equations of the mixing length theory (Kippenhahn & Weigert 1990), this yields a typical super-adiabatic gradient $\Delta\nabla T_0 = (\nabla - \nabla_{ad})_0 \approx 10^{-6}$ at this pressure level, i.e., a characteristic convective velocity $v_{conv} \lesssim 1$ m s$^{-1}$. Because the angular velocity of Jupiter is $\Omega_J = v_{rot}/R_J = 1.76 \times 10^{-4}$ rad s$^{-1}$, the ratio of inertial to Coriolis forces, known as the Rossby number, at 10 bar is thus $R_0 = v_{conv}/(\Omega_J) \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\Omega^2 \ll g$, which allows performing a linear stability analysis of the MLT equations in the presence of a Coriolis force, $2\Omega \times v_{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 suradiabatic gradient in the presence of rotation in the region probed by Galileo (Stevenson 1979)

$$ \Delta\nabla T_{\sim} (\Delta\nabla T_0)^{3/5} \left(\frac{\Omega^2 l}{g}\right)^{2/5} \simeq 6 \times 10^{-6}. \quad (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.

## 7. Conclusion

In this paper, we have examined models of Jupiter aimed at matching both the most recent Juno observations and the Galileo constraints. Our calculations were carried out with the concentric MacLaurin spheroid method with all the necessary
mathematical and numerical constraints (Debras & Chabrier 2018). Because of the tension that is due to the high observed abundances of helium and heavy elements in the external envelope and the low values of the high-order gravitational moments, the number of possible interior models is very limited. We verified that even if the 1 bar temperature observed by *Galileo* departs from an adiabat because of the impact of rotation, the departure remains modest enough to take the *P* = 1 bar, *T* = 165 K observed values as the external isentrope conditions.

We first showed that the new data from *Juno* cannot be reproduced with conventional two- or three-isentropic homogeneous layer models. These latter are not able to match both the values of the gravitational moments and the external abundance of metals, which confirms the analysis of Wahl et al. (2017). The first conclusion is that there must be at least two regions of compositional gradient within the planet interior.

Our thorough analysis suggests that the planet should be composed of at least four main regions, moving inward from the surface: (I) the external isentropic, homogeneous molecular/atomic H₂/He/Z envelope, extending downward to about 93% of the planet radius, (II) an inhomogeneous domain around *P* ~ 0.1–2 Mbar, encompassing the region of hydrogen pressure ionization, with a size of about ~10%–15% of the radius, characterized by a gradient of composition (∇X, ∇Y, ∇Z), and an inward positive entropy change, Δ*S* > 0 (i.e., a locally negative entropy gradient, ∇*S* = (*∂S*/∂*r*) < 0), (III) a second inner isentropic, homogeneous, metallic hydrogen envelope, extending from the bottom of region II down to 60%–70% of the radius, lying on a hotter isentrope than the outer envelope (*S*<sub>III</sub> > *S*<sub>I</sub>) with, most likely, a smaller metal mass fraction than in the outer homogeneous envelope (*Z*<sub>III</sub> < *Z*<sub>I</sub>), and (IV) a diluted Z-rich core composed of volatiles that exhibits a compositional gradient. Potentially, a central compact seed can be present, which is essentially composed of solid iron and silicates.

A main result of this study is that a substantial entropy increase, ΔS ≳ 0.3 k<sub>B</sub>/proton, is necessary in the inhomogeneous region II to fulfill all the observational constraints. If this is not the case, one needs to invoke very strong differential rotation to explain the values of the high-order gravitational moments, at odds with the estimate of Kaspi et al. (2018; and even Kaspi et al. 2017 if ΔS = 0). This suggests the occurrence of either superadiabatic layered convection and/or a first-order phase transition, be it hydrogen pressure ionization or H/He phase separation. If this entropy increase lies in the range 0.3 ≲ ΔS ≲ 1 k<sub>B</sub>/proton, which seems to be inferred from various relevant physical processes, the abundance Z of heavy elements in region II must be locally decreasing, i.e., exhibiting a positive gradient of composition, ∇Z > 0, but an increasing molecular weight, i.e., a negative molecular weight gradient, ∇μ < 0, because the helium fraction at the bottom of region II is much greater. In case of a strong entropy increase in region II, ΔS > 0, it is possible to strongly reduce the value of the high-order gravitational moments while still fulfilling the external metal abundance reported by *Galileo* by invoking the presence of a central compact core. The first impact of this latter is to restrain the mass domain of Jupiter's interior impacting the moments. Although it is possible to find models with an inward-increasing metal abundance within region II, ∇Z < 0, compatible with *Juno* and *Galileo*, they require such a

large ΔS or amount of differential rotation that it seems hardly possible to justify them on physical grounds.

Note that there is a degeneracy of solutions between the change of entropy ΔS in region II and the outer differential rotation. The stronger ΔS, the shallower and weaker the differential rotation, which eventually enables values consistent with the estimate of Kaspi et al. (2018). In contrast, if differential rotation extends deeper inward and/or is stronger than suggested by these authors, the change in entropy across the boundary region can be significantly lowered. According to the study of Cao & Stevenson (2017), however, the differential rotation cannot extend too deep, as magnetic reconnection eventually occurs deep in the envelope, leading to rigid rotation.

The entropy jump ΔS in region II is also related to and can be constrained by other conditions, namely, (i) the mass of the central dense core M<sub>c</sub>. Indeed, as shown in the study, the mass of the central core is directly correlated with ΔS (the larger ΔS, the larger M<sub>c</sub>) and then anticorrelated with the amplitude of the high-order gravitational moments. (ii) The gradient of helium and heavy elements within the boundary region II: the larger the increase in *Y* and *Z* between regions I and III (most probably an increase in *Y* and a decrease in *Z*), the larger the ΔS required to reproduce the *Juno* data. (iii) At last, ΔS is constrained by the physics of dense matter, namely the nature of hydrogen pressure ionization (critical temperature and pressure, and entropy discontinuity) and by the miscibility diagram not only of H/He, but also of the various dominant heavy elements in metallic hydrogen. Finally, it is worth pointing out that even in the absence of a first-order transition, region II, characterized by a strong compositional change, is prone to layered convection. As examined in Section 6, the inferred profile is indeed consistent with conditions derived in Leconte & Chabrier (2012) for the presence of layered convection.

The inward decrease in mean heavy element mass fraction in region II, and thus the oversolar value in the upper envelope inferred from *Galileo*, can have different explanations. When the local temperature at the H<sub>2</sub>-H<sup>+</sup> metallization pressure P<sub>c</sub> is lower than the critical (PPT) temperature (*T* < *T*<sub>c</sub>), nucleation of H<sup>+</sup>-rich bubbles can occur, which is associated with some heavy elements, and move upward across the critical line. This continuously enriches the upper envelope I in (some) heavy elements. Because hydrogen ionization immediately triggers H/He phase separation, with the formation of drowning He-rich droplets, this process at the same time yields an enrichment of He and associated species in the lower envelope III. When the above temperature condition is not fulfilled, hydrogen pressure ionization occurs smoothly, there is no bubble nucleation. In that case, in order to enrich the upper envelope I in heavy elements, there must be either an immiscibility of some species in the H/He/Z mixture at the relevant temperature and pressure, yielding a high equilibrium concentration of these species in the low-pressure, low-temperature molecular phase, or persistent layered convection. If the real enrichment in heavy elements is largely oversolar (i.e., much higher than the *Galileo* value), the occurrence of external impacts during Jupiter's history seems to be inevitable to explain it. In all cases, it seems difficult to avoid the presence of a first-order transition or persistent superadiabatic layered convection in Jupiter's gaseous envelope around the ~Mbar region. Accurate calculations of the long-term evolution of the
planet are definitely needed to assess or reject the viability of some of these static models. As mentioned in Section 6.5, however, properly handling such calculations appears to be a task of enormous complexity.

In conclusion, we have derived in this study interior models of Jupiter that are able to reproduce all the observed stringent gravitational constraints from the *Juno* mission and the outer helium and heavy element abundances from *Galileo*. These models differ appreciably from all Jupiter models that have been derived previously, which ignored either the *Juno* or the *Galileo* constraints, making these models (and related papers or reviews) obsolete. As shown above, however, because of the lack of a precise characterization of the main physical processes, there is still a degeneracy of possible models. Indeed, neither experimental nor numerical explorations of these processes have provided definitive information about the related fundamental questions so far. This illustrates the tight link between fundamental physics and astrophysics. Additional constraints also arise from the differential rotation in the planet. Indeed, high-order gravitational moments are essentially only sensitive to the outermost region of the planet, constraining the available level of differential rotation (see notably Hubbard 1999). As explained above, more constraints on differential rotation will help to constrain the change in entropy in the pressure ionization boundary domain, and subsequently, the mass or even the presence of the central compact core. This issue will be explored in a forthcoming paper. Finally, it is worth mentioning that the most favored models that can fulfill the constraints of both *Galileo* and *Juno* according to this study are basically those, or are at least among those, that have been intuited and explored in great detail by Stevenson & Salpeter 40 years ago in their two seminal papers and by Stevenson (1985)!

We wish to thank T. Fouchet, E. Jaupart, G. Laibe, and D. Stevenson for helpful discussions. G.C. acknowledges the warm hospitality of the OWL Institute at the University of Santa Cruz, where part of this work was accomplished. This work was supported by the Programme National de Planéto- logie (PNP) of CNRS-INSU cofunded by CNES.

## ORCID iDs

Florian Debras https://orcid.org/0000-0002-3113-4840

Gilles Chabrier https://orcid.org/0000-0002-8342-9149

## References

* Anders, E., & Grevesse, N. 1989, GeCoA, 53, 197

* Archinal, B. A., A'Hearn, M. F., Bowell, E., et al. 2011, CeMDA, 109, 101

* Ashcroft, N. W. 2004, PhRvL, 92, 187002

* Berardo, D., & Cumming, A. 2017, ApJL, 846, L17

* Bézard, B., Lellouch, E., Strobel, D., Maillard, J.-P., & Drossart, P. 2002, Icar, 159, 95

* Bjoraker, G. L., Wong, M. H., de Pater, I., et al. 2018, AJ, 156, 101

* Blackburn, R. 1966, Mettalurgical Review, 11, 159

* Bolton, S. J., Adriani, A., Adumitroaie, V., et al. 2017, Sci, 356, 821

* Campbell, J. K., & Synott, S. P. 1985, AJ, 90, 364

* Cao, H., & Stevenson, D. J. 2017, Icar, 296, 59

* Celliers, P. M., Millot, M., Brygoo, S., et al. 2018, Sci, 361, 677

* Chabrier, G., & Baraffe, I. 2007, ApJL, 661, L81

* Chabrier, G., Gallardo, J., & Baraffe, I. 2007, A&A, 472, L17

* Chabrier, G., Mazevet, S., & Soubiran, F. 2019, ApJ, in press

* Chabrier, G., Saumon, D., Hubbard, W. B., & Lunine, J. I. 1992, ApJ, 391, 817

* Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon)

* Christensen, U. R., & Aubert, J. 2006, GeoJI, 166, 97

* Christensen, U. R., & Yuen, D. A. 1985, JGR, 90, 10291

* Cohen, E. R., & Taylor, B. N. 1987, Rev. Mod. Phys., 59, 1121

* Connerney, J. E. P., Kotsiaros, S., Oliversen, R. J., et al. 2018, GeoRL, 45, 2590

* Debras, F., & Chabrier, G. 2018, A&A, 609, A97

* Duarte, L. D. V., Wicht, J., & Gastine, T. 2018, Icar, 299, 206

* Fisher, B. M., Orton, G. S., Liu, J., et al. 2016, Icar, 280, 268

* Folkner, W. M., Iess, L., Anderson, J. D., et al. 2017, GeoRL, 44, 4694

* Fortney, J. J., & Hubbard, W. B. 2003, Icar, 164, 228

* Gastine, T., Wicht, J., Duarte, L. D. V., Heimpel, M., & Becker, A. 2014, GeoRL, 41, 5410

* González-Cataldo, F., Wilson, H. F., & Militzer, B. 2014, ApJ, 787, 79

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

* Guillot, T., Stevenson, D. J., Hubbard, W. B., & Saumon, D. 2004, in Jupiter. The Planet, Satellites and Magnetosphere, ed. F. Bagenal, T. E. Dowling, & W. B. McKinnon (Cambridge: Cambridge Univ. Press), 35

* Hansen, C. J., & Kawaler, S. D. 1994, Stellar Interiors, Physical Principles, Structure, and Evolution (Berlin: Springer)

* Helled, R., & Stevenson, D. 2017, ApJL, 840, L4

* Hubbard, W. B. 1999, Icar, 137, 357

* Hubbard, W. B. 2012, ApJL, 756, L15

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

* Hyman, H. H. 1964, JChEd, 41, 174

* Iaroslavitz, E., & Podolak, M. 2007, Icar, 187, 600

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

* Kaspi, Y., Galanti, E., Hubbard, W. B., & Iess, L. 2018, Natur, 555, 223

* Kaspi, Y., Guillot, T., Galanti, E., et al. 2017, GeoRL, 44, 5960

* Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution (Berlin: Springer)

* Knudson, M. D., Desjarlais, M. P., Becker, A., et al. 2015, Sci, 348, 1455

* Knudson, M. D., Desjarlais, M. P., Lemke, R. W., et al. 2012, PhRvL, 108, 091102

* Kurokawa, H., & Inutsuka, S.-i. 2015, ApJ, 815, 78

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

* Leconte, J., & Chabrier, G. 2013, NatGe, 6, 347

* Li, C., Ingersoll, A., Janssen, M., et al. 2017, GeoRL, 44, 5317

* Licari, A. 2016, PhD thesis, ENS de Lyon

* Lorenzen, W., Holst, B., & Redmer, R. 2009, PhRvL, 102, 115701

* Lorenzen, W., Holst, B., & Redmer, R. 2011, PhRvB, 84, 235109

* Loubeyre, P., Le Toullec, R., & Pinceaux, J.-P. 1985, PhRvB, 32, 7611

* Lyon, S. P., & Johnson, J. D. 1992, ANL Rep. LA-UR-92-3407 (Los Alamos, NM: LANL)

* Mazevet, S., Licari, A., Chabrier, G., & Potekhin, A. Y. 2019, A&A, 621, A128

* Mazzola, G., Helled, R., & Sorella, S. 2018, PhRvL, 120, 025701

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

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

* Mirouh, G. M., Garaud, P., Stellmach, S., Traxler, A. L., & Wood, T. S. 2012, ApJ, 750, 61

* Moll, R., Garaud, P., Mankovich, C., & Fortney, J. J. 2017, ApJ, 849, 24

* Moll, R., Garaud, P., & Stellmach, S. 2016, ApJ, 823, 33

* Moore, K. M., Yadav, R. K., Kulowski, L., et al. 2018, Natur, 561, 76

* Morales, M. A., Hamel, S., Caspersen, K., & Schwegler, E. 2013a, PhRvB, 87, 174105

* Morales, M. A., McMahon, J. M., Pierleoni, C., & Ceperley, D. M. 2013b, PhRvL, 110, 065702

* Morales, M. A., Pierleoni, C., & Ceperley, D. M. 2010, PhRvE, 81, 021202

* Morales, M. A., Schwegler, E., Ceperley, D., et al. 2009, PNAS, 106, 1324

* Moreno, R., Marten, A., Matthews, H. E., & Biraud, Y. 2003, P&SS, 51, 591

* Musella, R., Mazevet, S., & Guyot, F. 2018, arXiv:1805.12439

* Nellis, W. J., Ross, M., & Holmes, N. C. 1995, Sci, 269, 1249

* Nettelmann, N., Becker, A., Holst, B., & Redmer, R. 2012, ApJ, 750, 52

* Pépin, C. M., Geneste, G., Dewaele, A., Mezouar, M., & Loubeyre, P. 2017, Sci, 357, 382

* Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62

* Ramsay, W., & Travers, R. W. 1897, RSPS, 61, 267

* Riddle, A. C., & Warwick, J. W. 1976, Icar, 27, 457

* Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S. 2011, ApJ, 731, 66

* Salpeter, E. E. 1973, ApJL, 181, L83

* Saumon, D., & Chabrier, G. 1992, PhRvA, 46, 2084

* Saumon, D., & Guillot, T. 2004, ApJ, 609, 1170

* Schneider, T., & Liu, J. 2009, JAtS, 66, 579

* Schöttler, M., & Redmer, G. 2018, PhRvL, 120, 115703

* Soubiran, F., Mazevet, S., Winisdoerffer, C., & Chabrier, G. 2013, PhRvB, 87, 165114
Soubiran, F., & Militzer, B. 1985, Icar, 62, 4

Soubiran, F., & Militzer, B. 2015, ApJ, 806, 228

Soubiran, F., & Militzer, B. 2016, ApJ, 829, 14

Stern, M. E. 1960, TellA, 12, 172

Stevenson, D. J. 1979, GApFD, 12, 139

Stevenson, D. J. 1985, Icar, 62, 4

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

Stevenson, D. J., & Salpeter, E. E. 1977b, ApJS, 35, 221

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

von Zahn, U., Hunten, D. M., & Lehmacher, G. 1998, JGR, 103, 22815

Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, GeoRL, 44, 4649

Wahl, S. M., Wilson, H. F., & Militzer, B. 2013, ApJ, 773, 95

Wilson, H. F., & Militzer, B. 2010, PhRvL, 104, 121101

Wilson, H. F., & Militzer, B. 2012a, PhRvL, 108, 111101

Wilson, H. F., & Militzer, B. 2012b, ApJ, 745, 54

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

Wood, T. S., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 157