# Study of Jupiter’s Interior: Comparison of 2, 3, 4, 5, and 6 Layer Models

Burkhard Militzer  Department of Earth and Planetary Science, Department of
Astronomy,
University of California, Berkeley, CA, 94720, USA  William B. Hubbard  Lunar
and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721, USA

###### Abstract

With the goal of matching spacecraft measurements from Juno and Galileo
missions, we construct ensembles of 2, 3, 4, 5, and 6 layer models for
Jupiter’s interior. All except our two layer models can match the planet’s
gravity field as measured by the Juno spacecraft. We find, however, that some
model types are more plausible than others. In the best three layer models,
for example, the transition from molecular to metallic hydrogen needs to be at
$\sim$500 GPa while theory and experiments place this transition at $\sim$100
GPa. Four layer models with a single sharp boundary between core and mantle
would be short-lived due to rapid convective core erosion. For this reason, we
favor our five layer models that include a dilute core surrounded by a stably
stratified core transition layer. Six layer models with a small compact core
are also possible but with an upper limit of 3 Earth masses for such a compact
core. All models assume a 1 bar temperature of 166.1 K, employ physical
equations of state, and are constructed with the nonperturbative Concentric
Maclaurin Spheroid (CMS) method. We analyze the convergence of this method and
describe technical steps that are needed to make this technique so efficient
that ensembles of models can be generated.

Giant planets, Jupiter’s interior, gravity science

##  1 Introduction

The space missions Juno (Bolton et al., 2017) and Cassini (Spilker, 2019) have
provided us with a wealth of new data for Jupiter and Saturn, the two largest
planets in our solar system. Multiple close flybys have determined the gravity
fields of these planets with exquisite precision, far exceeding the earlier
measurements by the Pioneer and Voyager missions (Campbell & Synnott, 1985;
Campbell & Anderson, 1989). The earlier missions determined Jupiter gravity
harmonics to be $J_{2}=14697\pm 1$, $J_{4}=-584\pm 5$, and $J_{6}=31\pm 20$
(all multiplied by 106 and adjusted to an equatorial radius of 71492 km) while
today these coefficients are known to three orders of magnitude higher
precision: $J_{2}=14696.5735\pm 0.0017$, $J_{4}=‐586.6085\pm 0.0024$, and
$J_{6}=34.2007\pm 0.0067$ (Durante et al., 2020). This improvement has led to
revision in the models of Jupiter’s interior required to match these gravity
measurements. When interior models were only constrained by Jupiter’s mass,
equatorial radius, $J_{2}$ and $J_{4}$, traditional three layer models for
Jupiter’s interior sufficed (Stevenson, 1982; Guillot, 2005). These models
typically included a compact core comprising up to 100% of heavy elements,
presumably left over from the planet’s formation by core accretion
(Bodenheimer & Pollack, 1986). All these three layer models assumed Jupiter’s
envelope to be inhomogeneous (Saumon & Guillot, 2004; Militzer et al., 2008;
Nettelmann et al., 2012), introducing a discontinuity in the helium or in the
heavy elements fraction at some radius or pressure in the envelope. The
motivation for the discontinuity was in part empirical because one needed an
additional free parameter to match $J_{4}$ while adjusting the core mass and
heavy element fraction in the envelope enabled one to match the planet’s mass
and $J_{2}$. An exception were the homogeneous models reported by Militzer et
al. (2008) who invoked deep differential rotation to match $J_{4}$ instead.
However, the required wind depth of $\sim$10 000 km far exceeded values
implied by Jupiter’s odd gravity harmonics (Kaspi et al., 2018) and typical
magnetic field models (Christensen et al., 2020). Both these latter approaches
favor a wind depth of $\sim$3000 km although an agreement has yet to be
reached on precisely how the wind speeds decay with depth.

A variety of model assumptions and techniques have been invoked to model
Jupiter’s interior structure in the past. Ni (2018) combined the theory of
figures (Zharkov & Trubitsyn, 1978) with the approach by Anderson & Schubert
(2007) to adjust coefficients of an empirical polynomial function that
represents the Jupiter’s interior density profile in order to match the Juno
gravity measurements. Nettelmann et al. (2021) extended the ToF to technique
to seventh order. Militzer & Hubbard (2023) constructed physical and abstract
models of Jupiter’s interior structure to demonstrate that the gravity
coefficients $J_{2}$, $J_{4}$ and $J_{6}$ are sufficient to constrain
Jupiter’s moment of inertia very tightly. The physical models rely on
realistic equations of state and assumptions for the compositions while
abstract models represent the planet’s interior by a small number of constant-
density. Howard et al. (2023) studied the effects of various published
equation of state (EOS) tables and also introduced an EOS modification. Moll
et al. (2017), Müller et al. (2020), and Helled et al. (2022) studied the
Jupiter’s interior convection and the evolution of a primordial, compact core
that was originally composed to 100% of heavy elements. Liu et al. (2019)
investigated whether such a core could be diluted by a giant impact. Core
dilution is plausible because ab initio computer simulations have shown that
all typical core materials such as water, silicates and iron are soluble in
metallic hydrogen at megabar pressures (Wilson & Militzer, 2012a, b; Wahl et
al., 2013; Gonzalez-Cataldo et al., 2014).

Kerley (2004) first constructed a homogeneous baseline model to demonstrate
that the pre-Juno gravity data cannot be matched unless the envelope is
represented by at least two distinct layers. The need for Jupiter’s envelope
to be inhomogeneous was recently reconfirmed by Miguel et al. (2022). The
original justification for the inhomogeneity was a first-order phase
transition from molecular to metallic hydrogen. However, ab initio computer
simulations (Vorberger et al., 2010, 2012) have since shown that such a
transition likely exists only at low temperature in hydrogen (Morales et al.,
2010) while the temperatures in giant planet interiors are everywhere above
the critical point so that molecular dissociation there occurs gradually. So
this process by itself is insufficient to cause an inhomogeneity.
Alternatively, the phase separation of hydrogen-helium mixtures (Morales et
al., 2013) has been considered but so far there is insufficient evidence that
the falling helium droplets will sequester any heavy elements besides neon
(Wilson & Militzer, 2010).

In this article, we incorporate the nonperturbative concentric Maclaurin
spheroid (CMS) method (Hubbard, 2013) into a streamlined modeling procedure
which achieves the high precision needed to exploit Juno constraints, while
automatically matching basic constraints such as mass, equatorial radius, and
low-degree zonal gravity harmonics. With Monte Carlo sampling (Militzer,
2023a) we then efficiently generate large numbers of models with increasing
degrees of complexity. We construct models for Jupiter’s interior assuming
different compositions and layered structures. We invoke one consistent set of
physical assumptions. We adopt one equation of state derived with ab initio
simulations (Militzer & Hubbard, 2013) without introducing any corrections to
lower its density, nor do we vary the temperature at 1 bar, 166.1 K measured
by the Galileo entry probe (Niemann et al., 1998; Wong et al., 2004). This
temperature value is also consistent with earlier radio occultation
measurements by the Voyager spacecraft that yielded a value of 165 $\pm$ 5 K
(Lindal et al., 1981). These remote observations very recently re-analyzed by
Gupta et al. (2022) who determined higher temperatures of 167$\pm$4 and
170$\pm$4 K for latitudes of 6∘S and 12∘N respectively.

We construct models that have between two and six layers as we illustrate in
Fig. 1. We start from a reference model with five layers (Militzer et al.,
2022) that contain an outer layer of molecular hydrogen, a helium rain layer,
a layer of metallic hydrogen, a core transition layer, and a prominent dilute
core in the center that may extend up to 60% of the planet’s radius. We
generate two types of four layer models by either eliminating the helium rain
layer or the core transition layer from our five layer reference models. Six
layer models are derived by inserting a compact core into the reference
models. Under our model assumptions, we will demonstrate that such a compact
core cannot be heavier than 3 Earth masses for the models to be compatible
with Jupiter’s gravity measurements. We also generate alternate types of five
layer models by removing either the helium rain layer, metallic hydrogen
layer, or the core transition layer.

Finally we revisit the traditional three layer models and even try to
construct two layer models with just a core and an homogeneous envelope. These
two layer models do not fit the gravity field, while our three, four, five,
and six layer models all do. Still some of these models are more plausible
than others. For example, matching the gravity field with three layer models
requires us to place a transition between a helium-depleted outer layer and an
helium-enriched inner layer at a pressure of 500 GPa, which is incompatible
with the physics of hydrogen-helium mixtures (Morales et al., 2013).

This paper is organized as follows. In section 2, we describe our model
assumptions and explain how we efficiently converge on model solutions for the
interior structure and winds that match the planet’s radius, mass, gravity
harmonic $J_{2}$, and protosolar helium abundance. In section 3, we discuss
the results for the different model types. We show that two layer models are
unable to match $J_{4}$ and explain why three layer models need to have a
molecular-to-metallic transition at $\sim$500 GPa. We then analyze the
properties of 4, 5, and 6 layer models. We describe how we place constraints
on the mass of the compact core and the helium distribution within the helium
rain layer. We conclude with section 4.

Figure 1: Illustration of our 2, 3, 4, 5 and 6 layer models in Tab. 1. At
every layer boundary we specify the pressure, temperature, and fractional
volumetric radius. The collection of thin solid lines represent the core
transition layer that may surround the dilute core. In the center is our
reference model with five layers. The two types of four layer models are
generated by removing either the helium rain or the core transition layers.
The five layer models B, C, and D were obtained by removing one layer one from
the reference model while adding a compact core. Such a core was included in
our final six layer model. It also distinguishes our simpler two and three
layer models.

##  2 Methods and Model Assumptions

###  2.1 Methods for the Interior

In this paper, we employ the CMS method to construct models of rapidly
rotating giant planets that meet certain observational constraints such as
planet mass, equatorial radius, helium content, etc. In principle one has two
choices of how to meet such constraints because the CMS method essentially
provides a way to converge onto a hydrostatic solution for a given pressure-
density relationship by adjusting the shapes of all spheroids. One can then
wait until the CMS method has converged, evaluate how well the constraints
have been met, then vary certain model parameters, and restart the CMS
procedure. Alternatively, one can adjust certain model parameters during the
CMS convergence procedure so that certain criteria are met automatically as a
hydrostatic solution is found. The latter method is obviously much more
efficient because it typically does not increase the cost of the CMS procedure
and more importantly, it reduces the number of dimensions of the parameter
space that model optimization and Monte Carlo calculation need to explore. We
now discuss how we converge onto a valid planetary interior model that
satisfies a number of plausible constraints:

(1) For that reason, almost all CMS calculations fix the equatorial radius
$a$, assumed to have been measured with such high precision that its remaining
uncertainty has a negligible effect on all computed planet properties. Here we
use $a=71492$ km (n.b. the error bar on $a$ is $\pm 4$ km) and associate it
with a pressure of 1 bar. (Conversely, in Wahl et al. (2021), the equatorial
radius was varied during the CMS convergence in order to match the planet mass
as models of specific exoplanets with different rotation periods and
oblatenesses were constructed.)

(2a) Next we adjust the helium concentration in the interior so that the
planet overall matches the protosolar hydrogen-helium fraction of
$\tilde{Y}_{0}\equiv Y_{0}/(X_{0}+Y_{0})=0.2777$ derived from $X_{0}=0.7112$
and $Y_{0}=0.2735$ (Lodders, 2010). In the outer molecular layer, we set
$Y_{1}/(X_{1}+Y_{1})=0.2369$ to be consistent with measurements of the Galileo
entry probe, $Y/(X+Y)=0.238\pm 0.005$ (von Zahn et al., 1998). We assume the
helium mass fraction $\tilde{Y}_{1}\equiv Y_{1}/(X_{1}+Y_{1})$ remains
constant in the molecular layer up to a pressure $P_{\rm rain,1}$ and from
$P_{\rm rain,2}$ on, it is again constant at $\tilde{Y}_{2}\equiv
Y_{2}/(X_{2}+Y_{2})$. In the helium rain layer from $P_{\rm rain,1}$ to
$P_{\rm rain,2}$, we gradually switch from one value to the other as a
function of $\log(P)$,

$\displaystyle\tilde{Y}(P)$ $\displaystyle=$
$\displaystyle\tilde{Y}_{1}+F_{2}(x)\left[\tilde{Y}_{2}-\tilde{Y}_{1}\right]$
(1) $\displaystyle F_{2}(x)$ $\displaystyle=$ $\displaystyle x^{\alpha}$ (2)
$\displaystyle x$ $\displaystyle=$ $\displaystyle\frac{\log{P}-\log{P_{\rm
rain,1}}}{\log{P_{\rm rain,2}}-\log{P_{\rm rain,1}}}$ (3)

The switching function, $F_{2}$, represents the fraction of $\tilde{Y}_{2}$ in
the mixture and satisfies the expected values at the boundaries,
$F_{2}(P=P_{\rm rain,1})=0$ and $F_{2}(P=P_{\rm rain,2})=1$. Besides
representing a linear change from $\tilde{Y}_{1}$ to $\tilde{Y}_{2}$ when
$\alpha=1$, it introduces the option of having a large helium fraction from
that layer being sequestered by setting $\alpha\gg 1$ or increasing the helium
fraction in that layer by specifying $\alpha\ll 1$. The simple form of Eq. 2
has the advantage that the average fraction across the helium rain layer is
given by,

$\left<F_{2}\right>=\int_{0}^{1}dx\,F_{2}(x)=1/(1+\alpha)\;\;.$ (4)

(2b) Next we construct a spline function to represent the cumulative mass as a
function of $\log(P)$ so that we can derive the mass fraction that is stored
in the layer between two pressure values. The global helium fraction,
$\left<\right.\tilde{Y}\left.\right>$, is given by

$\displaystyle\left<\right.\tilde{Y}\left.\right>$ $\displaystyle=$
$\displaystyle m_{1}\tilde{Y}_{1}\;\;+\;\;m_{\rm rain}\tilde{Y}_{\rm
rain}\;\;+\;\;m_{2}\tilde{Y}_{2}\;,$ (5) $\displaystyle=$
$\displaystyle\left[m_{1}+m_{\rm
rain}(1-\left<F_{2}\right>)\right]\;\tilde{Y}_{1}\;\;+\;\;\left[m_{\rm
rain}\left<F_{2}\right>+m_{2}\right]\;\tilde{Y}_{2}\;,$ (6)

where $m_{1}$, $m_{\rm rain}$, and $m_{2}$ represent the fractional masses of
the three respective layers without the contributions from heavy elements. Eq.
6 enables us to adjust $\tilde{Y}_{2}$ so that the global helium abundance
matches the protosolar value,
$\left<\right.\tilde{Y}\left.\right>=\tilde{Y}_{0}$.

(3) Finally, we discuss how the CMS procedure can converge directly onto the
exact mass and $J_{2}$ value of the planet. We decided to match
$J_{2}=(14696.5063\pm 0.0006)\times 10^{-6}$ exactly because it has been
measured with exquisite precision. (As recommend by Durante et al. (2020), we
do not include the tidal effects on $J_{2}$, estimated to be
6.72$\times$10-8.) The planet mass has previously been matched by introducing
an overall density scaling factor $\beta$ (Militzer & Hubbard, 2013) or by
adjusting the density of a compact core represented by a single innermost
spheroid (Wahl et al., 2017). Since our current set of models do not have a
compact core or its mass is very small, it is no longer practical to match the
planet’s mass by adjusting the core mass. So we instead match the planet’s
mass and $J_{2}$ by simultaneously adjusting $Z_{1}$ and $Z_{2}$ as the CMS
method converges to a hydrostatic solution. $Z_{1}$ is the mass fraction of
heavy elements in the outer layers of the planet that we assume to be constant
from the atmosphere to the beginning of the core transition layer. Here our
assumptions differ from those by Debras & Chabrier (2019) who matched $J_{4}$
and $J_{6}$ by assuming Jupiter’s outer enveloped is not fully convective and
that the concentration of heavy elements decreased with increasing pressure.
In this work, $Z_{2}$ is the mass fraction of heavy elements in the dilute
core. In the core transition layer, we assume $Z$ to increase linearly as a
function of $\log(P)$ from $Z_{1}$ and $Z_{2}$.

We find that adjusting $Z_{1}$ and $Z_{2}$ to match $M$ and $J_{2}$ is much
more efficient than first adjusting $Z_{1}$ to match the mass in an inner loop
and then adjusting $Z_{2}$ to match $J_{2}$ in an outer loop. For example if a
bisection procedure or Brent’s method (Press et al., 2001) are employed as an
outer loop to bring the deviation $f(Z_{2})=J_{2}(Z_{2})-J_{2}^{\rm target}$
to zero, one has to converge multiple CMS models, which is inefficient. To
start, one has to evaluate $f(Z_{2})$ for two different values of $Z_{2}$ that
yield $f$ values of opposite signs before $Z_{2}$ can be refined iteratively.
This means that even if one already has a reasonably good starting value for
$Z_{2}$, one has to step away from it and construct another hydrostatic CMS
solution for a second $Z_{2}$ value that is ultimately not needed. This is
unnecessary and inefficient. To avoid such unnecessary steps, we define two
following functions,

$\displaystyle f_{1}(Z_{1},Z_{2})$ $\displaystyle\equiv$ $\displaystyle
M(Z_{1},Z_{2})\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\,-M_{\rm planet}$ (7)
$\displaystyle f_{2}(Z_{1},Z_{2})$ $\displaystyle\equiv$ $\displaystyle
M(Z_{1},Z_{2})\times J_{2}(Z_{1},Z_{2})-M_{\rm planet}\times J_{2}^{\rm
target}$ (8)

and employ the Newton-Raphson method to adjust $Z_{1}$ and $Z_{2}$ while the
CMS method converges to a hydrostatic solution that then matches $M$ and
$J_{2}$ automatically. This approach requires us to compute the derivatives of
$M$ and $(MJ_{2})$ with respect to $Z_{1}$ and $Z_{2}$, which we will explain
in the following equations. (Alternatively, one may define
$f_{2}(Z_{1},Z_{2})\equiv J_{2}(Z_{1},Z_{2})-J_{2}^{\rm target}$, which leads
to iteration equations of similar complexity.)

In rare cases, the Jacobian derivative matrix (see. Eq. A1) for a particular
Newton-Raphson step cannot be inverted or the Newton-Raphson method starts
oscillating between two states. If this happens we insert a 2D regula falsi
step to depart from such a pathological situation. Details of our regula falsi
method are given in appendix A. We find that this combined approach works very
well except in cases when $P_{\rm rain,1}$ and $P_{\rm rain,2}$ exceed the
central pressure of the planet and the variable $Z_{2}$ becomes irrelevant so
that we no longer have two independent variables to match $M$ and $J_{2}$
simultaneously.

Within the CMS method (Hubbard, 2013), the planet’s interior is represented by
a series of nested spheroids indexed by $j$, each with volume $V_{j}$. The
total mass is derived from,

$\displaystyle M$ $\displaystyle=$
$\displaystyle\sum_{j}\rho_{j}(V_{j}-V_{j+1})=\sum_{j}\delta_{j}V_{j}\;,$ (9)
$\displaystyle V_{j}$ $\displaystyle=$
$\displaystyle\frac{2\pi}{3}\lambda_{j}^{3}\int_{-1}^{+1}d\mu\;\zeta_{j}^{3}(\mu)\;,$
(10)

where $\delta_{j}=\rho_{j}-\rho_{j-1}$ represents the density difference of
spheroid $j$ and next outer spheroid $j-1$. Following Hubbard (2013), the
shape function, $\zeta_{j}(\mu)\equiv\zeta_{j,k}$, is represented by a series
of discrete quadrature points, $\mu_{k}=\cos(\theta_{k})$, where $\theta_{k}$
are colatitude values. The distance of the spheroid surface from the origin is
given by $r(\mu)=\lambda_{j}\zeta_{j}(\mu)$ where $\lambda_{j}$ is the
equatorial radius of spheroid $j$.

In a similar fashion, one can derive an expression for the product $(M\times
J_{n})$,

$\displaystyle(MJ_{n})$ $\displaystyle=$
$\displaystyle\sum_{j}\delta_{j}\hat{J}_{j,n}=\sum_{j}\rho_{j}\left(\hat{J}_{j,n}-\hat{J}_{j+1,n}\right)$
(11) $\displaystyle\hat{J}_{j,n}$ $\displaystyle=$
$\displaystyle\frac{-2\pi}{3+n}\;\lambda_{j}^{n+3}\int_{-1}^{+1}d\mu\;P_{n}(\mu)\;\zeta_{j}^{n+3}(\mu)\;.$
(12)

The coefficients $\tilde{J}_{j,n}$ from Hubbard (2013) are recovered through,
$\tilde{J}_{j,n}=\delta_{j}\hat{J}_{j,n}/(\lambda_{j}^{n}M)$. The conventional
gravity harmonics are obtained by

$J_{n}=\frac{1}{M}\sum_{j}\delta_{j}\hat{J}_{j,n}=\sum_{j}\lambda_{j}^{n}\tilde{J}_{j,n}\;.$
(13)

Equations 9 and 11 would enable us to determine the derivatives of $f_{1}$ and
$f_{2}$ with respect to $Z_{1}$ and $Z_{2}$ if we knew how the density,
$\rho(\tilde{Y},Z)$, varied with the contents of helium and heavy elements.
Following Hubbard & Militzer (2016), we invoke the additive volume rule to
incorporate the heavy elements into our equation of state (Militzer & Hubbard,
2013; Militzer, 2013), $\rho^{*}$, that was computed for a single helium
fraction, $Y^{*}$, without any heavy elements ($Y^{*}\equiv\tilde{Y}^{*}$),

$\frac{\rho^{*}}{\rho(\tilde{Y},Z)}=(1-Z)A+Z\frac{\rho^{*}}{\rho_{Z}}\;\;\;\;\rm{with}\;\;\;\;A=\frac{(1-\tilde{Y})+(\tilde{Y}-Y^{*})\frac{\rho^{*}}{\rho_{\rm
He}}}{1-Y^{*}}\quad.$ (14)

$\rho_{\rm He}$ is an EOS for pure helium (Militzer, 2006, 2009; Militzer et
al., 2021) that enables us to perturb the helium fraction from our standard
value, $Y^{*}$. We incorporate heavy elements in the envelope by following
Hubbard & Militzer (2016) who conducted ab initio computer simulations of a
solar-proportion mixture of H2O:CH4:NH3 and then represented the results by a
function $\rho_{Z}/\rho^{*}$ with that varies slowly pressure, which was found
to be in agreement with later work by Soubiran & Militzer (2016).

Equation 14 not only allows us to compute, $\rho(\tilde{Y},Z)$, it also makes
the evaluation of $\frac{d\rho(\tilde{Y},Z)}{dZ}$ straightforward. It allows
us to determine how the density of every layer $\rho_{j}$ in Eqs. 9 and 11
depends on $Z$ and therefore gives us access to the derivatives of $f_{1}$ and
$f_{2}$ that we need when we employ the Newton-Raphson iterations to match the
planet’s mass and $J_{2}$ by adjusting $Z_{1}$ and $Z_{2}$ as the CMS
procedure converges to a hydrostatic solution.

Figure 2:  Our preferred technique to efficiently converge onto a
selfconsistent solution that includes a hydrostatic description of the
interior with the CMS method and a compatible wind model. It is novel that the
planet’s mass, $J_{2}$, and global helium abundance are matched automatically
as the calculation converges without a large computational overhead.
Calculations of spheroid shapes and the wind contributions to the $J_{n}$ are
performed only once in the main loop because they are rather expensive. In the
inner loop, the gravity coefficients, densities, pressures, the helium
fraction $Y_{2}$ and heavy element abundances $Z_{1}$ and $Z_{2}$ are updated
more often because they scale linearly with spheroid number. If the model
includes a compact core of a specific mass target, the $\lambda$ grid is
updated in the inner loop to match that. Besides converging the gravity
coefficients, $J_{n}$, other criteria for mass, $J_{2}$, or the overall helium
abundance may be employed to terminate the main loop.

In Fig. 2, we describe our preferred way to efficiently converge a CMS
calculation. We start by setting up an optimized $\lambda$ grid (Militzer et
al., 2019) that anchors the equatorial radii of all spheroids. All other
spheroid points are adjusted so that an equipotential emerges. Based on an
earlier CMS calculation for Jupiter, the $\lambda$ grid is chosen so that the
corresponding density values fall on a logarithmic grid such that
$\rho(\lambda_{j})/\rho(\lambda_{j+1})=$constant. This was shown to be a good
choice (Militzer et al., 2019) and CMS results do not depend on the specifics
of the earlier Jupiter model and there is no need to update the $\lambda$ grid
iteratively.

At the beginning of the main loop in Fig. 2, we use Newton’s method to update
the shapes of all spheroids $\zeta_{j,k}$ (Hubbard, 2013), which is the most
time consuming step besides the wind calculations. If the shapes are far from
equilibrium, it may happen that the updated shapes of two adjacent spheroids
cross, which is unphysical. We detect such rare events and correct spheroid
shapes. It may also happen that the shapes start oscillating between two
states. When we detect this, we insert a regula falsi step.

We perform the thermal wind calculations once in the main loop to update the
target value for $J_{2}=J_{2}^{\rm Juno}-J_{2}^{\rm wind}$. In the fast, inner
loop we update the gravity coefficients, $J_{n}$, spheroid densities
$\rho_{j}$ and pressures, $P_{j}$, because these calculations are inexpensive
and scale linearly with spheroid number. We also update $Z_{1}$, $Z_{2}$, and
$Y_{2}$ as described above.

If the model includes a compact core of prescribed mass, some additional steps
are needed. First we derive a scale factor, $f=(m_{\rm core}^{\rm
current}/m_{\rm core}^{\rm target})^{-1/3}$ and then scale the $\lambda$ grid
inside the core region with $\lambda_{j}=f\times r_{\rm
core}\times(N_{l}-j)/N_{l}^{\rm core}$ so that the equatorial radius of the
outermost core spheroid becomes $r_{\rm core}^{\rm(new)}=f\times r_{\rm
core}^{\rm(old)}$. $N_{l}^{\rm core}$ is the number of spheroids in the core
while $N_{l}$ is the total spheroid number. We also update the logarithmic
$\lambda$ grid in the envelope so that it ends at the revised $r_{\rm
core}^{\rm(new)}$ value.

We find that our model predictions are not very sensitive to the assumed
composition of a compact core because the density of typical core materials is
much higher than that of hydrogen and helium, so it is primarily the mass of a
compact core that matters. We still see some deviations when we compare models
with a rocky compact core with models that assume a 50:50 rock-ice mixture. We
follow Seager et al. (2007) when we adopt equations of state for the core
materials but then complement it with more results from more recent ab initio
simulations (Militzer & Wilson, 2010; Wilson & Militzer, 2014; Gonzalez-
Cataldo & Militzer, 2023). For rock, we assume a terrestrial mixture of MgSiO3
and iron with an iron mass fraction of 32.5%. We represent the ice in the
compact core by assuming a pure H2O composition.

###  2.2 Assumptions for Model Ensembles

We employ Monte Carlo methods (Militzer, 2023a, b) to construct ensembles of
Jupiter models by accepting and rejecting moves according to the
$\exp(-\chi^{2}/2)$ function that combines the terms,
$\chi^{2}=\chi^{2}_{J}+\chi^{2}_{\rm H-He}+\chi^{2}_{\rm wind}+\chi^{2}_{\rm
guide}$. The first and most important term measures the deviations of even and
odd gravity harmonics between model predictions and the Juno measurements
(Durante et al., 2020),

$\chi^{2}_{J}=\sum\limits_{i=1}^{10}\left[\frac{J_{i}^{\rm model}-J_{i}^{\rm
Juno}}{\delta J_{i}^{\rm Juno}}\right]^{2}\quad.$ (15)

$\delta J_{i}^{\rm Juno}$ are the 1-$\sigma$ uncertainties of the
measurements. Following Militzer et al. (2019), we added the term
$\chi^{2}_{\rm H-He}$ to generate models with helium rain parameters $P_{\rm
rain,1}$ and $P_{\rm rain,2}$ that are compatible with phase diagram of H-He
mixtures as derived by Morales et al. (2013). From the assumed entropy values
for the molecular and metallic layers, $S_{1}$ and $S_{2}$, one can infer the
temperatures $T_{1}=T(S_{1},P_{\rm rain,1})$ and $T_{2}=T(S_{2},P_{\rm
rain,2})$ from the isentrope that are given by the equation of state. For the
pairs $P_{\rm rain,1}$-$T_{1}$ and $P_{\rm rain,2}$-$T_{2}$, we find the
closest points on the immiscibility curve by Morales et al. (2013),
$P^{*}_{1}$-$T^{*}_{1}$ and $P^{*}_{2}$-$T^{*}_{2}$, and then define the
immiscibility penalty term,

$\chi^{2}_{\rm
H-He}=\sum\limits_{i=1}^{2}C_{P}\left|\frac{P^{*}_{i}-P_{i}}{P_{i}}\right|+C_{T}\left|\frac{T^{*}_{i}-T_{i}}{T_{i}}\right|\quad,$
(16)

where $C_{P}$ and $C_{T}$ are weights that need to be balanced with those in
other $\chi^{2}$ terms. We set $C_{T}/C_{P}=2$. We do not square the
individual terms because there is currently no agreement between ab initio
simulation and experimental prediction on the conditions where hydrogen and
helium become immiscible. Early path integral simulations showed that the two
materials are miscible at 15 000 K (Militzer, 2005). Vorberger et al. (2007)
showed with ab initio simulation that hydrogen and helium are miscible at 8000
K. With more careful ab initio Gibbs free energy calculations, Morales et al.
(2013) predicted hydrogen and helium to phase separate at approximately 6500 K
for a pressure of 150 GPa. However, recent shock wave experiments by Brygoo et
al. (2021) placed the onset of immiscibility at a much higher temperature of
10 200 K for 150 GPa. Because the deviations of the ab initio predictions are
large and these findings to not yet been reproduced with other laboratory
measurements, we will employ the Morales et al. (2013) results when we
evaluate the $\chi^{2}_{\rm H-He}$ term in Eq. 16.

We also add the wind term,

$\chi^{2}_{\rm
wind}=\frac{1}{m}\sum\limits_{i=1}^{m}\left\\{\begin{matrix}\left[H(\mu_{i})-H_{\rm
max}\right]^{2}&{\rm if\,\,}H(\mu_{i})>H_{\rm max}\\\
0&\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm if\,\,}H_{\rm min}\leq H(\mu_{i})\leq H_{\rm
max}\\\ \left[H_{\rm min}-H(\mu_{i})\right]^{2}&{\rm if\,\,}H(\mu_{i})<H_{\rm
min}\\\ \end{matrix}\right.\quad,$ (17)

that keeps the depth of our winds, $H$, within perscribed limits of $H_{\rm
min}=1500\,$km and $H_{\rm max}=4500\,$km so that the remain broadly
compatible with predictions by Guillot et al. (2018). We evaluate them at
$m=61$ equally spaced $\mu$ points between –1 and +1 with $\mu=\cos(\theta)$
and $\theta$ being the colatitude. We work directly with the observed cloud-
level winds from Tollefson et al. (2017) but then make the wind depth to be
latitude dependent. One may allow the surface winds to deviate from the
observed value and employ the same wind depth for all latitudes. Both types of
wind models are compared in Militzer et al. (2022).

We solve the thermal wind equation (Kaspi et al., 2016) to derive the density
perturbation, $\rho^{\prime}$,

$\frac{\partial\rho^{\prime}}{\partial
s}=\frac{2\omega}{g}\frac{\partial}{\partial z}\left[\rho u\right]\quad,$ (18)

for a rotating, oblate planet (Cao & Stevenson, 2017). Instead of
approximating the planet’s shape by a sphere, we employ the geometry that is
provided by the equi-potential surfaces from the CMS calculation. $z$ is the
vertical coordinate that is parallel to the axis of rotation. $s$ is the
length of a path that starts form the equatorial plane and follows an
equipotential surface that we derived with the CMS method. This method also
provides us with the static background density, $\rho$, and the local
acceleration, $g$, for a particular model. $u$ is the differential flow
velocity with respect to the uniform rotation rate, $\omega$. We represent $u$
as a product of the surface winds, $u_{s}$, from Tollefson et al. (2017) and a
decay function of $\sin^{2}(x)$ form (Militzer et al., 2019). This function
facilitates a rather sharp drop similar to functions employed in Galanti &
Kaspi (2021) and Dietrich et al. (2021). We integrate the density
perturbation, $\rho^{\prime}$, to determine the dynamic contributions to the
gravity harmonics before we combine them with the gravity harmonics from the
interior, $J_{n}^{\rm model}=J_{n}^{\rm interior}+J_{n}^{\rm wind}$ and enter
them into Eq. 15.

Finally we add a number of additional penalty terms,

$\chi^{2}_{\rm guide}=C\left\\{\begin{matrix}\left[p_{\rm
min}-p\right]^{2}&{\rm if\,\,}p<p_{\rm min}\\\ 0&{\rm otherwise}\\\
\end{matrix}\right.\quad,$ (19)

that help us guide the Monte Carlo ensemble to reach and remain in regions
where a certain parameter, $p$, satisfies the condition $p\geq p_{\rm min}$
that we consider physical. Similar terms can assure $p_{2}\geq p_{1}$ for two
model parameters. We choose a large value of $C$ like 1000 to assure
compliance. We verify that $\chi^{2}_{\rm guide}=0$ for models that we
publish. $Z_{1}\geq Z_{\rm protosolar}$ is an obvious condition to satisfy but
we also require $Z_{2}\geq Z_{1}$ and $S_{2}\geq S_{1}$.

##  3 Results and Discussion

###  3.1 Convergence tests

Table 1: Parameters of different interior models in Fig. 1.

Parameter | 2 layer | 3 layer | 4 layers | 4 layers | 5 layers | 5 layers | 5 layers | 5 layers | 6 layer
---|---|---|---|---|---|---|---|---|---
| model | model | type A | type B | type A | type B | type C | type D | model
$Z_{1}$ | 0.0206 | 0.0138 | 0.0161 | 0.0166 | 0.0156 | 0.0153 | 0.0153 | 0.0153 | 0.0160
$\alpha$ | - | - | 3.00 | - | 9.41 | 11.18 | - | 10.07 | 13.68
$P_{\rm rain,1}$ (GPa) | 900.00 | 550.00 | 93.32 | 442.93 | 93.08 | 91.64 | 515.00 | 91.81 | 93.09
$P_{\rm rain,2}$ (GPa) | 900.00 | 550.00 | 760.59 | 442.93 | 443.17 | 613.90 | 515.00 | 636.65 | 609.04
$P_{\rm core,1}$ (GPa) | 900.00 | 550.00 | 1189.65 | 760.90 | 786.17 | 613.90 | 622.18 | 736.75 | 731.92
$P_{\rm core,2}$ (GPa) | 900.00 | 550.00 | 1189.65 | 2016.15 | 2051.67 | 903.61 | 905.54 | 736.75 | 916.33
$Z_{2}$ | 0.1332 | 0.0931 | 0.1639 | 0.1780 | 0.1831 | 0.1095 | 0.1119 | 0.1083 | 0.1178
$M_{\rm core}^{\rm compact}(M_{E})$ | 0.0 | 4.0 | 0.0 | 0.0 | 0.0 | 2.5 | 2.5 | 2.5 | 2.0

Note. — We underlined repeated pressure values because they indicate that a
specific layer is not included in a model. Machine readable data files for
every model are available online Militzer (2024).

Figure 3:  Convergence of the gravity harmonics $J_{4}$ and $J_{6}$ with
increasing number of spheroids is demonstrated for selected 2, 3, and 5 layer
models that match mass, radius, and $J_{2}$ as measured by Juno or a slightly
modified value to account for wind effects. The 3 layer models have compact
cores of 4 Earth masses with and without an ice fraction. The 5 layer models
do not include any density discontinuities and thus converge much faster than
the 2 and 3 layer models. A comparison of the 5 layer models shows that a
correction of $J_{2}$ by 10-6 has a small effect on the computed $J_{4}$ and
$J_{6}$ values. The Juno measurements have been subtracted from both Y axes.
The results of the 2 layer models were shifted to reduce the range of Y axis.
All models were constructed without winds.

In Fig. 3, we compare how rapidly selected models of Jupiter’s interior
converge in $J_{4}$-$J_{6}$ space with increasing number of spheroids. All
models were constructed to agree with the protosolar helium abundance and to
exactly match Jupiter’s mass, radius, and, unless noted otherwise, also
$J_{2}$ value as measured by the Juno spacecraft. The 5 layer models converge
the fastest because they do not contain any discontinuities in density.
Already with 256 spheroids, the computed $J_{6}$ value differs only by
$2\times 10^{-10}$ from corresponding result with 65536 spheroids, a
difference that is smaller than Juno’s 1-$\sigma$ error bar of $2\times
10^{-9}$ (Durante et al., 2020).

Converging $J_{4}$ is a bit more challenging. With 1024 and 2048 spheroids,
$J_{4}$ differs by $2\times 10^{-9}$ and $3\times 10^{-10}$ from the
65536-spheroid result while the Juno $J_{4}$ error bar is $8\times 10^{-10}$.
In all following figures, we report results that were obtained with 2048
spheroids that we derived with the accelerated version of CMS method (Militzer
et al., 2019) using an acceleration factor of $n_{\rm int}=32$. For the models
in Fig. 3, we compared the results derived with $n_{\rm int}=32$ and those of
the not-accelerated version ($n_{\rm int}=1$) but the differences were much
too small to be discernible in Fig. 3. This underlines that the convergence of
$J_{4}$ and $J_{6}$ is foremost controlled by the spheroid number and model
type instead of the acceleration factor.

When the coefficients of Jupiter’s gravity field are reported (Durante et al.,
2020), the effects of the satellites have been removed from the forces acting
on the spacecraft. This includes not only their direct gravitational forces
but also the effect of the permanent tides that they introduce in Jupiter via
the Love number $k_{20}$. The measured $J_{2}$ was thus reduced by $6.72\times
10^{-8}$, a value that was derived with theoretical methods because this
effect cannot be measured directly. In Fig. 3, we include two 5 layer models
with altered $J_{2}$ values in order to determine which alteration magnitude
is required to significantly affect the $J_{4}$ and $J_{6}$ that we calculated
with a specific model. The results in this graph illustrate that $J_{2}$
modifications of $10^{-6}$ or less are negligible for the purposes of this
article, which means the tidal correction of $J_{2}$ is approximately one
order of magnitude too small to matter for our interior models.

In Fig. 3, we also show 2 and 3 layer models that have an abrupt change in
composition in the envelope. Because the density is no longer a continuous
function of pressure, many more spheroids are needed to converge $J_{4}$ and
$J_{6}$ with high precision. With 2048 spheroids, we are able to compute
$J_{4}$ and $J_{6}$ with a precision of $5\times 10^{-8}$ and $6\times
10^{-9}$ or better, respectively. Both values are sufficiently small for the
purpose of this study but the convergence could be improved with a more
careful treatment of the layer interfaces.

The magnitudes of $J_{4}$ and $J_{6}$ of our 2 layer model in Fig. 3 are
smaller than those of our 3 layer models. The 2 layer models better mimic the
dilute-core effects in Militzer et al. (2022) than the 3 layer models can
because they have 4 Earth masses stored in a compact core, which reduces the
dilute-core effect. The magnitudes of $J_{4}$ and $J_{6}$ are slightly larger
for an ice-free, rocky core composition than for a 50:50 rocky-ice mixture.
Ice-free, rocky cores have higher density and the surrounding H-He fluid is
thus exposed to a higher pressure, which increases its density and thereby
enhances the effect of the core. This is the same trend that we have seen when
modeling Saturn’s interior (Militzer et al., 2019). The masses of rocky cores
were found to be slightly smaller than those of rock-ice cores.

###  3.2 Two and three layer models

Figure 4:  Predictions from two and three layer models that have a sharp
transition in the planet’s envelope at a specified transition pressure (x axis
of upper diagram). Three layer models have a compact core of 2, 4, or 6 Earth
masses (blue, orange, black lines) that is either of rocky (solid lines) and
rock-icy (dash-dotted lines) composition. The green symbols and band indicate
the range of models with $J_{4}$ values equal to $J_{4}^{\rm Juno}\pm
0.7\times 10^{6}$ across all panels. The blue filled pentagons and circles
indicate models with transition pressures, $P_{T}$, of 500 and 800 GPa,
respectively. The arrow in the lower panel illustrates the shift in
$J_{4}$-$J_{6}$ space that occurs when a 6 $M_{E}$ compact core is added a
model with $P_{T}$=800 GPa. The open and filled triangles show models with
$Z_{1}$ values of protosolar or twice protosolar composition. The dark grey
region indicates negative and thus unphysical abundances, $Z_{1}<0$. The
medium grey regions show where $Z_{1}$ or $Z_{2}$ are less than the protosolar
value (0.0153 (Lodders, 2010)). The shaded triangular region in the middle
panel indicates models that are Ledoux unstable ($Z_{2}<Z_{1}$) (Ledoux,
1947). For reference, the yellow circle in the lower panel shows the Juno
measurements that are exactly reproduced by our preferred five layer models
that includes winds. Its prediction from the interior alone are indicated by
the yellow star.

The results from our two and three layer models are summarized in Fig. 4. All
models match mass, radius, $J_{2}$, and the two conditions for helium
abundances $Y_{1}$ and $Y_{2}$ so that Galileo helium measurement is
reproduced and the planet overall has a protosolar helium fraction. All models
have a sharp transition in the envelope at the prescribed transition pressure,
which is their main free parameter. For the three layer models, one can also
prescribe the mass of compact core, which can either have a rocky composition
with a terrestrial iron fraction of 32.5% or a rock-ice core with a Callisto-
type ice fraction of 50%. Switching between these two plausible core
compositions has only a modest effect on the other inferred properties. The
two upper panels of Fig. 4 show that, for given core mass and transition
pressure, ice-free cores increase $Z_{1}$ and decrease $Z_{2}$ slightly more
than the corresponding rock-ice cores if one compares the predictions with the
corresponding core-less two layer models. Also we found that a rock-ice core
can have a 13% larger mass than a ice-free core to yield the same shift in
$J_{4}$-$J_{6}$ space, which is consistent with predictions for Saturn
(Militzer et al., 2019) and is driven by the fact that the density contrast
between ice and H-He is less than that between rock and H-He.

The upper two panels of Fig. 4 shows how $Z_{1}$ and $Z_{2}$ respond to
variations in the transition pressure. We added symbols and shaded regions to
make it easier to follow different trends. The upper panel illustrates that
all two and three layer models require a unexpectedly large transition
pressure $\gtrsim 400$ GPa for $Z_{1}$ to be at least of solar abundance.

The upper panel also shows that, if models are required to match mass, radius,
and $J_{2}$ but not any $J_{n\geq 4}$, an increase in the core mass leads to
an increase in $Z_{1}$. At the same time, the middle panel shows that such an
increase of the core mass is accompanied by a reduction $Z_{2}$ so that the
mass in the central region of the planet is approximately preserved. Besides
these two main trends, the shape of the curves in the middle panel is
controlled by two competing effects. If we increase the transition pressure,
the mass of inner layer shrinks, so its heavy element fraction, $Z_{2}$, may
increase, as red curves for core-less models show. Conversely if the inner
region already has a sizable compact core ($\gtrsim$ 3 Earth masses), there is
less room for heavy elements in the metallic region and $Z_{2}$ shrinks with
increasing $Z_{1}$ (orange and black curves in the middle panel).

An important constraint on the core mass comes from the Juno measurements of
$J_{4}$ and $J_{6}$, which we illustrate in the lower panel of Fig. 4. The
green symbols and band represent models in the corridor around the Juno’s
measurements $J_{4}$ of size $\pm 0.7\times 10^{6}$ that was chosen because it
is difficult to construct wind models that yield a larger shift in $J_{4}$. We
replicate the green band in the upper two panels to further constrain the core
mass of plausible models. If one combines the $J_{4}$ interval with the
requirement that the outer envelope should have at least a protosolar
abundance of heavy elements ($Z_{0}=0.0153$ (Lodders, 2010)), models with
cores larger than 4 Earth masses can be ruled out. (In upper two panels, there
is no overlap between black curves and the green band in a region that has not
been greyed out.) The conditions on $J_{4}$ (green band) and $Z_{1}\geq Z_{0}$
also require the transition pressure to be 570 GPa or larger, as the upper
panel of Fig. 4 illustrates. Such a high transition pressure is not consistent
with the physics of hydrogen-helium mixtures. Ab initio computer simulations
have place the onset of the hydrogen-helium immiscibility region at
approximately 100 GPa (Morales et al., 2013). The discrepancy in pressure
makes all three layer models appear less plausible unless one increases
Jupiter’s interior temperature (Miguel et al., 2022) or changes the EOS to
reduce the density of H-He mixtures (Nettelmann et al., 2021). When Miguel et
al. (2022) employed the Militzer & Hubbard (2013) EOS, most models include a
helium transition from 300 to 450 GPa (see their Fig. A1). The trend towards
high values is the same that we find here but the pressure values in the work
by Miguel et al. (2022) are not as high as ours because they allow the 1 bar
temperature to vary and predict values of approximately 183 K. Miguel et al.
(2022) also reported a second, smaller group of interior models with a
transition pressure of approximately 200 GPa. We do not find models because we
do not include large jumps in temperature of up to 2000 K at the transition
pressure.

The bottom panel of Fig. 4 shows that the inclusion of 6 Earth mass cores
leads to large shifts in $J_{4}$-$J_{6}$ space away from the Juno measurements
because it diminishes the dilute-core effect that enabled us to match the Juno
measurements in the first place (Militzer et al., 2022). If one compares
rocky-core models without winds and a transition pressure of 800 GPa (blue
filled circles), the magnitude of the shift is $\Delta J_{4}\approx-8.5\times
10^{-6}$ and $\Delta J_{6}\approx+0.76\times 10^{-6}$ (see grey arrow). For a
lower transition pressure of 500 GPa (blue filled hexagons), the shift is
still half as large approximately because the dilute core region contains more
mass and the addition of compact core of 6 Earth masses does not have such a
large effect. Still it causes the predicted $J_{4}$ values to fall outside of
the green corridor where one could argue that wind models would to able to
bridge the remaining difference from the Juno measurements.

###  3.3 Four, five and six layer models

Figure 5:  Distributions of heavy elements and helium of models in Fig. 1 as
function of fractional cumulative mass that has been derived by adding the
contribution of different spheroids from the inside out. All models have
dilute cores that either end abruptly [models with two, three, and four (A)
layers] or gradually [four (B), five (A+B) and six layers]. The three, five
(B), and six layer models also have compact cores ($Z=1$) of 4.0, 2.5, and 2.0
Earth masses. In the outer envelope, all models match the helium abundance
that the Galileo entry probe measured. In the interior, they transition to a
layer with a higher helium abundance so that the planet has a protosolar
hydrogen-helium abundance overall. This transition may either be abrupt [two,
three, four (B) layer models] or occur over a extended helium rain region
[four (A), five (A+B) and six layer models].  Figure 6:  Temperature, density,
and pressure are plotted vs. fractional mass. For clarity only selected models
are shown. To better show the differences between models, we scale pressure by
density-2 because this function is a constant for a polytrope EOS of index 1
(Howard et al., 2023). This scaling is appropriate for the dense part of H-He
layers but not for the atmosphere nor inside of a compact core. With increase
mass (or radius), the two layer model includes a large decrease in $Z$ at the
layer boundary, which decreases its density (middle panel), increases the
$P/\rho^{2}$ term (lower panel) despite an decrease in temperature (upper
panel). The $P/\rho^{2}$ term of our five layer model increases throughout the
core transition layer ($M=0.2-0.5$) because $Z$ increases there. The
$P/\rho^{2}$ term is approximately constant inside the dilute core because it
is of constant composition. For comparison, we show an isentrope of protosolar
composition ($Y_{0}=0.274$, $Z_{0}=0.0153$) that starts from a 1 bar
temperature of 166.1 K.  Figure 7:  Predictions for the even $J_{n}$ from
different interior models and their wind contributions are plotted in
$J_{n}$-$J_{n+2}$ spaces. The Juno measurements are shown by the large red
circles with their 3-$\sigma$ error bars as reported in Durante et al. (2020).
The other symbols represent the seven models from Fig. 1. The large symbols
represent the contributions from the interior model only while the small
symbols represent interior and wind contributions combined. The combined
contributions fit the Juno measurements very well except for the two layer
model that fails to fit $J_{4}$.

In this section, we compare properties of the Jupiter models in Fig. 1 that
have at least four layers. The corresponding distributions of helium and heavy
elements are shown in Fig. 5. Temperature, density, and pressure are shown in
Fig. 6. The temperature profile is determined by an outer isentrope with
entropy, $S_{1}$, where hydrogen occur molecular form, and an inner isentrope
with $S_{2}$, and where hydrogen is metallic. See section 2.2 and Militzer et
al. (2022) for details. The bottom panel of Fig. 6 includes an isentrope of
protosolar composition for comparison. At lower pressures and high fractional
mass ($M>0.75$), its density is higher than that of all models because they
include a depletion of helium ($Y<Y_{0}$). At higher pressure and low fraction
mass ($M<0.5$), the protosolar density is lower because the models include
dilute cores ($Z\gg Z_{0}$).

All calculations including wind models were generated under consistent
assumptions (Militzer et al., 2022). The contributions from the interior and
winds are illustrated in Fig. 7. Except for our two layer models, all other
models can fit the Juno gravity measurements very well. Therefore the
following discussion will focus on the question of which models are more
plausible given the physics of hydrogen-helium mixtures and, to a lesser
degree, on magnetic field information.

With our two layer models, we were not able to match Jupiter’s $J_{4}$ and
$J_{6}$. Even when wind contributions were included, discrepancies of $\Delta
J_{4}\approx 3\times 10^{-7}$ and $\Delta J_{6}\approx 2\times 10^{-8}$
remained (see Fig. 7). For this reason, we do not discuss the two layer models
further. Still, the best models in this class have an abrupt change in
composition at $\sim$900 GPa and fractional radius $\sim$0.6.

Figure 8:  The symbols show the cost function $\chi^{2}$ of model ensembles
that were constructed with compact cores of a chosen mass under four different
assumptions. We added a rocky core to our reference five layer model (the red
symbols). The red line and shaded area were added to illustrate that for core
masses larger than 2.5 Earth masses, the cost function increases drastically,
and models that match the Juno data well can no longer be constructed. If
compact cores of a rock-ice composition are introduced instead, cores of 3.0
Earth masses still yield good models (black dashed lines). If the density is
lowered by 2% over the pressure range from 20 to 100 GPa, compact cores of
$\sim$5 Earth masses become possible (blue symbols, line, and shaded area). If
the density is lowered by 3% over the pressure range from 10 to 100 GPa,
compact cores of up to $\sim$6 Earth masses still match the Juno data quite
well (green symbols, line, and shaded area).

We begin the following discussion with our fiducial five layer model of type A
that include a helium-rain layer and a dilute core where the heavy elements
only contribute mass fraction of $Z_{2}\sim 0.18$ (see Fig. 5). Their fraction
decreases to $Z_{1}\approx 0.0153$ across the core transition layer where the
pressure decreases from 903 to 614 GPa and the fraction cumulative mass
increase from 0.20 to 0.54. By inserting a compact core with $Z=1$ into this
model, we obtain our six layer models. We are not able to insert a rocky core
heavier than 2.5 Earth masses without deviating from the Juno gravity
measurements as we illustrate in Fig. 8. If we assume a rock-ice rather than a
purely rocky composition for the compact core, the maximum mass increases to
3.0 Earth masses (see Fig. 8). Here we assumed the compact core to be
homogeneously mixed and to have a Callisto-type composition with 50%-50% rock-
ice mass fraction. Ab initio EOS tables were used to derive the density for
all materials in the compact core (Wilson & Militzer, 2014). The ice fraction
in the compact core reduces its density and therefore decreases the density
contrast to the hydrogen-helium mixture, which means the models can
accommodate slightly larger compact cores.

It is possible to accommodate larger compact cores but this would require a
change in the EOS of hydrogen and helium. If their density were reduced by 2%
below values of the Militzer & Hubbard (2013) EOS in the crucial pressure
interval from 20 to 100 GPa, a rocky compact core of $\sim$5 Earth masses
could be accommodated. If the density from 10 to 100 GPa is reduced by 3%
compact cores of $\sim$6 Earth masses become plausible as we illustrated in
Fig. 8.

The reason for a mass limit on the compact core is that it takes away from the
dilute-core effect that enabled us to match $J_{4}$ and $J_{6}$ in the first
place. Furthermore, Fig. 1 and 5 show that a compact core of 2.5 Earth masses
leads to a reduction of heavy elements in our dilute core from $Z_{2}\sim$18%,
to $\sim$11%. More importantly the presence of a compact core increases the
dilute core radius to $\sim$0.61 where it decays rather sharply until
$R\sim$0.65. It thereby “pushes” against the lower boundary of the helium rain
layer and thereby shrinks the metallic hydrogen layer. It is consistent with
the gravity measurements to eliminate the metallic hydrogen layer entirely,
which yields the five layer model of type B in Fig. 1. Alternatively one can
keep the metallic hydrogen layer and eliminate the core transition layer to
obtain the five layer model of type D. Model types B and D have extended
helium rain layers and a compact core, which make them similar to our six
layer models. Finally, one may remove the helium rain layer from our six layer
models to obtain our five layer models of type C. To match Juno’s $J_{4}$ and
$J_{6}$ measurements, the transition to metallic hydrogen must occur at a
rather high pressure of $\sim$500 GPa, which makes this model similar to our
four layer models of type B.

Additional information might be harnessed by analyzing Jupiter’s magnetic
field which was found to be unexpectedly complex (Moore et al., 2018). Recent
magnetohydrodynamic simulations compared the effects of stably stratified (not
convective) and fully convective layers in Jupiter’s interior and concluded
that a stably stratified interior structure would not be compatible with the
observed magnetic field (Moore et al., 2022). This favors models with thick
layers of metallic hydrogen and extended dilute cores because we assume these
regions to be homogeneous and convective. Still given the magnetic field data,
it is not possible to rule out any class of models rigorously.

There are two obvious choices for generating four layer models from our five
layer model of type A. Both are illustrated in Fig. 1. One can either
eliminate the core transition layer (model type A) or the helium rain layer
(model type B). In A-type models, the dilute core is terminated abruptly
($P_{\rm core,1}=P_{\rm core,2}$) while it decays gradually across the core
transition layer of our five layer models. B-type models have a sharp boundary
in composition at $P_{\rm rain,1}=P_{\rm rain,2}$ between the layer that is
dominated by molecular hydrogen and the layer below that contains mostly
metallic hydrogen.

Figure 9:  Correlation plot for parameters of our four and five layer models.
The four layer models of type B and five layer models of type C have no helium
rain layer, which implies $P_{\rm rain,1}=P_{\rm rain,2}$ and leads to the
straight diagonal line in panel (a). The other model types plot farther to
left in this panel because they have extended helium rain layers with $P_{\rm
rain,2}\gg P_{\rm rain,1}$. The four layer models of type A and five layer
models of type D have no core transition layer, which implies $P_{\rm
core,1}=P_{\rm core,2}$ and means that those models plot on a straight line in
panel (b). The other model types have extended core transition layers with
$P_{\rm core,2}>P_{\rm core,1}$. Panels (d) and (e) show that the heavy
element fraction of the core, $Z_{2}$, is positively correlated with the outer
pressures and in some cases also with the inner pressure of the core
transition layer. Panels (c) and (f) show that the total mount of heavy
elements varies between 23.5 and 25.5 Earth masses.

In Fig. 9, we plot correlation among various parameters of our four and five
layer models. Our four layer models of type B and five layer models of type C
have a sharp instead of gradual change in helium abundance. We find that this
change in composition must occur at a rather high pressure from 350 to 550
GPa. This is much higher than a pressure of $\sim$100 GPa for which ab initio
simulations have predicted of onset of the hydrogen-helium immiscibility
Morales et al. (2013). For that reason, we do not favor Jupiter models that do
not include an extended helium rain layer.

Fig. 9 also shows that the helium rain layers of our four layer models of type
A and our five layer models of types B and C have more extended helium rain
layer than our five layer models of type A (see also Fig. 1). The $P_{\rm
core,1}$-$P_{\rm core,2}$ correlation plot shows that our five layer models of
type A and four layer models of type B have thick core transition layers that
may extend approximately from 650 to 2400 GPa. In our four layer models of
type A, the dilute core ends abruptly at $P_{\rm core,1}=P_{\rm core,2}$ from
1000 to 1300 GPa. In our five layer models of type D, this sharp transition
occurs at lower pressures from 650 to 800 GPa. Our five layer models of types
B and D are not very differ but allow for larger $P_{\rm core,2}$ values of up
to 1100 GPa.

Furthermore Fig. 9 shows that $Z_{2}$, the heavy elements abundance of the
dilute core, is positively correlated with $P_{\rm core,2}$ and in some cases
also with $P_{\rm core,1}$. The reason is that for fixed $Z_{2}$ value with
$Z_{2}\gg Z_{1}$, an increase in $P_{\rm core,2}$ shrinks the dilute core and
this leads to a reduction of heavy element abundance of the core region
overall (see Fig. 5). To compensate for such a reduction, $Z_{2}$ need to rise
if $P_{\rm core,2}$ is increased. One also finds that five layer models of
types B, C, and D have lower $Z_{2}$ values because they have larger dilute
cores.

Finally panels (c) and (f) of Fig. 9 show that the total amount of heavy Z
elements, $M_{Z}$, in the planet varies between 23.5 and 25.5 Earth masses
(7.4–8.0%). Within the ensemble of a given model type, $M_{Z}$ is not strongly
correlated with $P_{\rm core,1}$ nor with $Z_{2}$ but it does depend on the
model type. Five layer models of types B, C, and D, that include a compact
core of 2.5 Earth masses, have a slightly lower total amount of heavy elements
of between 23.5 and 24.5 Earth masses. Four layer models of type A, that do
not have a core transition layer, have between 24 and 25 Earth masses worth of
heavy elements. The highest amounts of between 24.8 and 25.4 Earth mass are
predicted for five layer models of type A and for four layer models of type B.
Similar results can be expected because both model types differ only by the
helium rain layer. It should be noted, however, that one can construct models
with substantially higher amounts of heavy Z elements by assuming a
superadiabatic temperature exists in the core transition layer (Militzer et
al., 2022).

Figure 10:  Cost function $\chi^{2}$ vs. helium rain parameter, $\alpha$, in
Eq. 2 is shown for an ensemble of five layer models (blue diamonds) and one of
four layer models of type A with a sharp core boundary (red circles). With
increasing $\alpha$, more helium has been sequestered to deeper layers, which
reduces the density in the upper region of the helium rain layer and thereby
makes it easier to construct models that match the Juno and Galileo data. For
$\alpha$ less than 3, a good match becomes increasingly difficult to obtain,
which manifests itself in a sharp rise in the cost function for small $\alpha$
values. This increase is very similar for both type of ensembles as two shaded
regions illustrate.

Figure 10 was constructed to study whether Juno gravity measurements allow us
to constrain the helium distribution within the helium rain layer. For this
purpose, we constructed various MC ensembles for different exponents, $\alpha$
in Eq. 2, for four and five layer models of type A. We find that a linear
increase of helium fraction with $\log(P)$ is incompatible with the Juno
measurements and our modeling assumptions. Figure 10 shows that $\alpha$ needs
to be three or larger, which means a substantial amount of helium has been
sequestered from the upper part of the helium rain layer.

On the other hand, we were not able to obtain an upper limit of $\alpha$. With
increasing $\alpha$, the composition within the helium rain approaches that of
the molecular layer above and transition to the metallic layer becomes
increasingly sharp. In the limit of large $\alpha$, the five layer models of
type A approach the four layer models of type B that require a deep transition
($R\sim 0.71$) from 350 to 550 GPa.

We see the same trend in our three layer models that include a compact core
and a sharp transition from molecular to metallic hydrogen. This transition
needs to be at $\sim$550 GPa to match the gravity data. This is a pressure
value that is too high to be compatible with the physics of hydrogen-helium
mixtures. So in conclusion, we do not favor two layer, three layer, B-type
four layer models. On the other hand, four layer models of type A, five layer
models of types A, B and D but also six layer models with a very modest
compact core remain plausible.

###  3.4 $J_{n}$ weight functions

Equation 13 specifies which contribution a particular spheroids makes to the
gravity harmonics, $J_{n}$. Therefore, the individual terms,
$\lambda_{j}^{n}\tilde{J}_{j,n}=\delta_{j}\hat{J}_{j,n}/M$, have been
interpreted as weights, which have been plotted as function of spheroid radius
to illustrate that the most important contributions to higher order gravity
harmonics come from the outer layers of the planet (Militzer et al., 2016).
Here, we propose a slightly different normalization factor of $1/N_{j}$ for
these functions to remove the dependence on the choice for the $\lambda$ grid,

$w^{\rm
org}_{j,n}=\frac{\lambda_{j}^{n}\tilde{J}_{j,n}}{N_{j,n}}=\frac{\delta_{j}\hat{J}_{j,n}}{MN_{j,n}}\;\;{\rm
with}\;\;N_{j,n}=J_{n}[\lambda_{j}-\lambda_{j+1}]\;.$ (20)

We also divide by $J_{n}$ to remove their magnitudes from the weight
functions. We show these weight functions in the upper panel of Fig. 11 for
our three and five layer models. The density discontinuity in the three layer
model at 550 GPa introduces a $\delta$ function into the weights at $R\sim
0.7$. From a CMS perspective, this is very plausible. The discontinuity in
density is represented by one particular spheroid, $j$, whose outer boundary
coincides with the 550 GPa level. This spheroid adds a particularly large
density, $\delta_{j}$, to the planet, which is equal to the jump in density at
the discontinuity from $\rho_{j-1}$ to $\rho_{j}$. Therefore, the weights in
Eq. 20 resemble in a way the first derivative of the planet’s interior density
structure with respective radius, which explains why the weight functions are
not particularly smooth for $R\gtrsim 0.9$.

When one interprets the weight functions in the upper panel of Fig. 11, one
might be tempted to associate the value of $w^{\rm org}(\lambda_{i})$ with the
material that is stored between $\lambda_{j}$ and $\lambda_{j+1}$ while in
fact it represents an entire spheroid that has an equatorial radius of
$\lambda$ and includes material in the very center of the planet. However, one
can easily derive the contributions to $J_{n}$ that arise from material that
is stored in between the spheroids $j$ and $j+1$ and we use this contribution
to define a modified weight function,

$w^{\rm
mod}_{j,n}=\frac{\rho_{j}}{MN_{j,n}}\left[\hat{J}_{j,n}-\hat{J}_{j+1,n}\right]\;.$
(21)

The original and modified weight functions are normalized in the same way,

$1=\sum_{j}\,[\lambda_{j}-\lambda_{j+1}]\,w^{\rm org/mod}_{j,n}\approx\int
d\lambda\,w^{\rm org/mod}_{n}(\lambda)\;.$ (22)

We plot these modified weight functions for our three and five layer models in
the middle panel of Fig. 11. They still convey the message that the outer
parts of the planet contribute most to the higher $J_{n}$ but now it is easier
to see why a dilute core, that extends to $R\sim 0.5$, has an impact on
$J_{4}$ and $J_{6}$. Furthermore the modified weight functions are much
smoother than the original ones. A density discontinuity at 550 GPa now only
leads to a discontinuity in the weight functions rather than a $\delta$
function. In simple terms one can say that the modified weight functions
reflect the contributions from the density in each layer while the original
weight functions represent the contributions from the density increase in
every layer and are therefore closer to core concept of the CMS technique that
represents a rotating planet by a sum of nested, constant-density spheroids.
The modified weight functions are similar to those that Fortney et al. (2016)
derived with theory of figures calculations but they are not identical, nor
are the weight functions from different theory of figures calculations because
assumptions for the interior density structure matter. Nettelmann et al.
(2013) predicted the weight functions of $J_{2}$–$J_{8}$ for Saturn to
respectively have maxima at 0.82, 0.90, 0.92, and 0.94 fractional radii while
Guillot (2005) predicted all weight functions for Jupiter to peak within a
narrow range of only 0.95–0.97. With CMS calculations for Jupiter, we predict
the four weight functions to have maxima at 0.68, 0.83, 0.88, and 0.91
fractional radii respectively.

As illustrated in Eq. 22, one can also view the sum over spheroids as an
integral over the equatorial radius $\lambda$. This allows us to rewrite Eq.
11 and ask the question how a small density modification,
$\rho^{\prime}(\lambda)$, changes the computed product $(MJ_{n})$,

$(MJ_{n})+(MJ_{n})^{\prime}=\int
d\lambda\,\left[\rho(\lambda)+\rho^{\prime}(\lambda)\right]\tilde{w}_{n}\equiv\left<\rho+\rho^{\prime},\tilde{w}_{n}\right>\;\;{\rm
with}\;\;\tilde{w}_{n}=\frac{\hat{J}_{j,n}-\hat{J}_{j+1,n}}{\lambda_{j}-\lambda_{j+1}}$
(23)

$\left<\ldots\right>$ represents the integral over $\lambda$. The density
change can be expressed in terms of $P$ basis functions,
$\rho^{\prime}(\lambda)=\sum_{m=1}^{P}c_{m}f_{m}(\lambda)$, which yields,

$(MJ_{n})^{\prime}=\left<\rho^{\prime},\tilde{w}_{n}\right>=\sum_{m=1}^{P}c_{m}\left<f_{m},\tilde{w}_{n}\right>$
(24)

This linear equation opens up the possibility of adjusting the coefficients,
$c_{m}$, to trigger a change of just one specific $J_{n}$ while leaving the
mass ($M=-J_{0}$) and other $J_{n}$ unchanged to first order. In the lower
panel of Fig. 11, we give an example for a density modification function that
changes $J_{4}$ but not the mass nor $J_{2}$.

This function was derived by setting $P=3$ to consider changes in M, $J_{2}$,
and $J_{4}$. We opted to introduce $P$ different density modification
functions, $\rho_{p}^{\prime}(\lambda)=\sum_{m=1}^{P}C_{p,m}f_{m}(\lambda)$
and to represent their coefficients by the matrix $C_{p,m}$. One can solve for
this matrix by inverting the matrix
$F_{m,n}\equiv\left<f_{m},\tilde{w}_{n}\right>$ so that every density
modification function affects just one of the three variables M, $J_{2}$, and
$J_{4}$. This requires the matrix product of $\overset{\leftrightarrow}{C}$
and $\overset{\leftrightarrow}{F}$ to equal the identity matrix,

$\overset{\leftrightarrow}{I}=\overset{\leftrightarrow}{C}\overset{\leftrightarrow}{F}$
(25)

The function in the lower panel of Fig. 11 has positive and negative regions
so that mass can be preserved when the density is modified. The function also
peaks at large radii where the magnitude of $J_{4}$ weight function is higher
than that of other functions. This is just one specific example. One needs to
keep in mind the shape of the resulting density modification functions
strongly depends on the choice of basis functions, $f_{m}$, and on the
variable $P$ that specifies how many $J_{n}$ are considered.

Figure 11:  The upper and middle panels show the original and modified weight
functions vs. the volumetric radius. The colors label various even $J_{2}$ –
$J_{8}$. The thick dashed lines show the weight functions for the three layer
model in Figs. 1, 5, and 7 while the thin dotted lines represent our fiducial
five layer model. The lower panel illustrates a density modification in the
planet that changes $J_{4}$ but leaves the planet’s mass and $J_{2}$ unchanged
to first order. (All curves remain unchanged regardless whether they were
computed with or without acceleration, $n_{\rm int}=32$ and $n_{\rm int}=1$.)

##  4 Conclusions

We constructed models for Jupiter’s interior with a dilute core that match the
recent gravity measurements of the Juno spacecraft as well as the temperature
measurements of the earlier Galileo entry probe. We employ the nonperturbative
CMS method to derive a hydrostatic solution (without winds) for the interior
of the rapidly rotating planet and then insert the oblate structure into the
thermal wind equation to derive the gravity contributions from the winds.
These selfconsistent models allow us matches even and odd harmonics of
Jupiter’s gravity field under one set of physical assumptions. For example, we
use equations of state derived from ab initio computer simulations to derive
the density of given pressure, temperature, and composition. We do not change
ab initio EOS in ad hoc ways, increase the temperature in Jupiter’s outer
layer nor construct models that have a less than solar abundance of the heavy
elements.

We studied models with different numbers of layers in order to determine which
models match the Juno gravity measurements and identify those that are most
compatible with the physics of hydrogen-helium mixture and information on how
the a giant planet convects. We found that one cannot match the Jupiter’s
gravity field with two layer models. A good match is obtained with models that
have three layers or more. However, three layer models and four layer models,
that have an abrupt transition from molecular to metallic hydrogen, require
that this transitions occurs very deep in the planet at $\sim 500$ GPa
respectively while experiments and ab initio computer simulations indicate the
molecular-to-metallic transition occurs approximately at 100 GPa. For this
reason, we find these models less plausible than our four, five and six layer
models that include an extended helium rain layer that begins at $\sim$100
GPa. This provides further support for the hypothesis that there are at least
two changes in composition in Jupiter’s interior: a helium rain layer and a
transition to a dilute core.

Moll et al. (2017) studied the time scale for giant planet cores to be eroded
convectively under two assumptions. They showed that if there exists only a
single sharp boundary in composition between core and the H-He envelope, the
core erodes very quickly on a million-year timescale. In this case, the
envelope efficiently extracts heat and heavy elements. However, if the core is
surrounded by a stably stratified layer, the core will be preserved over
billions of years. For this reason we favor our five layer models with a
stably stratified core transition layer over our four layer models of type A.
They have a sharp core-envelope boundary, which implies the dilute core would
be rapidly eroded.

Finally we derived our six layer models by inserting a compact core into our
extended dilute core of our five layer models. If we assume a rocky
composition for the compact core, we find it mass cannot be larger than 2.5
Earth masses. If we assume a rock-ice composition this limit increases to 3.0
Earth masses, which is only about 1% of the planet’s mass. Making the compact
core more massive than that would require removing so much mass from the
dilute core that one can no longer match the gravity harmonics $J_{4}$ and
$J_{6}$. Based on these arguments, we currently favor our five layer model but
a six layer structure with a small compact core of up to 3.0 Earth masses are
equally plausible. We expect such models to be revised when the equation of
state of hydrogen is measured more precisely, additional ab initio computer
simulations are performed, interior and magnetic field calculations are
coupled more tightly or when the June spacecraft provides us with additional
information for the planet’s atmosphere or gravitational field.

##  Appendix A N-Dimensional Root Finder

The Newton-Raphson method (Press et al., 2001) is an elegant procedure to
iteratively solve $N$ equations, $\vec{f}(\vec{x})=0$, with $N$ unknowns that
are represented by the vector $\vec{x}$. Starting from an initial guess, the
vector $\vec{x}$ is improved step by step,

$\vec{x}^{\rm(new)}=\vec{x}^{\rm(old)}+\Delta\vec{x}\;\;\;\;{\rm
with}\;\;\;\;\overleftrightarrow{J}\Delta\vec{x}=-\vec{f}\left(\vec{x}^{\rm(old)}\right)\;.$
(A1)

This methods requires knowledge of the Jacobian derivative matrix,
$J_{ij}=\frac{\partial f_{i}}{\partial x_{j}}$, which is often be inaccessible
in practical applications that involve complex functions $\vec{f}$. So here we
propose to obtain these derivatives by linearly fitting the last N+1 points,
$\vec{x}^{(1)},\ldots,\vec{x}^{(N+1)}$ and the every function $f_{i}(\vec{x})$
separately for $i=1\ldots N$. To simplify the following equations, we set
$N=2$, but our approach is general. For every function $i=1$ and 2 separately,
we solve this system of linear equations:

$\left[\begin{array}[]{cc}x_{1}^{(1)}-x_{1}^{(3)}&\;\;x_{2}^{(1)}-x_{2}^{(3)}\\\
x_{1}^{(2)}-x_{1}^{(3)}&\;\;x_{2}^{(2)}-x_{2}^{(3)}\end{array}\right]\left(\begin{array}[]{c}\frac{\partial
f_{i}}{\partial x_{1}}\\\ \frac{\partial f_{i}}{\partial
x_{2}}\end{array}\right)=\left(\begin{array}[]{c}f_{i}\left(\vec{x}^{(1)}\right)-f_{i}\left(\vec{x}^{(3)}\right)\\\
f_{i}\left(\vec{x}^{(2)}\right)-f_{i}\left(\vec{x}^{(3)}\right)\end{array}\right)$
(A2)

to obtain approximate values for the Jacobian, $J_{ij}$. Then we perform the
standard Newton-Raphson step in Eq. A1 to obtain the step $\Delta\vec{x}$ and
next point $\vec{x}^{(4)}=\vec{x}^{(3)}+\Delta\vec{x}$ before dropping the
earliest point, $\vec{x}^{(1)}$. To start, one needs to set a simplex of $N+1$
points that may be constructed from a single point, $\vec{x}^{(1)}$, by
shifting the $N$ coordinates individually,
$x_{j}^{(1+i)}=x_{j}^{(1)}+\delta_{ij}\Delta x_{j}$ where $\Delta x_{j}$
represent a chosen step size and $\delta_{ij}$ is the Kronecker delta symbol.
This root-finding algorithm relies in $N+1$ points while the Newton-Raphson
method uses only one. So it is similar to switching from the 1D Newton method
that uses one point $x$ and the derivative $\frac{df}{dx}$ to the regula falsi
method (Press et al., 2001) that employs two points but does not rely on
derivatives.

This work was supported by NASA mission Juno and by the Center for Matter
under Extreme Conditions (CMEC) under a grant by the Department of Energy-
National Nuclear Security Administration.

## References

  * Anderson & Schubert (2007) Anderson, J. D., & Schubert, G. 2007, Science, 317, 1384, doi: 10.1126/science.1144835
  * Bodenheimer & Pollack (1986) Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391, doi: 10.1016/0019-1035(86)90122-3
  * Bolton et al. (2017) Bolton, S. J., Adriani, A., Adumitroaie, V., et al. 2017, Science, 356, 821, doi: 10.1126/science.aal2108
  * Brygoo et al. (2021) Brygoo, S., Loubeyre, P., Millot, M., et al. 2021, Nature, 593, doi: 10.1038/s41586-021-03516-0
  * Campbell & Anderson (1989) Campbell, J. K., & Anderson, J. D. 1989, Astron. J., 97, 1485
  * Campbell & Synnott (1985) Campbell, J. K., & Synnott, S. P. 1985, Astron. J., 90, 364
  * Cao & Stevenson (2017) Cao, H., & Stevenson, D. J. 2017, J. Geophys. Res. Planets, 122, 686, doi: 10.1002/2017JE005272
  * Christensen et al. (2020) Christensen, U. R., Wicht, J., & Dietrich, W. 2020, ApJ, 890, 61, doi: 10.3847/1538-4357/ab698c
  * Debras & Chabrier (2019) Debras, F., & Chabrier, G. 2019, The Astrophysical Journal, 872, 100, doi: 10.3847/1538-4357/aaff65
  * Dietrich et al. (2021) Dietrich, W., Wulff, P., Wicht, J., & Christensen, U. R. 2021, Monthly Notices of the Royal Astronomical Society, 505, 3177, doi: 10.1093/mnras/stab1566
  * Durante et al. (2020) Durante, D., Buccino, D. R., Tommei, G., et al. 2020, Geophys. Res. Lett., 47, e2019GL086572
  * Fortney et al. (2016) Fortney, J. J., Helled, R., Nettelmann, N., et al. 2016, in Saturn in the 21st Century, ed. B. Kevin, M. Flasar, N. Krupp, & S. Thomas (Cambridge University Press), 1–28.  https://arxiv.org/abs/1609.06324
  * Galanti & Kaspi (2021) Galanti, E., & Kaspi, Y. 2021, MNRAS, 501, 2352–2362
  * Gonzalez-Cataldo & Militzer (2023) Gonzalez-Cataldo, F., & Militzer, B. 2023, Phys. Rev. Res., in press
  * Gonzalez-Cataldo et al. (2014) Gonzalez-Cataldo, F., Wilson, H. F., & Militzer, B. 2014, Astrophys. J., 787, 79
  * Guillot (2005) Guillot, T. 2005, Ann. Rev. Earth Planet. Sci., 33, 493
  * Guillot et al. (2018) Guillot, T., Miguel, Y., Militzer, B., et al. 2018, Nature, 555, doi: 10.1038/nature25775
  * Gupta et al. (2022) Gupta, P., Atreya, S., Steffes, P. G., et al. 2022, arXiv:2205.12926
  * Helled et al. (2022) Helled, R., Stevenson, D. J., Lunine, J. I., et al. 2022, Icarus, 378, 114937, doi: https://doi.org/10.1016/j.icarus.2022.114937
  * Howard et al. (2023) Howard, S., Guillot, T., Bazot, M., et al. 2023, A&A, 672, A33, doi: 10.1051/0004-6361/202245625
  * Hubbard (2013) Hubbard, W. B. 2013, Astrophys. J., 768, 43, doi: 10.1088/0004-637X/768/1/43
  * Hubbard & Militzer (2016) Hubbard, W. B., & Militzer, B. 2016, Astrophys. J., 820, 80
  * Kaspi et al. (2016) Kaspi, Y., Davighi, J. E., Galanti, E., & Hubbard, W. B. 2016, Icarus, 276, 170
  * Kaspi et al. (2018) Kaspi, Y., Galanti, E., Hubbard, W., et al. 2018, Nature, 555, doi: 10.1038/nature25793
  * Kerley (2004) Kerley, G. I. 2004, Structures of the Planets Jupiter and Saturn
  * Ledoux (1947) Ledoux, P. 1947, Astrophys. J. Lett., 105, 305
  * Lindal et al. (1981) Lindal, G. F., Wood, G. E., Levy, G. S., et al. 1981, Journal of Geophysical Research: Space Physics, 86, 8721, doi: 10.1029/JA086iA10p08721
  * Liu et al. (2019) Liu, S.-F., Hori, Y., Muller, S., et al. 2019, Nature, 572, 355
  * Lodders (2010) Lodders, K. 2010, in Astrophysics and Space Science Proceedings, ed. A. Goswami & B. E. Reddy (Berlin: Springer-Verlag), 379–417
  * Miguel et al. (2022) Miguel, Y., Bazot, M., Guillot, T., et al. 2022, Astron. and Astrophys., 662, A18
  * Militzer (2005) Militzer, B. 2005, J. Low Temp. Phys., 139, 739
  * Militzer (2006) —. 2006, Phys. Rev. Lett., 97, 175501
  * Militzer (2009) —. 2009, Phys. Rev. B, 79, 155105
  * Militzer (2013) —. 2013, Phys. Rev. B, 87, 014202
  * Militzer (2023a) —. 2023a, The Astrophysical Journal, 953, 111, doi: 10.3847/1538-4357/ace1f1
  * Militzer (2023b) —. 2023b, Quadratic Monte Carlo, 03-08-23, Zenodo, doi: 10.5281/zenodo.8038144
  * Militzer (2024) —. 2024, Models for Jupiter’s Interior Structure, Zenodo, doi: 10.5281/zenodo.10471389
  * Militzer et al. (2021) Militzer, B., Gonzalez-Cataldo, F., Zhang, S., Driver, K. P., & Soubiran, F. 2021, Phys. Rev. E, 103, 013203
  * Militzer & Hubbard (2013) Militzer, B., & Hubbard, W. B. 2013, Astrophys. J., 774, 148
  * Militzer & Hubbard (2023) —. 2023, The Planetary Science Journal, 4, 95, doi: 10.3847/PSJ/acd2cd
  * Militzer et al. (2008) Militzer, B., Hubbard, W. B., Vorberger, J., Tamblyn, I., & Bonev, S. A. 2008, Astrophys. J. Lett., 688, 45
  * Militzer et al. (2016) Militzer, B., Soubiran, F., Wahl, S. M., & Hubbard, W. 2016, J. Geophys. Res. Planets, 121, 1552
  * Militzer et al. (2019) Militzer, B., Wahl, S., & Hubbard, W. B. 2019, The Astrophysical Journal, 879, 78, doi: 10.3847/1538-4357/ab23f0
  * Militzer & Wilson (2010) Militzer, B., & Wilson, H. F. 2010, Phys. Rev. Lett., 105, 195701
  * Militzer et al. (2022) Militzer, B., Hubbard, W. B., Wahl, S., et al. 2022, Planet. Sci. J., 3, 185
  * Moll et al. (2017) Moll, R., Garaud, P., Mankovich, C., & Fortney, J. J. 2017, ApJ, 849, 24, doi: 10.3847/1538-4357/aa8d74
  * Moore et al. (2018) Moore, K. M., Yadav, R. K., Kulowski, L., et al. 2018, Nature, 561, 76, doi: 10.1038/s41586-018-0468-5
  * Moore et al. (2022) Moore, K. M., Barik, A., Stanley, S., et al. 2022, Journal of Geophysical Research: Planets, 127, e2022JE007479, doi: https://doi.org/10.1029/2022JE007479
  * Morales et al. (2013) Morales, M. A., McMahon, J. M., Pierleonie, C., & Ceperley, D. M. 2013, Phys. Rev. Lett., 110, 065702
  * Morales et al. (2010) Morales, M. A., Pierleoni, C., Schwegler, E., & Ceperley, D. M. 2010, Proc. Nat. Acad. Sci., 107, 12799
  * Müller et al. (2020) Müller, S., Helled, R., & Cumming, A. 2020, Astronomy and Astrophysics, 638, id.A121
  * Nettelmann et al. (2012) Nettelmann, N., Becker, A., Holst, B., & Redmer, R. 2012, Astrophys. J., 750, 52
  * Nettelmann et al. (2013) Nettelmann, N., Püstow, R., & Redmer, R. 2013, Icarus, 225, 548, doi: 10.1016/j.icarus.2013.04.018
  * Nettelmann et al. (2021) Nettelmann, N., Movshovitz, N., Ni, D., et al. 2021, Planetary Science Journal, 2, 241, doi: 10.3847/psj/ac390a
  * Ni (2018) Ni, D. 2018, A&A, 613, A32, doi: 10.1051/0004-6361/201732183
  * Niemann et al. (1998) Niemann, H. B., Atreya, S. K., Carignan, G. R., et al. 1998, J. Geophys. Res., 103, 22831
  * Press et al. (2001) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2001, Numerical Recipes in C++ (Cambridge, UK: Cambridge University Press)
  * Saumon & Guillot (2004) Saumon, D., & Guillot, T. 2004, Astrop. J, 609, 1170
  * Seager et al. (2007) Seager, S., Kuchner, M. J., Hier-Majumder, C. A., & Militzer, B. 2007, The Astrophysical Journal, 669, 1279, doi: 10.1086/521346
  * Soubiran & Militzer (2016) Soubiran, F., & Militzer, B. 2016, Astrophys. J., in press, http://arxiv.org/abs/1606.04162
  * Spilker (2019) Spilker, T. 2019, Science, 364, 1046
  * Stevenson (1982) Stevenson, D. J. 1982, Annual Review of Earth and Planetary Sciences, 10, 257, doi: 10.1146/annurev.ea.10.050182.001353
  * Tollefson et al. (2017) Tollefson, J., Wong, M. H., de Pater, I., et al. 2017, Icarus, 296, 163, doi: https://doi.org/10.1016/j.icarus.2017.06.007
  * von Zahn et al. (1998) von Zahn, U., Hunten, D. M., & Lehmacher, G. 1998, J. Geophys. Res., 103, 22815
  * Vorberger et al. (2012) Vorberger, J., Donko, Z., Tkachenko, I. M., & Gericke, D. O. 2012, Phys. Rev. Lett., 109, 225001
  * Vorberger et al. (2010) Vorberger, J., Gericke, D. O., Bornath, T., & Schlanges, M. 2010, Phys. Rev. E, 81, 046404
  * Vorberger et al. (2007) Vorberger, J., Tamblyn, I., Militzer, B., & Bonev, S. 2007, Phys. Rev. B, 75, 024206
  * Wahl et al. (2021) Wahl, S. M., Thorngren, D., Lu, T., & Militzer, B. 2021, Astrophys. J., 921, 105
  * Wahl et al. (2013) Wahl, S. M., Wilson, H. F., & Militzer, B. 2013, Astrophys. J., 773, 95
  * Wahl et al. (2017) Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649, doi: 10.1002/2017GL073160
  * Wilson & Militzer (2010) Wilson, H. F., & Militzer, B. 2010, Phys. Rev. Lett., 104, 121101
  * Wilson & Militzer (2012a) —. 2012a, Astrophys. J., 745, 54
  * Wilson & Militzer (2012b) —. 2012b, Phys. Rev. Lett., 108, 111101
  * Wilson & Militzer (2014) —. 2014, Astrophys. J., 973, 34
  * Wong et al. (2004) Wong, M., Mahaffy, P. R., Atreya, S. K., Niemann, H. B., & Owen, T. C. 2004, Icarus, 171, 153
  * Zharkov & Trubitsyn (1978) Zharkov, V. N., & Trubitsyn, V. P. 1978, Physics of Planetary Interiors (Pachart), 388