Draft version October 18, 2024

Typeset using LATEX twocolumn style in AASTeX631

Jupiter Evolutionary Models Incorporating Stably Stratified Regions

Roberto Tejada Arevalo ,1 Ankan Sur ,1 Yubo Su ,1 and Adam Burrows 1, 2

1 Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08540, USA
2 Institute for Advanced Study, 1 Einstein Dr, Princeton, NJ 08540, USA

ABSTRACT
arXiv:2410.12899v1 [astro-ph.EP] 16 Oct 2024

We address the issue of which broad set of initial conditions for the planet Jupiter best matches
the current presence of a “fuzzy core” of heavy elements, while at the same time comporting with
measured parameters such as its effective temperature, atmospheric helium abundance, radius, and
atmospheric metallicity. Our focus is on the class of fuzzy cores that can survive convective mixing
to the present day and on the unique challenges of an inhomogeneous Jupiter with stably-stratified
regions now demanded by the Juno gravity data. Hence, using the new code APPLE, we attempt to put
a non-adiabatic Jupiter into an evolutionary context. This requires not only a mass density model,
the major relevant byproduct of the Juno data, but a thermal model that is subject to interior heat
transport, a realistic atmospheric flux boundary, a helium rain algorithm, and the latest equation of
state. The result is a good fit to most major thermal, compositional, and structural constraints that
still preserve a fuzzy core and that should inform future more detailed models of the current Jupiter
in the context of its evolution from birth.

Keywords: Solar system gas giant planets, Extrasolar gaseous giant planets, Planetary science, Plane-
tary structure

1. INTRODUCTION This ferment suggests that it is timely to attempt a

The Juno (Bolton et al. 2017) and Cassini (Matson deeper understanding of the structure and evolution of
et al. 2003) missions have significantly transformed our Solar System gas giant planets, in part to advance our
knowledge of the current interiors of Jupiter and Sat- knowledge of gas giant exoplanet formation and evolu-
urn, revealing the likely presence of convectively stable tion.
regions and extended heavy-element cores (Fuller et al. Recent discoveries of extended heavy element (Z) re-
2014; Fuller 2014; Wahl et al. 2017; Folkner et al. 2017; gions within Jupiter and Saturn may offer insight into
Nettelmann 2017; Iess et al. 2019; Militzer et al. 2019; their formation and evolution. The top panel of Figure 1
Debras & Chabrier 2019; Durante et al. 2020; Mankovich illustrates four examples of initial composition profiles
& Fuller 2021; Nettelmann et al. 2021; Militzer et al. from formation models of Jupiter, as presented by Vazan
2022; Militzer & Hubbard 2024). In addition, over the et al. (2016), Lozovsky et al. (2017), Helled & Steven-
past decade exoplanet research is starting to reveal the son (2017), and Stevenson et al. (2022). The bottom
masses, sizes, and demographic characteristics of thou- panel displays the current inferred heavy-element distri-
sands of exoplanets (e.g., see Winn & Fabrycky 2015; bution of Jupiter, as derived from Juno gravity data by
Fulton et al. 2017; Petigura et al. 2018; Rosenthal et al. Militzer & Hubbard (2024), Militzer et al. (2022), Net-
2021). Furthermore, missions such as JWST (Gard- telmann et al. (2021), Debras & Chabrier (2019), and
ner et al. 2006) and future ventures, including PLATO Wahl et al. (2017). These present-day heavy-element
(Rauer et al. 2014), and ARIEL (Tinetti et al. 2021), distributions of Jupiter encompass a wide range of pos-
are poised to launch a new era of exoplanet exploration. sibilities that agree with the Juno data, though only
a limited number of evolutionary models have incorpo-
rated formation-based initial conditions (Müller et al.
Corresponding author: Roberto Tejada Arevalo 2020; Stevenson et al. 2022). Evolutionary models in-
[email protected] formed by formation models have faced challenges in re-
producing the heavy-element profiles observed by Juno
2

(Müller et al. 2020). Particularly, the initial thermal In §2, we describe our evolutionary code, APPLE (Sur
profiles from formation models frequently lead to homo- et al. 2024), our updated boundary conditions (Chen
geneous composition profiles. The initial high interior et al. 2023), the EOSes used for the envelope and com-
entropies in the formation profiles used by Müller et al. pact core, and the H-He miscibility curves used for
(2020) led to extensive convective mixing throughout all models (Tejada Arevalo et al. 2024). Throughout
much of the extended core region over Gyr timescales. this work, we use the term “compact core” to refer
Thus far, evolutionary studies incorporating Jovian to the inner rocky, “Z = 1” region and use “fuzzy,”
extended cores have left open the questions of 1) which “extended,” or “dilute,” core interchangeably to indi-
broad set of initial conditions are best suited to match cate the region of the envelope with enhanced Z and/or
the current presence of a fuzzy core and measured pa- where dZ/dm < 0, where m is the interior mass. In §3,
rameters and 2) the character of helium phase separation we provide a set of models with varying initial condi-
expected to occur in Jupiter (Stevenson 1975; Wilson tions. We vary the total entropy profile and the sizes
& Militzer 2010; Hubbard & Militzer 2016; Mankovich of the extended cores and determine which initial con-
et al. 2016; Mankovich & Fortney 2020; Nettelmann ditions best preserve a fuzzy core, while approximately
et al. 2024). To begin to address the question of the pos- matching Jupiter’s current effective temperature and ra-
sible preservation of a fuzzy core from birth, Knierim & dius. Section §4 presents one favored Jupiter model that
Helled (2024) showed that cooler interior entropy (S) closely matches the observables in Table 1 and describes
profiles might maintain an extended initial Z profile the parameters necessary to obtain this model. We pro-
over evolutionary timescales. In contrast, hotter interior vide concluding remarks in §5.
S profiles lead to convective mixing and homogeniza-
tion over shorter timescales due to the expansion of the Table 1. Present-day Jupiter Observables
outer convective zone. Those authors concluded that the Considered
primordial entropy profile is the dominant factor shap-
ing the evolution of the composition profile and broader Parameter Value
planetary structure. Their conclusions inform formation Teff [Kelvin] 125.57 ± 0.07
models and evolutionary calculations but leave open the Yatm 0.234 ± 0.005
question of how Jupiter’s current effective temperature Z⊙ (out) 1.5–5 (∼3)
and atmospheric helium mass fraction can be explained. Equatorial Radius [km] 71, 492 ± 4
In this paper, we address these two aspects simul-
Note—The effective temperature comes
taneously. We conduct evolutionary models to deter- from Li et al. (2012), the atmospheric
mine general initial distributions of heavy elements that helium abundance from Von Zahn et al.
can best match Jupiter’s current effective temperature (1998), the outer metallicity from various
(Teff = 125.57 K; Li et al. 2012), the atmospheric he- measurements of the atmosphere (Li et al.
lium abundance as measured by the Galileo entry probe 2020), and the radius measurements come
(Yatm = 0.234 ± 0.005; Von Zahn et al. 1998), the outer from Seidelmann et al. (2007). The metal
abundances can range between 1.5 and 5
metallicity abundance of roughly three solar (∼3Z⊙ ) as
solar but have a concentration of roughly
measured by Voyager, Galileo, and Juno (see Helled 3 times solar (Figure 2 of Helled & Howard
& Howard 2024; Guillot et al. 2022; Li et al. 2020), 2024)
and the current equatorial radius of Jupiter (71,492
km; Seidelmann et al. 2007), while still preserving an
extended heavy-element core. These values are sum-
marized in Table 1. This work presents the first evo-
lutionary models of Jupiter with extended cores with
helium rain, updated atmospheric boundary conditions 2. METHODS
that account for ammonia clouds, irradiation, and he-
The evolutionary models presented throughout are
lium abundance dependence (Chen et al. 2023), the
calculated using APPLE (Sur et al. 2024) and include the
most current hydrogen-helium (H-He) and Z equations
effects of solid-body rotation via conservation of angu-
of state (EOS) (Chabrier & Debras 2021; Haldemann
lar momentum (see Equation 1 of Sur et al. 2024). A
et al. 2020), and the best available pressure, tempera-
final spin period of 9:55:29.704 hr is assumed (Yu &
ture, and Y -dependent H-He miscibility curves (Loren-
Russell 2009) and we neglect differential rotation. All
zen et al. 2009, 2011; Schöttler & Redmer 2018).
models calculate and compare to the equatorial radius
of Jupiter. Moreover, all models include the irradiated
 3

atmospheric boundary conditions of Chen et al. (2023).

These boundary conditions are improvements upon the
Fortney et al. (2011) boundary conditions in that they
include solar insolation in the detailed atmosphere cal-
culation, the effects of ammonia clouds, and naturally
obtain the updated albedo of Jupiter of 0.5 (Li et al.
2018). We use the H-He EOS of Chabrier & Debras
(2021, CD21) and mix it via the volume addition law
(see Tejada Arevalo et al. 2024; Chabrier et al. 2019;
Saumon et al. 1995, for details) with the latest multi-
phase water EOS tables from Haldemann et al. (2020,
AQUA), which includes the supercritical and superionic
water EOS of Mazevet et al. (2019).
The most recent H-He miscibility work is that of
Schöttler & Redmer (2018, herein SR18). Both SR18
and the older (and hotter) miscibility curves of Loren-
zen et al. (2009, 2011, herein LHR0911) are thus far the
most sophisticated Y -dependent ab initio H-He misci-
bility curves publicly available. We use both of these
miscibility curves in this work and adjust, as has been
done in past evolutionary work (Nettelmann et al. 2015;
Püstow et al. 2016; Mankovich et al. 2016; Mankovich
& Fortney 2020; Howard et al. 2024), the temperatures
of their miscibility curves at a given pressure and Y
fraction when trying to obtain the current helium abun-
dance of Jupiter’s atmosphere. In the absence of any de-
tailed knowledge of helium droplet formation (Stevenson
& Salpeter 1977; Mankovich et al. 2016) in regions of ex-
pected H-He phase separation (∼1–3 Mbar), APPLE uses
advection-diffusion methods to transport helium in the
miscibility region toward the inner convective regions
(see equation 49 in Sur et al. 2024), while conserving
total helium mass to machine precision. In all our evo-
lutionary calculations, we identify the corresponding he-
Figure 1. Top: Examples of the distribution of Z for forma- lium mass fraction at immiscibility (Ymisc (P, T )) at the
tion models from Stevenson et al. (2022), Vazan et al. (2016),
local temperatures and pressures of each Jupiter profile.
Lozovsky et al. (2017), and Helled & Stevenson (2017). The
profile names are also given since these authors explored To determine an immiscible region, we compare the ad-
multiple models. These formation distributions are taken justed helium mass fraction (Y ′ = Y /(Y + X), where
at the time when the proto-Jupiters reached final mass ac- X is the hydrogen mass fraction) with Ymisc (P, T ). If
cording to the authors. Bottom: Current composition gra- Ymisc (P, T ) < Y ′ , then APPLE deploys a modified advec-
dients as inferred by the Juno gravity data from Militzer & tion term in the diffusion equation to model helium rain
Hubbard (2024, MH24), Militzer et al. (2022, MH22), Net- (see section 5.4.2 of Sur et al. 2024).
telmann et al. (2021), Debras & Chabrier (2019, DC19), and
The regions unstable to convection are determined via
Wahl et al. (2017), as adapted by Helled et al. (2022). The
heavy element distributions inferred by various authors us-
the Ledoux criterion,
ing the Juno data employ different methods, equations of ∂S X ∂S
 
∂Xi
state, and assumptions concerning the current thermal pro- − < 0, (1)
file of Jupiter. These differences can potentially account for ∂r i
∂X i ρ,P ∂r
the differences seen in the inferred Z profiles. Hence, even
though Juno provided precise values of the J2n gravity mo- where the Xi are the helium (Y ) and heavy-element
ments, there remains a large range of potential Z distribu- (Z) mass fractions. A full derivation and thermody-
tions that might explain the Juno data. namic discussion of this criterion can be found in Te-
jada Arevalo et al. (2024). In convective regions, a
mixing-length flux procedure for both energy and com-
4

position, as described in Sur et al. (2024), is employed. are put in a compact core made up of 50% water and
In stably-stratified regions, the conductive thermal and 33.33% post-perovskite (MgSiO3 ), with the remainder
compositional fluxes are calculated using the French iron using the EOS used in Zhang & Rogers (2022),
et al. (2012) diffusivities, as modified by Sur et al. which is Keane’s EOS (Keane 1954; Stacey & Davis
(2024). Each timestep is determined to ensure the tem- 2004) expanded to high pressures and temperatures for
perature, helium fraction, and heavy-element fraction do this work. For simplicity, water is the only heavy species
not change in a zone by more than 1% and 750 spatial in the envelope. Furthermore, all initial models assume
zones are employed. an initial constant helium abundance relative to hydro-
gen at the protosolar value of Y ′ = Y /(X + Y ) = 0.277,
3. PRESERVATION OF INITIAL as calculated by Bahcall et al. (2006).
HEAVY-ELEMENT EXTENDED CORES 3.1. The Survival and Consequences of Flat-Top Fuzzy
This section explores the general contexts in which an Cores
extended “fuzzy” Jovian core of heavy elements might The initial heavy metal (Z) and entropy (S) profiles
survive until the present day. As illustrated in the bot- presented in this section follow the same general trends
tom panel of Figure 1, the range of fuzzy core struc- recently found by Knierim & Helled (2024). Namely,
tures inferred from the Juno data is quite broad. De- that initial Z profiles located at larger radii convectively
bras & Chabrier (2019), Militzer et al. (2022, MH22), mix1 more rapidly than those located deeper in the struc-
and Militzer & Hubbard (2024, MH24), infer flat Z pro- ture. From left to right, Figure 2 depicts what happens
files out to a steep Z gradient near ∼0.6 RJ , where RJ when the initial location of the steep gradient is moved
is the radius of Jupiter. If composition profile gradi- progressively outward. The convective zones move in-
ents are overrun by convection, they will then be mixed ward, and these zones convectively mix the extended
and will not survive over Gyr timescales. A narrow Z cores, as shown in the second row of Figure 2. Models
gradient, flanked by flat interior and exterior profiles, that mix throughout access all the heat of the envelope,
presents a physical challenge, since such a configuration leading to higher effective temperatures. In Figure 3,
would naturally arise only under conditions of homog- we show that the models with larger extended cores ex-
enizing convection. A fully convective, adiabatic initial perience lower effective temperatures due to inefficient
envelope, which would homogenize the planetary struc- heat transport from their interiors at early times. The
ture, is inconsistent with the presence of a fuzzy core. model with its initial fuzzy core extending out to 0.4
Therefore, we initialize the envelope as stably stratified, MJ (shown in blue) can hold its extended core over Gyr
imposing an initial positive entropy gradient to achieve timescales and shows higher effective temperatures at
this. Hence, we investigate the general conditions under early times since it can access more of this interior heat.
which convection can be suppressed, allowing a fuzzy The evolution of the equatorial radius is shown in the
core and its associated stabilizing Z gradients to per- bottom panel of Figure 3. These example models illus-
sist. trate the dependence of the effective temperature on the
In §3.1 we explore how a class of initial steep Z profile initial Z and S profiles. The 0.4 MJ model in Figure 3
gradients in the middle of the planet surrounded by flat reaches the current Jovian effective temperature at ∼3
Z distributions might survive Gyr timescales of ther- Gyrs, and its extended core is preserved only by narrow
mal, structural, and compositional evolution. In §3.2, stable regions enclosed by convective regions with no
we explore the dependence of the Z profile gradient on extended stable zone. The convective regions can pene-
internal entropy profiles by varying from cold to hot in- trate the Z profile gradient when the entropy profile of
ternal entropies. We present these models graphically in the outer region cools to lower values close to those of
Figures 2–5. This is followed in §3.3 by an investigation the innermost stable region of the model (see “staircase”
into how fuzzy cores with a different class of initial Z discussion in Vazan et al. 2018). These models are not
profiles which gradually decrease in Z abundance from representative of Jupiter’s evolution, but for demonstra-
the inside out might fare over time. This different class tion purposes we have used the unmodified LHR0911
of initial conditions resembles the more smoothly con- miscibility curve, shown in the bottom row of Figure 2.
tinuous present-day Jovian Z distributions inferred by
Nettelmann et al. (2021). We emphasize that we are 3.2. The Consequences of Initial Entropy
not directly addressing in this study planet formation
scenarios that might produce such profiles.
All the models explored in §3.1–§3.3 are given a total 1 Our working definition of “mixing” are regions where ∂Z/∂m =
heavy-element mass of 27 M⊕ . Of this total, 2.5 M⊕ 0.
 5

Initial fuzzy cores with higher (or “hotter”) S pro- lieve that this problem might be alleviated with higher
files mix more rapidly due to convection than those initial internal temperatures near the steep Z gradient,
with lower (or “colder”) S profiles (Müller et al. 2020; but the associated elevated internal entropies will mix
Knierim & Helled 2024). From left to right, Figure 4 the profile gradients, as demonstrated in Figure 4. A
illustrates this dependence of the Z profile on the S pro- negative Z gradient can be introduced interior to the
file. The location of the Z profile is maintained in Fig- steep portion of the profile to keep the extended core
ure 4 at 0.6 MJ to illustrate the S profile effects on the Z from fully mixing. This maintains the extended core out
profile. To guide the eye, a characteristic entropy at the to 0.5 MJ , but the effective temperatures remain below
position of the Z gradient (ranging from 6.5 to 8.5 kB 120 K at 4.56 Gyrs. The requirement to maintain the
baryon−1 ) is marked by black dash-dot lines in the first steep Z gradient at M ∼0.5–0.6 MJ conflicts with the
row. The inward advancement of the convective front is need to sustain high enough temperatures in the convec-
inhibited by the Z gradient in colder S profiles, causing tive region to achieve an adequate Teff for this kind of
the outer convective regions to cool more rapidly with- profile. Thus, while a steep Z profile gradient can persist
out heat replenishment from the interior. This outer over Gyr timescales under certain conditions, we disfa-
adiabatic cooling results in lower effective temperatures vor such initial configurations due to the requirement to
and smaller radii, as shown in Figure 5. As outer tem- match Jupiter’s effective temperature.
peratures drop, they intersect the shifted SR18 miscibil- Facing these challenges, we explored a more gradu-
ity curve, as shown in the left-most column. As such, ally sloping class of initial conditions. The exercise in
preserving these extended cores results in effective tem- §3.2 is repeated in Figures 6 and 7, where the initial
peratures much colder than the present-day Jovian effec- fuzzy cores are instead gradual throughout the inner en-
tive temperature. As in the previous section, and since velope. The unmodified miscibility curves of LHR0911
these are only example models, we arbitrarily increased have been used to illustrate helium rain in these exam-
the demixing temperatures of the SR18 miscibility curve ple models. As in Figure 4, Figure 6 depicts models for
by +550 K for these example calculations. which a characteristic entropy (shown in the first row)
ranges from 6.5 to 8.5 kB baryon−1 , demonstrating the
3.3. The Survival and Consequences of Gradually effect of the gradual increase of the S profile from left
Sloping Fuzzy Cores to right. The rightmost model described by the highest
entropy mixes completely by 0.4 Gyr, while the others
The shortcomings and challenges of the profiles ex-
maintain their extended fuzzy cores to varying degrees.
plored in §3.1 and §3.2 motivate the exploration of a
In Figure 7, we show that the gradual core evolutionary
different kind of initial distribution. The reader should
model characterized by 7.5 kB baryon−1 (third column),
keep in mind that although formation models have pro-
which indeed preserves a stable fuzzy core, can reach
vided initial conditions to be used in evolutionary mod-
an effective temperature closer to that of Jupiter. This
els, these have not in the past been shown to preserve
class of initial conditions allows for a more gradual evo-
fuzzy cores larger than ≳ 0.2 MJ (Müller et al. 2020).
lution of the effective temperature (see Figure 7). As
Based on our example evolution models described in §3.1
shown in Figure 6, internal entropies of between 7 and
and §3.2, initially steep Z profile gradients extending be-
7.5 kB baryon−1 appear adequate for this kind of initial
yond 0.4 MJ , as advocated in the models of MH24, do
condition.
not seem at the same time to easily explain the Jovian
effective temperature in an evolutionary context. Those
initial Z distributions which have an initially stable in- 4. EVOLUTION OF JUPITER
terface at interior masses greater than 0.5 MJ mix too To find a favored model with the class of initial condi-
rapidly to provide a stable Z region at the present age, tions discussed in §3.3, we conducted a parameter search
as shown in Figure 2. If such a profile is initialized with by varying the total mass of heavy elements, Mz , from
colder interior entropy profiles, then the stable Z region 25 M⊕ to 50 M⊕ . For each Mz , the compact core mass
will endure (Figure 4). The outer regions, though, will (Mc ) was varied from 2 to 15 M⊕ . For each of these mod-
cool rapidly and fail to account for Jupiter’s current ef- els, we varied the miscibility temperature modifications
fective temperature, as illustrated in Figure 5. These (∆T ) from -500 K to 500 K (for LHR0911) and from 0
colder models can be provided with higher exterior en- K to 2,500 K (for SR18). Figures 2–7 show that internal
tropy profiles to achieve higher initial effective temper- characteristic entropy values of 7.5 kB baryon−1 strike a
atures. Although such models experience higher initial balance between preserving an extended core and yield-
effective temperatures, they nonetheless cool to effective ing appropriate effective temperatures for Jupiter at the
temperatures below 120 K at later ages. One may be- present age. Thus, each model was initialized with a
6

characteristic interior entropy of 7.5 kB baryon−1 and an by flat Z profiles left over from formation will lead to
outside entropy of 8.2 kB baryon−1 . Although a wide va- effective temperatures that are too low for the present
riety of initial profile distributions could be constructed Jupiter. This is due to the lack of heat resupply from
and brought close to matching observed data, this par- the interior, since the stable regions located at 0.4 MJ
ticular class of initial Z distributions is naturally suited ≤ M ≤ 0.6 MJ inhibit heat transport.
to match the observables presented in Table 1. Specifi- Since the aim of this study was to investigate which
cally, it preserves an inner stable region better and pro- initial distributions 1) can retain extended stable cores
vides enough erosion of the extended core via convective and 2) can account for Jupiter’s effective temperature,
mixing to obtain a metallicity of 3Z⊙ in the outer re- atmospheric helium abundance, and equatorial radius,
gions. we chose initial distributions better suited for that pur-
Our preferred model for Jupiter’s evolution is pre- pose. Past evolutionary models informed by formation
sented in Figures 8 and 9. The observed parameters models have thus far been unable to reproduce the in-
of this model are compared to the data for Jupiter in ferred Juno Z profiles (Müller et al. 2020; Stevenson
Table 2. The relative errors of this model’s effective et al. 2022). The evolutionary models of Vazan et al.
temperature and radius are within 0.05% and 0.2% re- (2018) end with an extended core region out to 0.38
spectively, all while Yatm reaches the value of 0.234. MJ but did not attempt to impose the constraint of the
In this model, helium rain begins at 4 Gyrs. Impor- measured atmospheric helium abundance, used the older
tantly, the miscibility temperatures required to achieve SCvH EOS (Saumon et al. 1995), and did not employ
Yatm = 0.234 are +300 K higher than the original an up-to-date atmospheric boundary. The stable region
LHR0911 immiscibility temperatures. At the same time, predicted by the formation and evolutionary models of
a significant temperature modification of +2,200 K to Müller et al. (2020) extends to only 0.2 MJ . More re-
the SR18 miscibility curve is required. The required cently, Howard et al. (2023a) investigated the H-He EOS
miscibility curves of LHR0911 and SR18 are compared dependence of static Jupiter models, also favoring an
to the original model in Figure 10, with experimental extended core out to 0.2 MJ . The evolutionary models
H-He phase separation data from Brygoo et al. (2021) presented in our work begin with a completely stable
interpolated and plotted for comparison. The elevated structure but also end with an extended core region out
immiscibility temperatures are in part a result of the to ∼0.3 MJ .
higher profile temperatures in the convective regions Our work is distinguished from the works of Vazan
due to the presence of heavy metals in the envelope. et al. (2018), Müller et al. (2020), Mankovich & Fortney
The difference seen in Figure 10 between the shifted (2020), Howard et al. (2024) in the following ways:
LHR0911 and SR18 miscibility curves is a consequence
of their different shapes in the ∼1−4 Mbar pressure re- 1. The use of updated boundary conditions which ac-
gion of relevance to helium rain in Jupiter. A recent count for irradiation, the Bond albedo, and am-
evolutionary calculation that incorporates only H-He in monia clouds (Chen et al. 2023). These boundary
the interior of Jupiter’s envelope predicts colder demix- conditions are Y-dependent, adjusting to the outer
ing temperatures, by as much as -1,250 K concerning Y abundance during helium rain.
the original LHR0911 miscibility curve (Howard et al.
2. The use of an updated H-He EOS (Chabrier &
2024). The best-fit Jupiter model of Mankovich & Fort-
Debras 2021) and the incorporation of a multi-
ney (2020) seems to have a homogenous envelope metal-
phase water EOS (Haldemann et al. 2020).
licity of only ∼one solar, with the correspondingly lower
envelope temperatures and lower required shift to the 3. An exploration of different classes of initial condi-
SR18 miscibility curve. We conclude that the inclusion tion profiles for S and Z.
of heavy metals in the envelope results in higher temper-
atures and, therefore, warrants higher H-He demixing 4. The application of the two most recent H-He mis-
temperatures. cibility curves (LHR0911, SR18) used in the evo-
Our preferred Jupiter model is the first to incorpo- lution of Y in the presence of Z profile gradients.
rate the evolution of initially stable extended cores to These H-He miscibility curves are allowed to adapt
match the effective temperature, radius, outer metal and adjust to the local pressures, temperatures,
abundance, and atmospheric helium abundance with he- and helium abundance.
lium rain within reasonable uncertainties. One of the
5. The deployment of a new evolutionary code
subsidiary conclusions we draw from the explorations in
(APPLE ; Sur et al. 2024) with updated micro-
§3.1 and §3.2 is that steep Z profile gradients surrounded
physics, such as conductivities, helium rain mod-
 7

eling methods, and the inclusion of centrifugal ef- than to constrain the possible character of a fuzzy core
fects on the structures. at formation, if that measured extended core is primor-
dial.
5. CONCLUSIONS We have left to future work attempts to fit in the
context of evolutionary models the current Juno lower-
Generally, the Z profiles that retain stable regions are
moment (J2 and J4 ) gravity data. This is in part due
deeper in the planet (≤ 0.4 MJ ) compared to those in-
to the difficulty researchers have had fitting them in
ferred by DC19, MH22, and MH24 and must begin with
the context of Jupiter’s measured surface metallicity
a characteristic interior entropy of ∼7.5 kB baryon−1 .
(Nettelmann 2017; Nettelmann et al. 2021; Helled et al.
Lower interior entropies lead to colder outer tempera-
2022). Such static fits to Jupiter (see Figure 1) seem of-
tures, and higher interior entropies lead to more con-
ten to require low metallicities, and sometimes negative
vective mixing (see Figure 4). In addition, metallicity
metallicities (Nettelmann et al. 2021) in the outer enve-
profile gradients surrounded by flat profiles located at
lope. Debras & Chabrier (2019) has proposed a positive
≳ 0.4 MJ will either convectively mix or lead to effective
Z gradient in the outer regions of Jupiter (red line in the
temperatures significantly colder than 125 K. Including
bottom panel of Figure 1) while Howard et al. (2023b)
a stabilizing region interior to the steep profile gradient
find that a positive Z gradient of such magnitude re-
leads to a more stable core but aggravates the fit with
quires an excessive amount of late accretion inconsistent
effective temperature. As a result, initial Z distribu-
with collisional evolution models of the Solar System.
tions with steep stable profile gradients located in the
As proposed by Stevenson (1982), core erosion could
outer regions (M ≥ 0.6 MJ ) may not be the result of
be responsible for generating Z gradients and fuzzy cores
formation processes.
on evolutionary timescales, and core erosion will be a
Obtaining ∼3Z⊙ metallicity in the outer envelope
subject of future work. We have not in this paper in-
leads to higher temperatures in the outer convective
corporated the potential effects of semi-convection and
region since (generally) higher metallicity temperature
semi-convective interfaces (Mirouh et al. 2012; Wood
profiles are significantly hotter than pure H-He profiles.
et al. 2013; Moll et al. 2016) that have been shown
As a result, to deplete sufficient helium from the atmo-
to endure over evolutionary timescales in the presence
sphere, significantly higher H-He demixing temperatures
of rotation (Fuentes et al. 2023, 2024), though Müller
are required. In this work, the LHR0911 and SR18 mis-
et al. (2020) found they had no significant effect with-
cibility curves require higher temperatures by ∼300 K
out rotation. Clearly, the effects of both rotation and
and ∼2200 K, respectively.
semi-convection require further scrutiny.
While the observed metallicity of Jupiter’s atmosphere
Uncertainties exist in the hydrogen-helium and heavy-
may be due to recent impact events (see discussions in
element equations of state, the H-He demixing tempera-
Debras & Chabrier 2019; Helled et al. 2022), it remains
tures and pressures, and in other material properties at
unlikely that positive Z profile gradients will remain sta-
high pressures that render both evolutionary and static
ble over even megayear timescales. As such, we have as-
models provisional. In the context of the high Zs in-
sumed here that the atmospheric metallicity of Jupiter
ferred in current models, the lack of physically sophis-
is the same throughout the outer convective region.
ticated equations of state for hydrogen/helium/heavy-
We emphasize that the fits to the Juno gravity data
element mixtures for the range of compositional mix-
constrain only the mass density profile in the current
tures discussed in the literature is a particularly impor-
Jupiter. The interior thermal distributions are not con-
tant shortcoming of all current models of Jupiter (and
strained at all and are therefore assumed by all re-
Saturn).
searchers attempting to use the Juno gravity data. This
Despite all these uncertainties, evolutionary models
accounts for the wide range of inferred density distribu-
are crucial for determining the long-term stability of
tions and forcefully suggests that the interior thermal
these fuzzy cores in the presence of convective mixing,
profiles are very much unconstrained. This is where evo-
in informing formation models, and in clarify what evo-
lutionary models can play a useful role. Since they fun-
lutionary paths are feasible, given the various observa-
damentally depend upon the thermal interior and the
tional constraints. Their continued development and re-
associated energy transport model, they bring to bear
finement are key to the eventual understanding of this
other constraints, such as the atmospheric effective tem-
class of planet, whose importance extends beyond the
perature at the current epoch. They also can connect
Solar System giants to the giant exoplanets now being
the formation phase with the current phase, something
discovered in profusion and that await a definitive un-
Nature does and theorists should attempt. However, our
derstanding of our local variants.
current models shed little light on Jovian birth, other
8

Funding for this research was provided by the Center

for Matter at Atomic Pressures (CMAP), a National
Science Foundation (NSF) Physics Frontier Center, un-
der Award PHY-2020249. Any opinions, findings, con-
clusions, or recommendations expressed in this material
are those of the author(s) and do not necessarily reflect
those of the National Science Foundation. RTA is sup-
ported by the Ford Foundation Predoctoral Fellowship.
YS is supported by a Lyman Spitzer, Jr. Postdoctoral
Fellowship at Princeton University. We thank Dr. Matt
Coleman and Princeton Research Computing for their
computational and technical resource support. We fur-
ther thank Yi-Xian Chen, Dr. Jisheng Zhang, and Dr.
Akash Gupta for their useful discussions.
 9

Figure 2. Example evolutionary calculations of four mock models with varying initial sizes of the extended core. From left to
right, the models increase in extent from 0.4 to 0.8 MJ . All models have a total envelope heavy-element mass of 27 M⊕ with
2.5 M⊕ in the compact core (Z = 1). The four models in this study assume the same initial entropy of 7.5 kB baryon−1 at the
outer edge of the fuzzy core but shift in location to match the location of each stable region. The surrounding flat Z profile
regions have been initialized with stabilizing S profile gradients. None of these models should be taken as the preferred models
of Jupiter. Rather, these models are meant to explore the contexts for the survival of the given fuzzy cores to late times. Since
these are example models only, we have used the original LHR0911 miscibility curves to demonstrate helium rain. The first
row depicts the entropy profile and the inward extension of the convective zone as each structure cools. At later ages, when
the outer region has cooled more than the inner stable region, the convective zones surround the steep stable regions, leaving
behind so-called “staircases” (see also Vazan et al. 2018). The second row shows the evolution of the Z profile. Generally, these
figures show that shallower (i.e., located closer to the surface) extended cores experience more convective mixing than deeper
extended cores. The third row shows the evolution of the helium abundance profile. Due to the conservation of mass and the
setting of the hydrogen/helium mass ratio equal to the proto-solar value (Y /(X + Y ) = 0.277; Bahcall et al. 2006), the helium
abundance increases with decreasing heavy metal abundance. The last row shows the evolution of the temperature profile with
the helium rain region of the 0.4 MJ model shaded in red. Note that all but the 0.4 MJ model mixes entirely by 4.56 Gyrs. As
a result, the temperatures in the convective envelope of the 0.4 MJ model are colder than the rest at this age. At early ages,
however, the temperature profiles of this model are higher than the rest, since it can access more internal heat due to a deeper
stable region. Deeper extended cores allow more heat to be available to the outer convective zone and this maintains higher
outer envelope temperatures.
10

Table 2. Jupiter Data vs. Derived Evolutionary Model

Parameter Value Derived Model Values

Teff [K] 125.57 ± 0.07 125.44
Yatm 0.234 ± 0.005 0.234
Zout [Z⊙ ] ∼3 3.26
Equatorial Radius [km] 71, 492 71, 697
Note—Jupiter observable data compared to the values of our pre-
ferred model from §4 at 4.56 Gyrs. This particular model has a
total heavy element mass of Mz = 44 M⊕ and a compact core
mass of Mc = 14 M⊕ and uses the LHR0911 miscibility curve
with a demixing temperature modification of +350 K to reach the
observed helium mass fraction of 0.234.

Table 3. H-He Miscibility Temperature Shifts

Authors Misc. Curve (∆T )

This Work LHR0911, SR18 (+300 K, +2200 K)
Howard et al. (2024) LHR0911, SR18 (−1250 K, +350 K)
Mankovich & Fortney (2020) SR18 (+539 K)
Mankovich et al. (2016) LHR0911 (≲ −100 K)
Püstow et al. (2016) LHR0911 (+500 K)
Nettelmann et al. (2015) LHR0911 (+250 K)
Note—List of H-He miscibility temperature shifts in the literature
needed to match helium abundance observations of Jupiter and Sat-
urn. The temperature shifts of this work as a result of evolving stably
stratified Z profiles are provided in the first row. We note that the
experimental results of Brygoo et al. (2021) show a demixing tempera-
ture of H-He that is ∼5000 K higher than the SR18 demixing temper-
atures and ∼3000 K higher than the LHR0911 demixing temperatures
at megabar pressures.
 11

Figure 3. Effective temperature (top panel) and equatorial

radius (bottom panel) evolution of the models presented in
Figure 2. The models are distinguished by the initial extent
of the stable regions corresponding to the columns in Fig-
ure 2. The deepest fuzzy core model, 0.4 MJ , accesses more
heat than other models at early times. At later times, how-
ever, this model still cannot efficiently access interior heat,
while the other models have convectively mixed and, thus,
their interior heats are made available to their exteriors. The
effective temperature of this model, therefore, is higher than
those of other models at early times, but lower at later times.
The violet horizontal line in the top panel marks the observed
effective temperature of Jupiter of 125.57 K (Li et al. 2012)
and the equatorial radius in the bottom panel, as measured
by (72,492 km; Seidelmann et al. 2007). Due to the high
interior temperatures, the radii remain larger than the equa-
torial radius of Jupiter in these mock evolutionary models.
The purpose of Figure 2 and this figure is not to obtain the
correct radius (for that, see §4 and discussions thereafter),
but rather to illustrate that only extended stable regions lo-
cated at ≲ 0.4 MJ can maintain a semblance of a stable
region.
12

Figure 4. Example evolution models of four mock models with varying initial entropy profiles, increasing from low entropy
(left column) to high entropy (right column). These models were initialized with a steep Z profile gradient surrounded by
flat Z regions, but are provided with stabilizing S profile gradients on either side. To observe the effect of varying the overall
entropy profile, these models all begin with the same Z profile, with a steep gradient fixed initially for all models at 0.55-0.6
MJ . From left to right, this figure shows that higher internal entropies destroy the initial Z profile gradients. As a result of
their lower initial entropies, those that survive experience colder temperatures in the convective region. The first row shows the
inward-moving convective zones as each model cools. The second row shows the evolution of the Z profile. Even a modestly
cold initial interior entropy of 7.5 kB baryon−1 can convectively mix an initially stable region at this location (0.6 MJ ). Higher
interior entropies mix the Z gradients to varying degrees, with the last model on the left mixing completely in under 200 Myr.
The third row depicts the evolution of the helium abundance (Y ) profile, including helium rain, where we have modified the
demixing temperatures of the SR18 miscibility curves by +550 K for example purposes only. Due to conservation of mass and
the setting of the hydrogen/helium mass ratio equal to the protosolar value, regions of high Z correspond to regions of low Y ,
and vice-versa. The last row shows the temperature profiles with the helium rain region shown as the red shaded region. The
two coldest models shown in the left-most columns experience less convective mixing over evolutionary timescales, and their
outside regions are therefore allowed to cool without constant resupply of heat from the interior. This leads to even colder outer
temperature profiles, which eventually intersect the miscibility curve employed here (seen in gray). On the other hand, the two
hottest models shown in the rightmost columns mix completely and remain hotter throughout their interior leading to higher
temperature profiles at later ages.
 13

Figure 5. Evolution of the effective temperature (top panel)

and equatorial radius (bottom panel) of the four mock mod-
els presented in Figure 4. The four models are distinguished
by the entropy at the boundary of the fuzzy core correspond-
ing to the columns of Figure 4. The models initialized with
higher entropies result in hotter effective temperatures and
larger radii since their stable regions are convectively mixed
at earlier ages. For example, the model with interior charac-
teristic entropy of 8.5 kB baryon−1 that the black line repre-
sents mixes entirely in under 200 Myrs. The red line repre-
senting the 7.5 kB baryon−1 model mixes almost entirely by
1 Gyr, after which the effective temperature recovers from
its steep initial dive. The discontinuity jumps visible in the
red line are due to rapid mixing events of the Z profile. The
dotted gray line marks 4.56 Gyr. At later times, the ef-
fective temperature of the coldest model (blue line; 6.5 kB
baryon−1 ) has cooled to under 100 K with a smaller radius
than the rest of the models. This model was initialized with
colder interior temperatures and thus did not convectively
mix. Its convective region therefore evolved without heat re-
supply from the interior leading to low temperatures on the
outside. The violet horizontal line in the upper panel marks
the observed effective temperature of Jupiter of 125.57 K (Li
et al. 2012) and the violet line in the bottom panel the equa-
torial radius of Jupiter (71,492 km; Seidelmann et al. 2007).
14

Figure 6. Example evolutionary models with gradually extending cores. This figure follows the exercise conducted in Figure 4,
starting at a median low interior entropy of 6.5 kB baryon−1 (left) to high median interior entropy of 8 kB baryon−1 (right).
To show helium rain in these models and only for example purposes, we have used here an unmodified LHR0911 miscibility
curve. As in all prior figures, the models have a total heavy envelope mass of 27 M⊕ with a 2.5 M⊕ compact core. The first
row depicts the inward motion of the convective zones as the planet cools. The second row shows the evolution of the Z profile
and, importantly, showcases that higher interior entropies lead to rapid convective mixing of the composition gradient. Cooler
initial entropy profiles let their outer regions cool at a faster rate without access to heat from the interior, leading to encounters
with the LHR0911 miscibility curves at later ages. The third row shows the evolution of the total helium abundance, Y , and
its depletion due to helium rain. Given mass conservation, regions of higher heavy element abundance correspond to regions
of lower helium abundance. The last row shows the evolution of the temperature profile and the helium rain region is shaded
in red. Compared to the models presented in Figure 4, these initial conditions keep the convective zones from moving to the
deep interior regions, leading to stable regions at Gyr timescales, with higher metal abundances in the exterior without mixing
completely.
 15

Figure 7. Evolution of the effective temperature (top panel)

and radius (bottom panel) of the four models presented in
Figure 6, denoted by the corresponding median internal en-
tropies corresponding to the columns of Figure 6. The dotted
gray line marks 4.56 Gyrs. Unlike the models in Figure 4,
the initial distributions of S and Z allow for a more gradual
evolution of the effective temperature even for the 7.5 kB
baryon−1 model, represented by the red curve. This model
retained a stable region out to 0.35 MJ , and came closer to
Jupiter’s effective temperature. The magenta line in the top
panel marks the observed effective temperature of Jupiter of
125.57 K (Li et al. 2012) and the equatorial radius of Jupiter
(71,492 km; Seidelmann et al. 2007) in the bottom panel.
16

Figure 8. Evolutionary model to the present-day that best matches the effective temperature, radius, atmospheric helium
abundance, and outer envelope metallicity within observational uncertainties. The initial mean entropy per baryon is 7.5. The
top left shows the evolution of the entropy profile, the top right the evolution of the composition (Z) profile, the bottom left
the evolution of the helium (Y ) profile, and the bottom right the evolution of the temperature profile. This particular model
contains a total of 44 M⊕ of heavy elements with a 14 M⊕ compact core (Z = 1), marked by the shaded blue regions. As
shown in the upper right panel, the initial Z profile becomes convectively mixed from the outside toward the inside, and this
change is most pronounced at early times. The dark-shaded region represents the normalized radius as a function of normalized
mass, whose coordinates are on the second Y-axis and whose boundaries are color-coded according to the initial and final age.
The resulting extended core (dZ/dm < 0) extends out to 30% of the planet’s mass, corresponding to 42% of the radius at 4.56
Gyrs. The outer metallicity begins at 1.5Z⊙ and rises to 3.26 Z⊙ . Regions of higher metallicity correspond to regions of lower
helium mass fraction. Helium rain begins at 4 Gyrs in this model, for which we use the LHR0911 miscibility curves with a
temperature shift of +350 K to deplete the helium in the outer layers to 0.234. The latter is consistent with the Galileo entry
probe measurements (Von Zahn et al. 1998). The demixing temperatures encountered by the temperature profile are shown in
the last panel. The dotted lines in the last panel show the adjustment of the LHR0911 miscibility curves to their corresponding
temperature and pressure profiles. The long-shaded region corresponds to the helium rain region.
 17

Figure 10. Comparison of the original LHR0911 and SR18

miscibility curves (solid lines) with two Jupiter temperature
profiles at 4.56 Gyrs, including the effects of the inclusion
of heavy elements in the EOS. The low-pressure demixing
temperatures of the SR18 miscibility curves have been lim-
ited to 1.2 Mbar and those of LHR0911 to above 1 Mbar.
The dotted blue curve depicts our favored model using the
LHR0911 miscibility curve, the dotted red curve is a similar
model to the SR18 miscibility curve, and the dashed yellow
model is a fiducial H-He homogeneous Jupiter profile model.
At this age, the outer envelope has a Z of ∼3.3 Z⊙ , mak-
ing the temperatures significantly higher compared to those
of a 1Z⊙ model and more so when compared with the pure
H-He model (dashed yellow). These higher temperatures at
4.56 Gyrs motivate the need for higher miscibility temper-
ature modifications compared to the recent H-He models of
Figure 9. Evolution of the effective temperature (top panel) Mankovich & Fortney (2020, +539 K, SR18) and Howard
and equatorial radius (bottom panel) of the preferred evolu- et al. (2024, -1250 K, LHR0911; +350 K, SR18). The inter-
tionary model of Jupiter in Figure 8. In the top panel, the polated experimental result of Brygoo et al. (2021) is shown
model reaches the observed effective temperature of Jupiter in magenta for comparison. The darker and lighter shaded
(125.57 ± 0.07 K; Li et al. 2012) to within 0.06 K (0.05%) region highlights the lower 1σ and 3 σ uncertainties of that
of the observational uncertainty. The effective temperature experiment, respectively.
of Jupiter is shown as the horizontal magenta line in the top
panel. This model comes to within 0.2% of the equatorial
radius of Jupiter (71,492 km; Seidelmann et al. 2007) at 4.56
Gyrs. The values of Jupiter compared to those of this model
are provided in Table 2. The equatorial radius is shown as
the dash-dot magenta line in the bottom panel.
18

REFERENCES
Bahcall, J. N., Serenelli, A. M., & Basu, S. 2006, The Helled, R., Stevenson, D. J., Lunine, J. I., et al. 2022,
Astrophysical Journal Supplement Series, 165, 400, Icarus, 378, 114937, doi: 10.1016/j.icarus.2022.114937
doi: 10.1086/504043 Howard, S., Müller, S., & Helled, R. 2024, preprint, 1.
Bolton, S. J., Lunine, J., Stevenson, D., et al. 2017, Space http://arxiv.org/abs/2407.11120
Science Reviews, 213, 5, doi: 10.1007/s11214-017-0429-6 Howard, S., Guillot, T., Bazot, M., et al. 2023a, A&A,
Brygoo, S., Loubeyre, P., Millot, M., et al. 2021, Nature, doi: 10.5281/zenodo.7598377
593, doi: 10.1038/s41586-021-03516-0 Howard, S., Guillot, T., Markham, S., et al. 2023b, A&A.
Chabrier, G., & Debras, F. 2021, ApJ, 917, 6pp, https://arxiv.org/pdf/2311.07646.pdf
doi: 10.3847/1538-4357/abfc48 Hubbard, W. B., & Militzer, B. 2016, The Astrophysical
Chabrier, G., Mazevet, S., & Soubiran, F. 2019, ApJ, 872, Journal, 820, 13pp, doi: 10.3847/0004-637X/820/1/80
27pp, doi: 10.3847/1538-4357/aaf99f Iess, L., Militzer, B., Kaspi, Y., et al. 2019, Science, 364,
Chen, Y.-X., Burrows, A., Sur, A., & Arevalo, R. T. 2023, doi: 10.1126/science.aat2965
ApJ, 957, 8, doi: 10.3847/1538-4357/acf456 Keane, A. 1954, Australian Journal of Physics, 7, 322,
Debras, F., & Chabrier, G. 2019, ApJ, 872, 22pp, doi: 10.1071/PH540322
doi: 10.3847/1538-4357/aaff65 Knierim, H., & Helled, R. 2024, arXiv e-prints,
Durante, D., Parisi, M., Serra, D., et al. 2020, Geophysical arXiv:2407.09341, doi: 10.48550/arXiv.2407.09341
Research Letters, 47, 1, doi: 10.1029/2019GL086572
Li, C., Ingersoll, A., Bolton, S., et al. 2020, Nature
Folkner, W. M., Iess, L., Anderson, J. D., et al. 2017,
Astronomy, 4, 609, doi: 10.1038/s41550-020-1009-3
Geophysical Research Letters, 44, 4694,
Li, L., Baines, K. H., Smith, M. A., et al. 2012, Journal of
doi: 10.1002/2017GL073140
Geophysical Research: Planets, 117, 1,
Fortney, J. J., Ikoma, M., Nettelmann, N., Guillot, T., &
doi: 10.1029/2012JE004191
Marley, M. S. 2011, The Astrophysical Journal, 729, 32,
Li, L., Jiang, X., West, R. A., et al. 2018, Nature
doi: 10.1088/0004-637X/729/1/32
Communications, 9, 3709,
French, M., Becker, A., Lorenzen, W., et al. 2012, The
doi: 10.1038/s41467-018-06107-2
Astrophysical Journal Supplement Series, 202, 11pp,
Lorenzen, W., Holst, B., & Redmer, R. 2009, Physical
doi: 10.1088/0067-0049/202/1/5
Review Letters, 102,
Fuentes, J. R., Anders, E. H., Cumming, A., & Hindman,
doi: 10.1103/PhysRevLett.102.115701
B. W. 2023, ApJ Preprint, 18.
—. 2011, Physical Review B, 84,
https://arxiv.org/pdf/2305.09921.pdf
doi: 10.1103/PhysRevB.84.235109
Fuentes, J. R., Hindman, B. W., Fraser, A. E., & Anders,
Lozovsky, M., Helled, R., Rosenberg, E. D., &
E. H. 2024, preprint
Bodenheimer, P. 2017, The Astrophysical Journal, 836,
Fuller, J. 2014, Icarus, 242, 283,
doi: 10.3847/1538-4357/836/2/227
doi: 10.1016/j.icarus.2014.08.006
Fuller, J., Lai, D., & Storch, N. I. 2014, Icarus, 231, 34, Mankovich, C., Fortney, J. J., & Moore, K. L. 2016, ApJ,
doi: 10.1016/j.icarus.2013.11.022 832, 13pp, doi: 10.3847/0004-637X/832/2/113
Fulton, B. J., Petigura, E. A., Howard, A. W., et al. 2017, Mankovich, C. R., & Fortney, J. J. 2020, The Astrophysical
The Astronomical Journal, 154, 109, Journal, 889, 51, doi: 10.3847/1538-4357/ab6210
doi: 10.3847/1538-3881/aa80eb Mankovich, C. R., & Fuller, J. 2021, Nature Astronomy, 5,
Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, 1103, doi: 10.1038/s41550-021-01448-3
SSRv, 123, 485, doi: 10.1007/s11214-006-8315-7 Matson, D. L., Spilker, L. J., & Lebreton, J.-P. 2003, Space
Guillot, T., Fletcher, L. N., Helled, R., et al. 2022. Science Reviews, 104, 1
http://arxiv.org/abs/2205.04100 Mazevet, S., Licari, A., Chabrier, G., & Potekhin, A. Y.
Haldemann, J., Alibert, Y., Mordasini, C., & Benz, W. 2019, A&A, 621, A128,
2020, Astrophysics A&A, 643, 105, doi: 10.1051/0004-6361/201833963
doi: 10.1051/0004-6361/202038367 Militzer, B., & Hubbard, W. B. 2024, Icarus, 411, 115955,
Helled, R., & Howard, S. 2024, preprint. doi: 10.1016/j.icarus.2024.115955
https://ssd.jpl.nasa.gov/planets/phys par.html Militzer, B., Wahl, S., & Hubbard, W. B. 2019, The
Helled, R., & Stevenson, D. 2017, The Astrophysical Astrophysical Journal, 879, 78,
Journal Letters, 840, L4, doi: 10.3847/2041-8213/aa6d08 doi: 10.3847/1538-4357/ab23f0
 19

Militzer, B., Hubbard, W. B., Wahl, S., et al. 2022, The Stevenson, D. J. 1975, Physical Review B, 12. https:
Planetary Science Journal, 3, 185, //journals.aps.org/prb/pdf/10.1103/PhysRevB.12.3999
doi: 10.3847/PSJ/ac7ec8 —. 1982, Planetary and Space Science, 30, 755,
Mirouh, G. M., Garaud, P., Stellmach, S., Traxler, A. L., & doi: 10.1016/0032-0633(82)90108-8
Wood, T. S. 2012, The Astrophysical Journal, 750, 61, Stevenson, D. J., Bodenheimer, P., Lissauer, J. J., &
doi: 10.1088/0004-637X/750/1/61 Angelo, G. D. . 2022, The Planetary Science Journal, 3,
15pp, doi: 10.3847/PSJ/ac5c44
Moll, R., Garaud, P., & Stellmach, S. 2016, The
Stevenson, D. J., & Salpeter, E. E. 1977, The Astrophysical
Astrophysical Journal, 823,
Journal Supplement Series, 35, 221. https:
doi: 10.3847/0004-637X/823/1/33
//articles.adsabs.harvard.edu/pdf/1977ApJS...35..221S
Müller, S., Helled, R., & Cumming, A. 2020, A&A, 638,
Sur, A., Su, Y., Arevalo, R. T., Chen, Y.-X., & Burrows, A.
doi: 10.1051/0004-6361/201937376
2024, The Astrophysical Journal, 971, 104,
Nettelmann, N. 2017, Astronomy and Astrophysics, 606, 1,
doi: 10.3847/1538-4357/ad57c3
doi: 10.1051/0004-6361/201731550
Tejada Arevalo, R., Su, Y., Sur, A., & Burrows, A. 2024,
Nettelmann, N., Amoros, M. C., Tosi, N., Helled, R., & The Astrophysical Journal Supplement Series, 274, 34,
Fortney, J. J. 2024 doi: 10.3847/1538-4365/ad6cd7
Nettelmann, N., Fortney, J. J., Moore, K., & Mankovich, C. Tinetti, G., Eccleston, P., Haswell, C., et al. 2021, arXiv
2015, MNRAS, 447, 3422, doi: 10.1093/mnras/stu2634 e-prints, arXiv:2104.04824,
Nettelmann, N., Movshovitz, N., Ni, D., et al. 2021, The doi: 10.48550/arXiv.2104.04824
Planetary Science Journal, 2, 18pp, Vazan, A., Helled, R., & Guillot, T. 2018, A&A, 610, 14,
doi: 10.3847/PSJ/ac390a doi: 10.1051/0004-6361/201732522
Petigura, E. A., Marcy, G. W., Winn, J. N., et al. 2018, Vazan, A., Helled, R., Podolak, M., & Kovetz, A. 2016,
The Astronomical Journal, 155, 89, ApJ, 829, doi: 10.3847/0004-637X/829/2/118
doi: 10.3847/1538-3881/aaa54c Von Zahn, U., Hunten, D. M., & Lehmacher, G. 1998,
Journal of Geophysical Research: Planets, 103, 22815,
Püstow, R., Nettelmann, N., Lorenzen, W., & Redmer, R.
doi: 10.1029/98JE00695
2016, Icarus, 267, 323, doi: 10.1016/j.icarus.2015.12.009
Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017,
Rauer, H., Catala, C., Aerts, C., et al. 2014, Exp Astron,
Geophysical Research Letters, 44, 4649,
38, 249, doi: 10.1007/s10686-014-9383-4
doi: 10.1002/2017GL073160
Rosenthal, L. J., Knutson, H. A., Chachan, Y., et al. 2021. Wilson, H. F., & Militzer, B. 2010, Physical Review
https://arxiv.org/abs/2112.03399 Letters, 104, 1, doi: 10.1103/PhysRevLett.104.121101
Saumon, D., Chabrier, G., & Van Horn, H. M. 1995, The Winn, J. N., & Fabrycky, D. C. 2015, Annual Review of
Astrophysical Journal Supplement Series, 99, 713 Astronomy and Astrophysics, 53, 409,
Schöttler, M., & Redmer, R. 2018, Physical Review Letters, doi: 10.1146/annurev-astro-082214-122246
120, doi: 10.1103/PhysRevLett.120.115703 Wood, T. S., Garaud, P., & Stellmach, S. 2013, The
Seidelmann, P. K., Archinal, B. A., A’Hearn, M. F., et al. Astrophysical Journal, 768, 157,
2007, Celestial Mechanics and Dynamical Astronomy, 98, doi: 10.1088/0004-637X/768/2/157
155, doi: 10.1007/s10569-007-9072-y Yu, Z. J., & Russell, C. T. 2009, Geophysical Research
Stacey, F. D., & Davis, P. M. 2004, Physics of the Earth Letters, 36, 26, doi: 10.1029/2009GL040094
Zhang, J., & Rogers, L. A. 2022, The Astrophysical
and Planetary Interiors, 142, 137,
Journal, 938, 131, doi: 10.3847/1538-4357/ac8e65
doi: 10.1016/j.pepi.2004.02.003