Accepted for publication in the Astrophysical Journal

Preprint typeset using LATEX style emulateapj v. 5/2/11

BAYESIAN EVOLUTION MODELS FOR JUPITER WITH HELIUM RAIN

AND DOUBLE-DIFFUSIVE CONVECTION
C HRISTOPHER M ANKOVICH , J ONATHAN J. F ORTNEY ,

AND

K EVIN L. M OORE

arXiv:1609.09070v1 [astro-ph.EP] 28 Sep 2016

Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA 95064, USA
Accepted for publication in the Astrophysical Journal

ABSTRACT
Hydrogen and helium demix when sufficiently cool, and this bears on the evolution of all giant planets at
large separations at or below roughly a Jupiter mass. We model the thermal evolution of Jupiter, including
its evolving helium distribution following results of ab initio simulations for helium immiscibility in metallic
hydrogen. After 4 Gyr of homogeneous evolution, differentiation establishes a thin helium gradient below 1
Mbar that dynamically stabilizes the fluid to convection. The region undergoes overstable double-diffusive convection (ODDC), whose weak heat transport maintains a superadiabatic temperature gradient. With a generic
parameterization for the ODDC efficiency, the models can reconcile Jupiters intrinsix flux, atmospheric helium content, and radius at the age of the solar system if the Lorenzen et al. H-He phase diagram is translated
to lower temperatures. We cast the evolutionary models in an MCMC framework to explore tens of thousands
of evolutionary sequences, retrieving probability distributions for the total heavy element mass, the superadiabaticity of the temperature gradient due to ODDC, and the phase diagram perturbation. The adopted SCvH-I
equation of state favors inefficient ODDC such that a thermal boundary layer is formed, allowing the molecular
envelope to cool rapidly while the deeper interior actually heats up over time. If the overall cooling time is
modulated with an additional free parameter to imitate the effect of a colder or warmer EOS, the models favor
those that are colder than SCvH-I. In this case the superadiabaticity is modest and a warming or cooling deep
interior are equally likely.
Subject headings: planets and satellites: physical evolution planets and satellites: interiors planets and
satellites: individual (Jupiter) methods: statistical
1. INTRODUCTION
Cool giant planets are relics of the protoplanetary systems
from which they formed in the sense that they do not fuse protons, and they are well-bound enough that even hydrogen does
not escape appreciably over tens of billions of years. Their
thermal evolution is thus relatively simple, and understanding
it empowers us to use the present states of giant planets to
learn about their history and formation. The open questions
about planet formation thus motivate a comprehensive theory
of giant planet evolution, which will continue to be driven
heavily by our own, well-studied giant planets, Jupiter and
Saturn.
A Henyey-type stellar evolution calculation for a Jupitermass object was first performed by Graboske et al. (1975),
who showed that a convective, homogeneous sphere of fluid
hydrogen and helium could cool to Jupiters observed luminosity over roughly the right timescale, and noted that among
all model inputs, the equation of state (EOS) and superadiabaticity of the temperature gradient have the strongest influence on the overall cooling time.
These two fundamental physical inputs are closely related.
The EOS (paired with a hydrostatic model) is necessary to
translate the planets tangible properties (surface temperature
and composition; external gravity field, size and shape) into
an interior density distribution. Knowledge of the thermodynamic state of matter in these regimes includes understanding any phase transitions that can operate in a Jovian-mass
planets interior, the two most important of which are (i) the
transition from molecular hydrogen to its denser, pressureionized liquid metallic phase, and (ii) the limited solubility of neutral helium in that liquid metallic hydrogen once it
cmankovich@ucsc.edu

cools below a critical temperature (Stevenson 1975). The latter of these two effects has observable ramifications because
the helium-rich phase tends to sink, releasing gravitational energy (constituting a power source beyond mere contraction)
and depleting the outer envelope in helium (Salpeter 1973;
Stevenson & Salpeter 1977). Ultimately a robust theory of
giant planet evolution must reconcile the atmospheric helium
mass fraction Yatm with the helium content of the protosolar
nebula, and this demand constrains the plausible EOS and HHe phase diagram.
Since the critical temperature for H-He phase separation increases with pressure more slowly than the temperature along
a planetary adiabat, the equilibrium helium abundance increases toward the center the planet. Thus in the limit that
the hydrogen-helium mixing ratio is equal to its equilibrium
value throughout the liquid metallic hydrogen part of the mantle, there exists a stabilizing helium gradient that acts to mitigate the convectively unstable temperature gradient. In this
case the dynamics of the fluid (and the degree of macroscopic
vertical heat transport that ensues) are dictated by the competing microscopic diffusion of heat and solute; the fluid is
in the double-diffusive regime. In such a region the temperature gradient can be significantly larger than the adiabatic
gradient, leading to potentially dramatic modifications to the
planets cooling time. For example, double-diffusive convection has been invoked in recent years to explain Saturns
luminosity excess (Leconte & Chabrier 2013; the case of a
global heavy-element gradient), the inflation of hot Jupiters
(Kurokawa & Inutsuka 2015), and Jupiters late thermal evolution including helium rain (Nettelmann et al. 2015), which
we are revisiting in this paper. Although differentiation alone
contributes additional luminosity, extending the overall cooling time, any superadiabatic temperature structure associated

Mankovich, Fortney and Moore

with double-diffusive convection generally cools the surface

more quickly. For models undergoing helium rain, cases with
adiabatic P T profiles thus give an upper limit to the cooling
time (Fortney & Hubbard 2003, Pstow et al. 2016). The inclusion of double-diffusive convection offers a continuum of
shorter cooling times, modulated by the efficiency of the heat
transport through the double-diffusive region.
Nettelmann et al. (2015) sought a solution for Jupiters
evolution to its current state assuming a superadiabatic temperature profile in the framework of layered double-diffusive
convection (LDDC; Mirouh et al. 2012, Wood et al. 2013)
and found that a suitable combination of LDD layer height
and modifications to the H-He phase diagram could match
Jupiters observed 1-bar temperature and Yatm . However, it is
likely that a quasi-stable turbulent state like the layered structures characterized by the direct hydrodynamics simulations
of Wood et al. (2013) would look quite different in the presence of a phase transition and rainout of a main component.
For example, in the context of helium phase separation, homogeneous layers themselvesfinite volumes of (P, T ) space, at
effectively uniform helium abundanceare intrinsically unstable to H-He phase separation, and the influence that droplet
formation and rainout has on the merging or bifurcation of
convective layers, or the transport of solute between layers,
has yet to be assessed from the hydrodynamical perspective.
The present work thus relaxes the assumption of layered convection, opting instead to treat the superadiabaticity with a
generic parameterization, the only physical content of which
is the demand that the temperature gradient lie somewhere
between the minimum value for overstable double-diffusion
and the upper limit imposed by the Ledoux criterion. This
amounts to the criterion that gravity waves be linearly overstable, so that thermal transport is enhanced relative to the
purely diffusive case by some degree of turbulence.
Despite growing confidence that helium has begun differentiating in Jupiters recent past (and billions of years ago in
the case of Saturn; see Fortney & Hubbard 2003 and Pstow
et al. 2016), it is not known whether helium rain alone can resolve the gaps in our understanding of giant planet evolution
given Jupiter and Saturns luminosities and the helium content of their molecular envelopes at the present day. Although
the thermodynamic conditions for phase separation of helium
from liquid metallic hydrogen have been evaluated since the
early work of Stevenson (1975) and Straus et al. (1977), quantitative knowledge covering the relevant pressures and H-He
mixing ratios has only become available over the past several years as a result of ab initio simulations making use of
density functional theory molecular dynamics (Morales et al.
2009; Lorenzen et al. 2009, 2011; Morales et al. 2013). The
present work demonstrates that using ab initio results for the
H-He phase diagram, a differentiating non-adiabatic Jupiter
comprised of hydrogen and helium surrounding a dense core
of heavy elements explain Jupiters evolutionary state at the
solar age.
To assess the viability of the evolution models, we formulate the problem in a Bayesian framework, using Jupiters observed Teff , Yatm , and volumetric mean radius Rvol to derive
posterior probability distributions for the model parameters
using a Markov chain Monte Carlo sampling algorithm. Most
importantly, we make a probabilistic determination of the superadiabatic temperature gradient to be expected in the deep
interior, and simultaneously estimate the temperature correction that must be applied to the Lorenzen et al. (2011) phase
diagram to satisfy the Galileo entry probe measurement of

Yatm (Nettelmann et al. 2015). The present work thus extends

the basic approach of Fortney & Hubbard (2003)using forward thermal evolution models to infer a most likely H-He
phase diagramwith the power of a Bayesian parameter estimation method and the treatment of non-adiabatic P T profiles.
In 2 we describe our modeling approach using Modules
for Experiments in Stellar Astrophysics (MESA; Paxton et al.
2011, 2013, 2015), including our atmospheric boundary condition, treatment of helium phase separation, thermal transport, and other modifications that were necessary for our application. We describe the three free parameters in our inhomogeneous, non-adiabatic models, namely the heavy-element
core mass Mc , the double-diffusive superadiabaticity (or density ratio) R , and the phase diagram temperature offset
Tphase . In 3 we present results of evolutionary calculations, first validating our models for the case of homogeneous composition, then discussing in detail the late inhomogeneous, non-adiabatic evolution as a result of helium rain.
We then repeat these calculations, but treating the planets
equilibrium temperature Teq as a fourth free parameter controlling the overall cooling time, mimicking the influence of
a colder or warmer EOS than the adopted Saumon et al.
(1995) EOS. Marginalizing over this parameter allows us in
an indirect sense to marginalize over the plausible H-He equations of state and thus obtain the most general estimates for
the remaining three parameters. In 4 we summarize and contextualize our findings.
2. PLANETARY EVOLUTION MODELS
Our evolution models are computed using Modules for
Experiments in Stellar Astrophysics (MESA; Paxton et al.
2011, 2013, 2015). The models are hydrostatic, nonrotating
spheres with envelopes consisting of binary mixtures of 1 H
and 4 He surrounding dense inert cores of heavy elements. In
the density-temperature regime relevant to giant planets with
M . MJ , MESA employs the Saumon et al. (1995) equation of
state, interpolated over hydrogens molecular-metallic phase
transition such that the density varies smoothly between the
two phases (SCvH-I). This EOS is advantageous for studies
of helium phase separation because it provides the necessary
state variables for arbitrary mixtures of hydrogen and helium,
which is critical for solving the energy equation throughout
the interior of the differentiating planet. This is one of our
principal motivations for using SCvH-I over more recent HHe EOSs obtained with ab initio methods (e.g., Militzer &
Hubbard 2013) in spite of the latter class comparing more favorably with shock experiments. Rosseland mean radiative
opacities are taken from Freedman et al. (2008) and a privately communicated 2011 update. Electron conduction opacities are based on Cassisi et al. (2007).
Jupiter is oblate as a result of its rapid rotation, while our
present models are perfect spheres. A suitable mean planet radius with which to compare our model radii is the volumetric
radius (Seidelmann et al. 2007)
1/3

Rvol R2/3
eq Rpol = 69911 6 km,

(1)

where Req and Rpol are Jupiters equatorial and polar radii at
1 bar; Rvol is the radius of sphere enclosing the same volume as does Jupiters 1-bar surface. Because the overall compactness of the planet is steeply sensitive to its heavy element mass, the high precision of this radius measurement
translates into an extremely narrow range of allowed core

masses. As an example, fitting the radius of our homogeneous, adiabatic Jupiter modelfor which Mc is the sole free
parameterto Jupiters Rvol at the solar age using MCMC
produced Mc = (25.33 0.03) ME (the quoted value corresponding to the median and the error to the 68% confidence
interval). Our models prefer large core masses because of the
assumption that all heavy elements are relegated to a dense
core; in reality, at least half of Jupiters heavy element mass
probably resides in the hydrogen-dominated mantle and envelope (Saumon & Guillot 2004). Thus in our findings (3),
Mc should be interpreted as a total heavy element mass. Indeed, the fact that our adopted equation of state is limited to
hydrogen and helium is the reason why we make no effort to
calculate the oblateness and associated gravitational multipole
moments (J2 , J4 , . . .) for our models.
The presence of heavy elements in the envelope and the action of rotation would both modify the hydrostatic structure
and thus the total cooling time obtained for an evolutionary
model. To assess the sensitivity to the heavy element distribution, we computed models with fixed total heavy element
mass MZ = Mc + MZ,env = 28 ME and varying combinations of
the core mass Mc and envelope heavy element content MZ,env ,
modelling the heavy elements in the H-He envelope simply by
taking Y Y + Z for the purposes of this diagnostic. We find
that models with more of their heavy elements mixed throughout the envelope cool more quickly (e.g., Baraffe et al. 2008),
with a characteristic spread of 70 Myr in the time for models with no H-He phase separation to cool to Jupiters volumetric radius. This is a relatively short time compared to the
2 5 108 yr spread in cooling times obtained by varying the total heavy element content or atmospheric boundary
condition, as discussed in 3. The centrifugal support provided by rotation would also modify the radius evolution, but
as there are a number of complexities associated with rotation
(e.g., the likely differential rotation as a function of radius or
latitude, and the details of the planetary figure as a result of
even rigid rotation), we do not model rotation here. It should
be noted that the effect of rotation is to some degree degenerate with the total amount and distribution of heavy elements,
since the tendency for centrifugal acceleration to prefer larger
radii can be offset by incorporating more heavy elements.

. 3000 K, which on a present-day Jupiter (or Saturn) adiabat

is well outside the transition from molecular to metallic hydrogen predicted by either group. Although it appears from
their Figure 1 that the phase diagram of Morales et al. (2013)
predicts that Jupiter has not yet cooled into the immiscibility
region of (P, T ) space for a protosolar mixture, the favored
model of Hubbard & Militzer (2016) has a cooler deep interior such that it has already undergone H-He phase separation
according to both phase diagrams. These cooler deep layers
are a result that Militzer et al. (2008) attributed to the important inclusion of the nonideal entropy of mixing between H
and He.
Apart from the major differences in the two modern phase
diagrams behavior at T . 3000 K, P . 1 Mbar, their overall
shapes at the warmer conditions relevant to Jupiter and Saturn
are in rough agreement, with the major effect being an temperature offset (up to 1000 K) between the two. It is possible
that this discrepancy stems from including versus excluding
the nonideal entropy of mixing between the two components.
It was demonstrated by Nettelmann et al. (2015) that applying the raw phase diagram of Lorenzen et al. (2011) to Jupiter
resulted in too large a loss of helium from the molecular envelope, motivating a global downward offset in the demixing temperature such that the onset of phase separation takes
place later in the planets history. It is reassuring that an offset
motivated by Jupiters observed atmospheric helium content
brings the two phase diagrams into closer agreement. In this
paper, we also allow for a global temperature offset Tphase
to this phase diagram and estimate its value from available
data; we confirm the result of Nettelmann et al. (2015) that
for models computed with SCvH-I, the necessary downward
offset is of order 200 K. Shifts to the published phase diagram
in pressure are also conceivable, but as the phase diagrams of
Lorenzen et al. (2011) and Morales et al. (2013) agree fairly
closely on the minimum pressure for phase separation at temperatures typical of Jupiter and Saturns molecular-metallic
transitions, the temperature offset appears to be the most important correction.

2.1. Initial conditions

Adiabatic, extended (1-bar radius R = 2 RJ ) initial models of
mass M = MJ = 1.89861 1030 g were created using MESAs
create_initial_model capability described in Paxton
et al. (2013). For each model, a fraction of the total mass was
converted into a dense, nonevolving core of specified mass Mc
and density c . All models in the present work adopt a constant core density c = 15 g cm3 for simplicity. The evolution
is insensitive to the particular choice of core density; the total
radius R is set by Mc nearly independent of c .

T (K)

Helium Rain in Jupiter

2.2. Hydrogen-helium phase separation

We couple our evolutionary calculations to the H-He phase
diagram of Lorenzen et al. (2011), whose calculations include
the simplifying assumption of ideal entropy of mixing between the two species. The group of Morales et al. (2013)
instead performed direct thermodynamic integrations, thus including nonideal contributions to the entropy of mixing, but
at the expense of subtantially more sparse sampling in helium
number fraction xHe . The results of the two groups are in reasonable agreement, diverging most noticeably at temperatures

12000
10000
8000
6000
4000

1 Mbar
2 Mbar
4 Mbar
10 Mbar

2000
0.00 0.05 0.10 0.15 0.20

He number fraction xHe

immiscibility region
for xHe =0.0866
0

1 2

3 4

5 6

P (Mbar)

F IG . 1. H-He phase diagram of Lorenzen et al. (2011), illustrating immiscibility regions as a function of helium number fraction xHe for four pressures
(left panel), and as a function of pressure for the protosolar mixture (right
panel). Immiscibility regions (shaded regions bounded by thick curves) are
precluded in terms of equilibrium thermodynamics, i.e., by the criterion that
the Gibbs free energy be stable with respect to perturbations in the helium
concentration. The vertical dashed line in the left panel designates the protosolar mixture, xHe = 0.0866 (Y = 0.275). Open circles in the right panel are
raw data and the curve is a linear interpolation in log P. A solar-age profile of
a representative differentiated Jupiter (Mc = 30 ME , R = 0.25, Tphase = 0)
is shown in the thin black solid curves.

The phase diagram of Lorenzen et al. (2011) adopted for

the models in this paper is illustrated in Figure 1, wherein
phase curves are shown in (xHe , T ) space (left panel) and
in (P, T ) space (right panel). The data are finely sampled

Mankovich, Fortney and Moore

in xHe but more coarsely in P, where data are available at

P = (1, 2, 4, 10, 24) Mbar. Our model assumes no phase
separation at pressures lower than 1 Mbar where no data
are available. This assumption is reasonable given the nearvertical (P, T ) phase curves found by both Lorenzen et al.
(2011) and Morales et al. (2013) for the relevant temperatures
4 kK < T < 6 kK, although the latter phase diagram situates
the boundary at slightly lower pressure P 0.8 Mbar. To
compute the equilibrium abundances for zones at arbitrary T
and P, we interpolate in the tables linearly in xHe and log P to
obtain Tphase , and numerically solve the equation
Tphase (xHe , P) + Tphase T = 0

(2)

for the equilibrium helium concentration xHe .

Differentiation begins once the planetary P T profile cools
into the immiscibility region for the protosolar mixture, at
which time one or more grid points are supersaturated in helium. For this and each subsequent timestep, we assume the
thermodynamic equilibrium distribution of helium throughout
the interior, which amounts to the assumption that all helium
excess is delivered efficiently by gravitational settling to lower
depths such that it can exist in chemical equilibrium with its
surroundings. Here efficient delivery means that the local
supersaturation can be reduced to zero faster than the other
relevant timescales, namely the planets thermal timescale,
the large-scale convective circulation time, and the simulation
timestep.
The thermal timescale th = Egr /L, with Egr denoting the
total binding energy, is roughly 109 yr. Timesteps in our simulations are typically 106 to 107 yr. The convective circulation
time presents the strongest condition: it can be estimated with
mixing length theory as (Guillot et al. 2004) MLT 108 s, or
3 yr. Stevenson & Salpeter (1977) derived the minimum size
a He-rich droplet would need to attain in order for its terminal
speed to exceed the average speed of convective motions of
the ambient fluid. Balancing the terminal speed from a turbulent drag law with an average convective speed estimated
from mixing length theory, those authors obtained a minimum
droplet length scale on the order of 1 cm. Repeating their calculation, we balance the (downward) buoyancy force with a
drag force through a convective plume:
d 3 g Cd 2 v2 ,

(3)

where is the density excess compared to the surrounding

H-dominated fluid, d is the length scale of the droplet, C is
a nondimensional drag coefficient, and v is the relative speed
between the droplet and the medium. A minimum droplet size
dmin corresponds to a droplet having v = vMLT so that
dmin

Cv2MLT
,
g

(4)

where vMLT 10 cm s1 , , and g 3 103 cm s2 . If

as in Stevenson & Salpeter (1977) we assume turbulent flow
over the scale of a droplet (Re vd/ & 103 ), then C = 0.5
(Landau & Lifshitz 1959) and we arrive at a smaller minimum
size dmin 102 cm; the discrepancy may be the result of a
dimensional error in the earlier calculation. However, given
a typical kinematic viscosity = 4 103 cm2 s1 (French
et al. 2012), the Reynolds number Re ud/ 25 so that a
turbulent drag force is not appropriate. If instead we balance
buoyancy with a Stokes drag appropriate for low Reynolds

number, then (Landau & Lifshitz 1959)

3
dmin
g 6dmin vMLT

(5)

and we obtain

dmin =

6vMLT
g

1/2
102 cm,

(6)

incidentally the same estimate as in the turbulent case.

Again following Stevenson & Salpeter (1977), this droplet
size can be translated into a droplet formation timescale if we
suppose that the droplets grow by microscopic diffusion of He
nuclei into the He-rich pockets. Assuming the microscopic
diffusivity DHe to be of order the He-He self-diffusion coefficient obtained in the recent ab initio simulations by French
et al. (2012), the timescale to form a sinkable droplet is
sedimentation =

2
dmin
(102 cm)2
101 s.
= 3
DHe 10 cm2 s1

(7)

The strong hierarchy in timescales

sedimentation  MLT  th

(8)

implies that the transport of excess helium toward the center of the planet is probably efficient in spite of convection
(or double-diffusive convection, which mixes material over
much longer timescales), so assuming the equilibrium mixture throughout the planet at each timestep is an adequate
starting point for evolutionary models.
For all timesteps in which a cell has been cooled below
its critical temperature, we lower the abundance of the outermost supersaturated grid point (i.e., the first grid point with
P > 1 Mbar) to its equilibrium value. We apply this same
equilibrium abundance throughout the molecular envelope,
reflecting the fact that outside of regions where phase separation is taking place, the species are being rapidly mixed by
convection. Since a single timestep corresponds to millions of
large-scale convection cycles, the entire molecular envelope
acts as a reservoir of helium for the phase separation and rainout taking place near the molecular-metallic transition. The
molecular envelope thus depletes uniformly in a given step.
Iterating inward over grid points, we enforce the local equilibrium abundance (a monotonically increasing function of
depth) in each cell, and propose that same abundance as
the tentative (He-enriched) mixture to be applied to the remainder of the metallic interior. Eventually, a grid point is
reached whose equilibrium abundance is greater than or equal
to its proposed abundance. At this point the inward iteration
ceases and the homogeneous, He-enriched interior is stable to
phase separation; all that remains is to enforce conservation
of helium overall. The proposed profile always has a helium
deficit, which is resolved by the constraint that all helium nuclei that have rained out from above are mixed into the homogeneous interior; we raise the helium abundance of the homogeneous interior accordingly. Since this adjustment typically leaves the uppermost layers in the homogeneous interior
marginally supersaturated, we repeat the iteration over all grid
points as many times as necessary to achieve the equilibrium
profile. In practice this takes one or two more iterations.
2.3. Modes of heat transport
The pioneering work of Hubbard (1968, 1969, 1970) made
the case for Jupiters envelope as a convective fluid of primarily hydrogen and helium. The short mean free path of

Helium Rain in Jupiter

photons and electrons in the interior imply a small thermal

conductivity, with the result that convection, rather than radiation or conduction, carries nearly all the intrinsic flux.
The temperature gradient is thus marginally superadiabatic;
estimates from mixing length theory yield ad . 106
throughout the envelope and ad . 109 within hydrogens molecular-metallic phase transition at megabar pressures. Our homogeneous Jupiter models indicate that convection maintains these small superadiabaticities to within an
order of magnitude over the history of the planet with the exception of an early, brief radiative window similar to that described by Guillot et al. (1995). As demonstrated by Guillot et al. (2004), the inclusion of alkali metals enhances the
opacity enough to ensure convection at all depths within a
present-day Jupiter. Our models, which include modern opacities and self-consistently allow for radiative transport wherever rad < ad , confirm this result: the intermediate radiative window is only mildly subadiabatic ( ad 101 )
and vanishes before t 2 108 yr (see also Fortney et al.
2011), and has an insignificant effect on the overall cooling
time. Thus for the homogeneous phases of the evolution, our
models may be directly compared to models constructed assuming = ad always.
In the more general case allowing for stratification in the
mean molecular weight , the Schwarzschild & Hrm (1958)
criterion for dynamical convection

instability reduces to Equation 9 and marginally superadiabatic temperature gradients can be sustained by convection.
Stevenson (1979) came to the same conclusion, arguing that
the otherwise weak vertical heat transport provided by these
overstable oscillations is mitigated by the occasional breaking
of waves (Rosenblum et al. 2011), redistributing solute such
that = ad to a good approximation.
Metallic hydrogen environments in cool giant planets differ from the stellar case for two reasons: (i) the formation,
rainout, and deeper redissolution of helium droplets tends to
enforce a persistent, stabilizing composition gradient, and (ii)
since the (conduction-limited) diffusion of heat is quite inefficient, overstable gravity waves have relatively slow growth
rates so that wave-breaking events are rare and the fluid is only
weakly turbulent; vertical heat transport is thus enhanced relative to the purely diffusive case but is still much weaker than
in the case of overturning convection. The resulting temperature gradient is substantially superadiabatic, possibly closer
to the Ledoux limit = ad + (/) (Stevenson 1979;
Mirouh et al. 2012).

> ad

ad
= ad + R R =
.

(/)

(9)

must be replaced by the Ledoux (1947) criterion

> ad + ,
(10)

where
d ln

(11)
d ln P
is the slope of the mean molecular weight along the planetary profile, and and are two thermodynamic derivatives
defined by




ln
ln
, =
.
(12)
=
ln P,T
ln T P,
(Details on the novel calculation of (/) in MESA are
given in Paxton et al. 2013 3.3). If mean molecular weight
increases toward the planets center, as is the case for a differentiated planet, then > 0 and the regime

ad < < ad +
(13)

corresponds to the situation wherein a superadiabatic temperature profile is dynamically stabilized by the chemical stratification. The mixing that ensues in this regime is termed semiconvection in the stellar context (Schwarzschild & Hrm
1958; Sweigart & Gross 1974) and double-diffusive convection in the hydrodynamic context (Turner 1974). In this case
the Brnt-Visl frequency N of the fluid is real-valued, admitting gravity waves. These modes are in general overstable
if fluid parcels can exchange a significant fraction of their
heat with the environment over an oscillation period. The
linear stability analysis of Kato (1966) demonstrated that in
the stellar case, where radiative diffusion is efficient, the thermal diffusion timescale tends to be short compared to buoyant
oscillation periods N 1/2 so that the criterion for convective

2.4. A parametric model for double-diffusive convection

We assume the temperature gradient to be adiabatic unless there exists a stabilizing composition gradient > 0, in
which case the temperature gradient is steeper than the adiabat
by an amount proportional to , i.e.,
(14)

Equation 14 defines the density ratio R , which describes the

relative (and competing) contributions that the temperature
and composition stratifications make to the overall density
stratification. We take R as a free parameter that we seek
to estimate. Although the commonly cited criterion of Equation 13 is a necessary condition for semiconvection, it is not
sufficient; the linear instability demands the somewhat more
strict criterion (Walin 1964; Kato 1966)
Pr +
Rcrit .
(15)
Pr + 1
Thus for the semiconvective instability to grow, the temperature gradient must be superadiabatic by a nonvanishing
amount determined by the Prandtl number Pr and diffusivity
ratio defined by

Pr =
,
=
,
(16)
T
T
1 > R >

where is the fluids kinematic viscosity, T is its thermal

diffusivity, and is the diffusivity of solute, in this case the
diffusivity of helium atoms in a mixture which is predominantly metallic hydrogen. In Figure 2 we show values of Pr,
, and Rcrit , derived from the ab initio transport properties obtained by French et al. (2012) for a Jupiter adiabat. Here the
calculation of (and thus Rcrit ) assumes an effective composition diffusivity that is of order the He-He self-diffusion
coefficient reported in that paper. The values of Rcrit indicate
that overstable modes can grow at 1 Mbar and deeper for density ratios R > 101 .
2.5. Energetics of evolving compositions
Enriching or depleting a Lagrangian fluid element in helium generally modifies its internal energy per gram u and

Mankovich, Fortney and Moore

100

10

10

10

Rcrit
Pr

overstable
g-modes

100

101

102

P (Mbar)
F IG . 2. Estimates of the dimensionless quantities Pr and (Equation 16), as well as the critical density ratio Rcrit for overstable doublediffusive convection (Equation 15). Quantities are derived from the ab initio
transport properties of French et al. (2012) for the metallic hydrogen part of
Jupiters interior. The shaded region is the intersection of Rcrit < R < 1 and
P > 1 Mbar, within which the stable stratification from helium rain admits
growing-amplitude gravity waves and thus some degree of double-diffusive
convection.

does work modifying its density following the fundamental

thermodynamic relation
T

ds du
d(1/)
=
+P
,
dt dt
dt

(17)

where s is the entropy per gram and d/dt denotes a time

derivative. This change in heat content is commensurate with
the energy gained or lost by the fluid element via photons:
ds
dL
= T
(18a)
dm
dt
du
d(1/)
= P
.
(18b)
dt
dt
Here L is the local luminosity and m is the Lagrangian coordinate; other energy sources and sinks (nuclear fusion or
fission, tidal dissipation, neutrino cooling) are negligible for
our application. While the two forms of the energy equation Eqs. 18a and 18b are fundamental, MESAs solvers do
not work with the entropy directly, and also do not work with
(u, ) by default (although the latter option exists). In either of
the two standard thermodynamic bases (, T ) or (P, T ), these
two variables (along with m, r, and Xi ) are solved for simultaneously and then s or u are computed from the EOS post-hoc.
Because s or u are not solved for directly, their finite differences over time are subject to numerical noise, and direct finite differences of either form 18a or 18b thus yield noisy
luminosity profiles. The approach MESA takes by default is to
instead recast Equation 18 into time derivatives of the quantities comprising the adopted thermodynamic basis. Since the
phase diagram of Lorenzen et al. (2011) provides the equilibrium helium abundance over (P, T ) space, and these are also
the independent variables in the Saumon et al. (1995) EOS,
we choose to adopt (ln P, ln T ) as the thermodynamic basis
for our MESA calculations. In this basis the energy equation
takes the form (Paxton et al. 2013)



dL 
d ln T
d ln P
=
c
T

,
(19)
P
ad
dm Xi
dt
dt
where we have assumed radiation pressure is negligible as is
appropriate for T . 104 K.

The standard transformation of Equation 18 into Equation 19 (e.g., Kippenhahn & Weigert 1990) ignores the fact
that entropy and internal energy depend not only on P and T
but also on the composition vector Xi , and hence Equation 19
is only accurate at fixed composition. This poses no substantial problem for energy conservation in stellar models, where
large abundance changes typically only happen as a result of
fusion, in which case nuclear energy generation overwhelms
the T (ds/dt) term in the energy equation; one important exception is the accretion of material with a composition different from the stellar surface. For our application, it is necessary to add to Equation 19 the component of dL/dm that
arises from composition changes at fixed P and T :

dL 
u dXi
(1/) dXi
=
P
.
(20)
dm P, T
Xi dt
Xi dt
Here the repeated indices denote summation over species
i = 1, . . . , N 1 where N is the total number of species in the
model. (Since all N mass fractions sum to unity, only N 1
mass fractions are independent.) All models in this work assume a two-component mixture of 1 H and 4 He, so that Equation 20 reduces to just

u dY
(1/) dY
dL 
=
P
.
(21)

dm P, T
Y dt
Y dt
In practice we calculate this term for each cell as

dL 
u(P, T,Y0 ) u(P, T,Y1 )
=
dm P,T
t
 1

(P, T,Y0 ) 1 (P, T,Y1 )
P
t

(22)

where Y0 and Y1 denote the helium mass fractions before and

after our helium redistribution step, and t denotes the (finite)
timestep. Obtaining u(P, T,Y1 ) and 1 (P, T,Y1 ) requires one
additional call to the EOS module per zone per timestep.
2.6. Model atmospheres
The overall cooling time of an evolutionary giant planet
model depends strongly on the boundary condition applied at
the atmosphere, since that condition determines how rapidly
the planet can radiate. We apply the self-consistent model
atmospheres of Fortney et al. (2011) as fit analytically by
Leconte & Chabrier (2013), which provide the temperature
at the 10 bar level T10 as a function of surface gravity g and
intrinsic temperature Tint . The planets effective temperature
Teff in a given timestep is given by
4
4
Teff
= Tint
+ Teq4 ,

(23)

where the equilibrium temperature Teq describes the orbitaveraged temperature of a body radiating as much energy as
it absorbs from the Sun following
(1 A)L
,
(24)
16a2
with L the instantaneous luminosity of the star and a the
planets orbital semimajor axis. Effective temperatures for our
models are calculated assuming the value for A determined
by Voyager measurements in the infrared (Hanel et al. 1981;
Pearl & Conrath 1991). Although the planets albedo is certainly a function of time as a result of changing atmospheric
Teq4 =

Helium Rain in Jupiter

Teff (K)

138
Mc =20 ME
136
Teq =109.0 1.4 K
134
132
130
128
126
124
122
3.0 3.5 4.0 4.5 5.0 5.5

Age (Gyr)

Teq =109.0 K
no core
Mc =10 ME
Mc =20 ME
Mc =30 ME

3.0 3.5 4.0 4.5 5.0 5.5

Age (Gyr)

F IG . 3. Sensitivity of the cooling time to the surface boundary condition

(left panel) and heavy-element mass (right panel) for homogeneous, adiabatic
Jupiter models. The dashed curve corresponds to Teq = 109.0 K, and the blue
(magenta) curve corresponds to plus (minus) Teq = 1.4 K.

3. RESULTS OF EVOLUTIONARY CALCULATIONS

We first validate our general modeling approach and implementation of the model atmospheres by computing homogeneous, adiabatic evolutionary sequences. Figure 3 shows
evolution in the age-Teff plane for models with core masses
between 0 and 30 ME and a range of assumed equilibrium
temperatures reprensenting the uncertainty in the Voyager determination of Jupiters Bond albedo. Models with higher
equilibrium temperatures generally take longer to cool because they absorb more stellar flux, and models with greater
heavy element content generally cool faster because they are
more compact. Figure 4 provides a summary of cooling times
attainable by homogeneous models across Mc Teq parameter
space, and demonstrates that over the cooling time depends
more steeply on the atmospheric boundary condition than on
the assumed heavy element mass. The total cooling times
agree closely with published results also using the SCvHI EOS (Fortney & Hubbard 2003; Saumon & Guillot 2004;
Fortney et al. 2011).
The SCvH-I EOS generally leads to slow cooling for the homogeneous models. These were evolved to Jupiters observed
Teff = 124.4 K, a temperature notably not reached within the
age of the solar system (4.56 Gyr) for any of the models,

110

Voyager

106
104
102
100
0

5.8
5.6
5.4
5.2
5.0
4.8
4.6
4.4
4.2
4.0

cool (Gyr)

108

Teq (K)

dynamics and chemistry (Hubbard & Smoluchowski 1973),

assuming the present-day value of the albedo throughout the
evolution is an acceptable (and at this stage necessary) approximation.
While Voyager 1 measured a Bond albedo A = 0.343
0.032, corresponding to Teq = 109.9 1.3 K, Pearl & Conrath (1991) noted that the value of A determined using Voyager 2 radiometry is larger than that determined by Voyager 1
by roughly 12%. As those authors suggested, if this discrepancy is the result of unidentified systematic error such that
the true value is 6 6% larger than the Voyager 1 measurement, then a revised estimate for Jupiters Bond albedo is A =
0.366 0.035. This corresponds to a modestly smaller equilibrium temperature Teq = 108.5 1.4 K. In the present work
we adopt the median equilibrium temperature Teq = 109.0 K
of the two Voyager determinations as ground truth. In 3.2
we also show our full calculations repeated with Teq as a
free parameter to address the broad modeling uncertainties
principally those associated with the EOSthat contribute
to the uncertain overall cooling time for a contracting giant
planet.

10

15

Mc (ME)

20

25

30

F IG . 4. Time for homogeneous 1.0 MJ models to cool to Jupiters Teff as

a function of their heavy-element mass Mc and equilibrium temperature Teq .
The horizontal lines designate the Voyager measurement of Teq (see 2.6).
The color scale is piecewise linear such that yellow corresponds to the solar
age 4.56 Gyr.

even those with vanishing heavy element mass. It is clear that

whatever superadiabaticity arises from helium rain in the inhomogeneous case must act to accelerate the planets cooling
in spite of any additional luminosity associated with differentiation. This is the first indication that inhomogeneous models
computed with the SCvH-I equation of state will tend to favor
fairly weak heat transport in the helium gradient region such
that the temperature distribution is strongly superadiabatic
and a thermal boundary layer is established. As described by
Stevenson & Salpeter (1977), in this case the molecular envelope can cool rapidly while the cooling of the deeper interior
is stalled or even reversed. Indeed, sections 3.1 and 3.2 below
demonstrate that best-fitting models have a steep enough temperature gradient in the stable region that the cooling of the
molecular envelope is accelerated while the helium-enriched
metallic interior is heated over time.
Some more recent equations of state based on ab initio methods (e.g., Militzer & Hubbard 2013) yield adiabats
which are colder at Mbar pressures and Jupiter- and Saturnlike entropies than the adiabats obtained from the semianalytic model of Saumon et al. (1995). These colder equations of state compare more favorably with shock experiment,
and generally predict shorter cooling times for homogeneous
Jupiters because their adiabats have less total internal energy
for a given global entropy. With such an EOS, if inhomogeneous evolution is to help with addressing Jupiters luminosity
constraint, then any double-diffusive convection in the deep
interior must act to lengthen the cooling relative to the homogeneous case. Thus for a colder equation of state, more
modest superadiabaticities should be expected such that the
differentiation luminosity overwhelms the accelerated cooling of the envelope due to the double-diffusive bottleneck at
1-2 Mbar. This situation is closer to the luminosity problem for Saturn, where the drastic underluminosity of homogeneous models is robust with respect to the assumed equation
of state. In 3.1 we illustrate the central role that the efficiency of heat transport by ODDC plays in determining the
thermal evolution. In 3.2 we retrieve strong superadiabaticities for our nominal SCvH-I case. Then to address the systematic modelling uncertainty associated with the H-He equation
of state, we repeat our calculations with Teq taken as a fourth

Mankovich, Fortney and Moore

free parameter to modulate the overall cooling time, with high

(low) Teq mimicking the effect of a warmer (colder) equation
of state.
3.1. Inhomogeneous evolution
Interior profiles for Jupiter models undergoing helium rain
are illustrated in Figure 5 for four different values of R ,
including the adiabatic case R = 0. (In the case of evolving composition profiles, ad = ad (P, T, Xi ) where the Xi
are now functions of depth. We refer to the case where
= ad everywhere as adiabatic, although the profiles are
emphatically not isentropic.) For this figure and all others
in this section, unless otherwise indicated, the two remaining free parameters are arbitrarily chosen as Mc = 30 ME and
Tphase = 0 K for illustrative purposes. The first column shows
profiles at 3.5 Gyr, at which time helium rain has not yet commenced and the models are thus still in an identical state. The
remaining three panels show profiles shortly after the onset of
helium rain (3.8 Gyr), then at a more typical time after helium rain onset (4.0 Gyr), and finally at the solar age (4.56
Gyr). The inset in the center left panel plots Teff as a function
of age for the four models.
Evolution in the adiabatic case R = 0 is simplest because
there is minimal feedback between the evolving composition and temperature profiles, and thus the shape of L(m)
can be understood by inspecting just the composition term of
dL/dm (Equation 21); the thermal term retains essentially the
same smooth profile as before helium rain sets in. The layers with decreasing helium abundancefrom the planetary
surface down to the lower boundary of the helium gradient
regionhave an increasing internal energy u and specific volume 1/, and thus dL/dm < 0 there. Similarly, layers deeper
than the bottom of the gradient region have a uniformly increasing helium abundance and thus dL/dm > 0 there. Hence,
a global maximum in the luminosity is attained at the base
of the gradient region. The adiabatic model also exhibits the
most extended helium gradient region among the models considered (spanning 1 to 2.4 Mbar at the solar age), simply because a shallower T (P) profile intersects the immiscibility gap
over a broader range in pressure. Likewise, larger values of
R represent steeper T (P) profiles, which generally intersect
the immiscibility gap over a more narrow range in pressure.
This can be seen by comparing the two solar-age profiles in
the right panel of Figure 1, and is manifested in the relative
widths of the features in Y and in Figure 5.
The evolution is more complex in the case of nonzero superadiabaticity R > 0, owing to the feedback between the
composition profile and the temperature profile. Since in this
case the settling of helium directly modifies T (P) via Equation 14, it contributes to both the thermal (Equation 19) and
composition (Equation 21) components of dL/dm. For cases
with R & 0.20, double-diffusive convection poses a sufficiently strong thermal barrier that the homogeneous, adiabatic
deep interior actually heats up with time, and the reduced heat
flux impinging on the bottom of the molecular envelope allows the envelope to cool relatively quickly. The temperature evolution in the vicinity of the helium gradient is illustrated for these same four values of R in Figure 6, and Figure 7 shows the evolution of the core-mantle boundary in T
space. Indeed since a more pronounced temperature contrast
over the helium gradient region drives a steeper composition
gradient as per the phase diagram, a runaway effect ensues,
with Teff and Yatm plummeting as stronger stratifications are
realized. In the most extreme case shown here (R = 0.30),

the effective temperature decreases by 8 K over 108 yr in this

phase, to be compared to the roughly 1 K per 108 yr cooling
rate before the onset of helium rain. After roughly a thermal
time for the homogeneous interior, the gradient region feels
significant heating from below, quenching the runaway and
once again assuming a state of slow evolution in which the
surface cools by roughly 1 K per Gyr.
The effect of translating the phase curve in temperature
is summarized in Figures 7 and 8, which show evolutionary
tracks for two families of models, one with the phase diagram unmodified (solid curves) and one with a representative
offset of Tphase = 230 K for illustration purposes (dashed
curves). The crossing of the two families of curves in Figure 8 demonstrates the fundamental anticorrelation between
R and Tphase . For instance, an effective temperature of 124
K at 4.56 Gyr can be realized either by a model with relatively modest superadiabaticity using the unmodified phase
diagram, or by a model with a more extreme superadiabaticity and a delayed helium rain onset. However, the offset of
Tphase = 230 K delays the onset of helium rain by nearly
800 Myr and consequently leads to a more modest depletion of helium from the molecular envelope by the time the
model reaches the solar age. As demonstrated by Nettelmann
et al. (2015) and discussed in 2 above, a downward offset
of roughly this magnitude must be applied to the Lorenzen
et al. (2011) phase diagram to yield values of Yatm at the solar
age which are consistent with the Galileo entry probe measurement (Table 1). Lastly, we note that translating the phase
diagram to lower temperatures also leads to a more localized
helium gradient region, the Tphase = 230 K case yielding a
gradient roughly 1/3 the geometric thickness of the gradient
established in the Tphase = 0 case for each of the R values
considered here. In the following section, we leave Tphase
as a free parameter alongside Mc and R and systematically
estimate all three simultaneously by fitting Jupiters Teff and
Yatm along with its volumetric radius Rvol .
3.2. Bayesian Parameter Estimation
We estimate model parameters and their statistical uncertainties using Markov chain Monte Carlo (MCMC). In particular, given our nonlinear three-parameter model
{Mc , R , Tphase }

(25)

and fundamental Jupiter data (see Table 1)

D {Teff , Yatm , Rvol },

(26)

we calculate the posterior probability distribution from

Bayes theorem
P(|D) P(D|)P().

(27)

In the likelihood P(D|) we assume Gaussian errors for the

data D:
m 2
1h
(Teff Teff
)
ln P(D|) = ln(2T2eff ) +
2
2
Teff
+ ln(2Y2atm ) +
+ ln(2R2 vol ) +

m 2
(Yatm Yatm
)
2
Yatm
(R Rm )2 i
vol

R2 vol

(28)

where a superscript m denotes the model outcome at the solar

age; each sample from the posterior distribution thus entails

Helium Rain in Jupiter

0.30
=0.30
=0.25
=0.20
=0.10
adiabatic
homogeneous
R
R
R
R

0.25
0.20

Teff (K)

135
130
125
120

3.5 4.0 4.5 5.0

Age (Gyr)

1.0

1.5

2.0

2.5

1.0

3.5 Gyr

L (10 9 L )

0.15
2.0
1.5
1.0
0.5
0.0
2.5
2.0

1.5

2.0

2.5

1.0

3.8 Gyr

1.5

2.0

2.5

1.0

4.0 Gyr

1.5

2.0

2.5

4.56 Gyr

1.5

1.0

0.5
10

100

101 10

P (Mbar)

100

101 10

P (Mbar)

100

101 10

P (Mbar)

100

101

P (Mbar)

T (K)

R =0.00

R =0.25

R =0.20

R =0.30

0.8 1.0 1.2 1.4 1.6 1.8 2.0 0.8 1.0 1.2 1.4 1.6 1.8 2.0

P (Mbar)

4.55
4.45
4.35
4.25
4.15
4.05
3.95
3.85
3.75
3.65

P (Mbar)

F IG . 6. Evolution of the interior temperature profile for differentiating

Jupiter models (the same models as in Figure 5). Each panel plots a time
sequence of profiles for the model indicated. Color maps to model age.

a full evolutionary calculation. Samples are drawn from the

posterior distribution Equation 27 using the Python MCMC
implementation emcee (Foreman-Mackey et al. 2013). Our
prior probability distribution P() is chosen to be uniform
over Mc 0, all Tphase , and 0.1 R 1 following the criteria for linear instability (see 2.3, in particular Equation 15);
elsewhere it is zero.
The late evolution of Teff , Yatm and R for the a subset of

21000

R =0.30

0.25
0.20
0.15
0.10
adiabatic

20500

TCMB (K)

7500
7000
6500
6000
5500
5000
7500
7000
6500
6000
5500
5000

Age (Gyr)

T (K)

F IG . 5. Interior profiles for differentiating Jupiter models. Shown are the helium mass fraction Y (top row), the temperature gradient d ln T /d ln P
(middle row; see Equation 14) and local luminosity L (bottom row) as functions of pressure for 1.0 MJ models with Mc = 30 ME and Tphase = 0 K, for four
different values of the fractional superadiabaticity R : 0.20 (light orange), 0.25 (dark orange), 0.30 (red), and the adiabatic case R = 0 (blue). Thin black curves
show a model with no phase separation for reference. The four columns correspond to four points in model age as labeled in text in the lower panels, and as
indicated by the open circles in the inset in the center left panel. To emphasize detail, the first two rows show Y and over a different pressure scale than the
luminosity panels.

20000
19500
19000 Mc =30 ME
18500
4.25

Tphase =0
Tphase = 230 K

4.30

4.35

4.40

CMB (g cm

3)

4.45

F IG . 7. Evolution of the temperature and density at the core-mantle

boundary (the innermost grid point in our simulations) for six different values
of the superadiabaticity R of the temperature profile in the helium gradient
region. The two families of curves are for two representative temperature offsets applied to the H-He phase diagram: the solid tracks assume Tphase = 0;
the dashed tracks assume Tphase = 230 K. For each model with Tphase = 0
that reaches the solar age, that point is indicated with a filled circle.

the evolutionary sequences in the resulting chain are shown in

Figure 9, which color codes tracks by their value of R . For
the duration of the initial homogeneous phase, the differences
in Teff and R between tracks stem solely from the heavy element core mass Mc , since that parameter adjusts the mean

Mankovich, Fortney and Moore

130
125
0.25
120 Mc =30 ME R =0.30
3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0

Yatm

Teff (K)

135

140
135
130
125
120
0.28

a
b
c
c

Age (yr)

Teq (K)

Teff (K)

Yatm

Rvol (km)

4.56 109

109.0a

124.4 0.3b

0.234 0.005c

69911 6d

2.6
Hanel et al. (1981)
von Zahn et al. (1998)
Seidelmann et al. (2007)

density and hence the total radius of the planet. As described

in 3.1, the time at which phase separation sets in (and Yatm
first diverges from the protosolar value) is set by Tphase , and
the trend toward later phase separation onset with increasing
R is the result of the two parameters substantial covariance.
Timesteps spanning the solar age are typically 107 yr, and the
values of Teff , Yatm and R were interpolated linearly within
m
m
, and Rm . Proposed steps
those timesteps to obtain Teff
, Yatm
in which the model terminated before the solar age (as would
be expected for a model with simultaneously large R and
Tphase , for example) were rejected, since a calculation of the
likelihood (Equation 28) would not be possible.
The models that terminate before 5 Gyr did so either because (i) the model cooled to Teff = 120 K, at which point we
stop to avoid unnecessary computation time, or (ii) upon the
onset of phase separation, the luminosity inversion described
in the previous subsection (and evident in the R = 0.30 case
in Figure 5) grew to the degree that a negative luminosity was
realized in the interior, typically just outside the helium gradient region. This behavior follows from attempting to enforce
large values of R such that the the large luminosity generated
by deposition of helium into the metallic interior cannot be
communicated upward through the weakly turbulent doublediffusive layer. Nettelmann et al. (2015) identified the same
effect in their models with LDDC, noting that there exists a
minimum layer height such that the luminosity is still positive
throughout the interior. In our models this translates to an effective upper limit on the values of R attainable by models
with strictly positive luminosity profiles.
The outcome of our Bayesian parameter estimation is the
posterior probability distribution shown in Figure 10, wherein
each panel plots the full joint distribution marginalized over
the other two parameters. The medians of each distribution
are indicated, as are the central 68% confidence regions. Typ-

R (109 cm)

TABLE 1
G LOBAL J UPITER DATA

0.24
0.20

Age (Gyr)

F IG . 8. As in Figure 7, but showing the effective temperatures as a function of age. The marker shows Jupiters observed Teff at the solar age.

Teq =109 K

0.16
7.10
7.05
7.00
6.95
6.90
3.0

0.2
0
0.2
5
0.3
0
0.3
5
0.4
0
0.4
5

Tphase =0
Tphase = 230 K

Teff (K)

10

3.5

4.0

Age (Gyr)

4.5

5.0

F IG . 9. Late evolution of Jupiter models undergoing helium rain and

overstable double-diffusive convection, sampled using MCMC. Open circles
with error bars designate observations (Table 1); the error on Teff is smaller
than the marker. The color of an evolutionary a track encodes its R value
as specified by the colors of the bars in the histogram (inset, center panel),
which is a coarsely binned version of the marginalized posterior probability
distribution in Figure 10. Only a random subset of tracks are shown, and
more likely tracks are plotted on top of less likely ones.

ical models (as characterized by the medians) have massive

cores (Mc = 27.7 ME ), are strongly superadiabatic in the compositionally stratified region (R = 0.31), and have substantial
downward offsets to the phase diagram (Tphase = 235 K.)
The posterior distribution of Tphase is narrowly peaked at
these large negative offsets, and quite robustly rules out the
unmodified phase diagram. This result is driven by the requirement that the models homogeneous molecular envelope
retains enough helium to match the modest depletion measured by the Galileo entry probe (Table 1).
The single largest modeling uncertainty for the cooling
of the giant planets is associated with the equation of state.
Depending on the assumed EOS, Saumon & Guillot (2004)
found that cooling times for homogeneous Jupiter models
spanned 3.1 to 5.4 Gyr. (For comparison, our homogeneous
models with SCvH-I cool in 5 to 5.5 Gyr for realistic heavy
element masses and equilibrium temperatures; see Figures 3
and 4.) The range in cooling times obtained from applying
different equations of state owes mostly to the fact that for
different P(T ) relations, models with a given entropy possess
different total thermal energy content and hence take more or
less time to cool to Jupiters present-day luminosity.
As a means of exploring how our results would be affected
by the application of a different equation of state, we repeat
our full calculations with a modified atmospheric boundary
condition. Although the model atmosphere and equation of
state are not directly related, they are degenerate in that they
both dictate the overall timescale for the thermal evolution.
This can be made explicit by first integrating the energy equa-

11

that the cooling is extended to the age of the solar system.

140
135
130
125
120
0.28

Teq as free parameter

Teff (K)

Mc (ME)

Yatm

0.24

24
242
230
238
236
234
232
220
8

R (109 cm)

0.20

Tphase (K)
F IG . 10. Posterior probability distributions for the heavy element mass
Mc , density ratio R in the helium gradient region, and phase diagram offset
Tphase based on the evolutionary sequences shown in Figure 9. Each distribution is the full three-dimensional joint distribution marginalized over the
other two parameters. The vertical dashed line near each peak designates the
median, with the flanking vertical dashed lines enclosing the central 68% of
cumulative probability.

tion (equation 18a) over the mass of the planet to yield

Z M


ds
4
Lint =
T dm = 4R2 SB Teff
Teq4 .
dt
0

(29)

For the simplified example of a planet cooling through a sequence of isentropes, ds/dt is independent of m and the second equality can be integrated to yield the total cooling time
!
RM
Z cool
Z s0
T (m, s) dm
0

ds,
cool =
dt =
4 (s) T 4
4R2 (s)SB Teff
0
scool
eq
(30)
where s0 designates an arbitrary large starting entropy, scool
designates the planets current entropy, and other symbols
have their usual meanings. Equation 30, while emphatically
not how our evolutionary sequences are calculated, serves as a
heuristic tool to demonstrate that choosing a colder EOS (such
that the mean temperature along a given adiabat is lower) and
reducing the solar input Teq have a similar effect.
Our modification of the boundary condition as a proxy for
a different H-He EOS is motivated by the lack of other realistic EOS options presently available in MESA at the relevant densities and temperatures. Varying the parameter Teq
offers a simple means of producing a different total cooling
time, the relation between the two being illustrated in Figure 4. As an example, we find that a homogeneous, adiabatic
model with Mc = 30 ME and Teq reduced to 100 K cools to
Jupiters Teff in just 4.2 Gyr. Since very large superadiabaticities tend to reduce the cooling (see Figure 8), differentiating
models that satisfy the basic constraints of Table 1 in spite of
a cold boundary condition must have small values for R such

0.16
7.10
7.05
7.00
6.95
6.90
3.0

0.1
0
0.1
5
0.2
0
0.2
5
0.3
0
0.3
5

0.2
0.24
0. 6
0.328
0. 0
0.332
0.34
0. 6
0.438
0

27 Probability density
27.0
27 .2
27.4
27.6
28.8
28.0
28 .2
.4

Helium Rain in Jupiter

3.5

4.0

Age (Gyr)

4.5

5.0

F IG . 11. As in Figure 9, but including Teq as a additional free parameter

to mimic the effect of an EOS predicting a warmer or colder interior.

Figure 11 shows the late evolution of models computed

with this artificially free boundary condition, with a uniform
prior chosen for Teq ; the posterior distributions of Mc , R ,
Tphase and Teq are shown in Figure 12. In this case the extra freedom in the boundary condition leads to a wide variety of total times spent in the homogeneous phase of evolution, and consequently the other parameters Mc , R and
Tphase have markedly wider posterior probability distributions. Most likely models have equilibrium temperatures several K cooler than the measured value Teq = 109 K such that
the homogeneous cooling is more short-lived, and in the inhomogeneous evolution that follows, double-diffusive convection can proceed with temperature gradients closer to the adiabat and still satisfy the Teff constraint at the solar age. Indeed, the models display an overall preference for the lowest possible density ratios; the marginalized posterior probability density increases uniformly toward lower values of R
and peaks at the lower boundary imposed by the step-function
prior R > Rcrit = 0.1, which was imposed in light of the linear
criterion for the double-diffusive instability (see 2.4). As a
result, the most likely evolution during the inhomogeneous
phase is secular cooling of the envelope, with no dropoff
of the surface temperature or helium abundance over short
timescales. These models might be considered preferable to
the most likely models obtained in the three-parameter case
of Figures 9 and 10 because if Teff , Yatm , and Rvol were undergoing drastic changes on a 108 year timescale, then observing
Jupiter in its present state would be somewhat serendipitous.
Freeing Teq also allows much more modest corrections to
the phase diagram, the 68% credible interval now spanning
200 to 100 K. Importantly, despite the freedom in the overall cooling time, the unmodified phase diagram is still found

12

Mankovich, Fortney and Moore

96
98
10
100
102
104
106
118
0

0.1
0.10
0.12
0.14
0. 6
0.218
0. 0
0.222
0.24
6

.0
30

.0
29
.5

29

Probability density

Mc (ME)

25
200
150
100
0
50
0
50
10
0

28

.0
28
.5

Probability density

to be extremely unlikely with 95% of probability lying at offsets Tphase < 50 K.

The posterior distribution of Teq values in Figure 12 should
not be interpreted as a new determination of Jupiters Teq , as
that is a measured quantity. Rather, it contains information
about a most likely equation of state: it is significant that the
distribution excludes the measured value almost completely,
with only 0.4% of cumulative probability within the error of
the measurement (Table 1). If we suppose that our model contains the essential physics, this can be taken as evidence that
the real EOS for hydrogen and helium predicts substantially
colder interiors than does SCvH-I. For reference, a homogeneous evolutionary sequence computed with the median values Mc = 28.8 ME and Teq = 102.9 K from the distributions in
Figure 12 cools to Jupiters Teff in 4.46 Gyr (see Figure 4).

Tphase (K)

Teq (K)

F IG . 12. As in Figure 10, but including Teq as an additional free parameter

to mimic the effect of an EOS predicting a warmer or colder interior.

4. DISCUSSION

The framework developed here consists of a Python class

for creating instances of MESA work directories, modifying
MESA inlists, executing the evolution program, and processing its output, all as part of a single likelihood calculation
called by an emcee sampler. The method renders it straightforward to add additional parameters to the model or incorporate new or updated constraints in any quantity output by
MESA directly. It is readily adaptable to a host of different problems, most obviously non-adiabatic thermal evolution
models for Saturn, where the same fundamental physics operate. Beyond just the Jovian planets and H-He immiscibility,
our technique has broad applicability for deriving properties
of objects from giant planets to brown dwarfs and stars, e.g.,
retrieving the age and composition for an object with measured mass and radius. Performing these retrievals with a

code as mature as MESA means that our knowledge of stellar/planetary evolution is built in, including complexities such
as self-consistently determining mixing boundaries, modeling
double-diffusive transport processes, or calculating nuclear
energy generation rates with state of the art nuclear networks.
Thus in the example of retrieving an objects composition and
age from its measured mass and radius, meaningful inferences
can be made about not just the bulk composition, but the composition profile, and indeed the composition profiles possible origins and evolution. The Bayesian approach automatically provides meaningful error bars for model parameters,
and combining it with the open source MESA package offers
more flexibility than traditional grid-based isochrone fitting
because new parametersand indeed new physicscan be
added at will.
This work builds on that of Nettelmann et al. (2015) principally in two ways: first, it makes the weakest possible
assumption about the temperature gradient resulting from
double-diffusive convection in the deep interior, abandoning
the assumption of layered convection following the flux laws
derived by Wood et al. (2013) in favor of a generic model
wherein any temperature gradient can be attained as long as
it is consistent with the criterion for linearly overstable gravity waves. Second, performing the calculations in an MCMC
framework allows a probabilistic determination of all model
parameters simultaneously, and we find a multitude of models that satisfy the imposed constraints (Table 1). We demonstrated that SCvH-I predicts strongly superadiabatic temperature profiles in Jupiters helium gradient region, such that the
planets surface cools rapidly as most of the metallic hydrogen
interior heats up over time. Repeating the calculations with a
variable boundary condition to probe the effects of using a
different EOS, we found that more modest superadiabaticities
are preferred, although the distribution of allowed values is
still broad. We found in all cases that the unperturbed phase
diagram of Lorenzen et al. (2011) is highly unlikely.
That such a diversity of models meeting the imposed constraints were obtained in Figures 11 and 12 underscores the
severe uncertainties that persist in modelling the evolution of
giant planets. Admittedly, the present work does not exploit
all the available data. Most importantly, our models make
no use of of Jupiters gravitational harmonics or its axial moment of inertia, both of which constrain the interior density
profile. As discussed in 2, a calculation of oblateness and
the associated non-spherical components of the gravity field
is beyond our scope because the only EOS currently available
for modeling giant planets in MESA, SCvH-I, is limited to hydrogen and helium. All heavy elements are in an inert core
rather than partially distributed through the envelope, and as
such the density profiles in our models are somewhat unrealistic and are not suited for fitting to J2 or any higher-order
moments. Nonetheless, we view these models as complimentary to the detailed static models computing using more realistic equations of state (e.g., Hubbard & Militzer 2016) in
that we use a forward thermal evolution model to derive estimates for Jupiters deep superadiabatic temperature stratification and corrections to the H-He phase diagram, both of which
should be taken into consideration for improving static models of the Jovian planets. Our findings support the existing
body of evidence indicating that a realistic H-He equation of
state departs significantly from SCvH-I.
We thank Nadine Nettelmann for thoughtful comments

Helium Rain in Jupiter

that helped to significantly improve the manuscript. C.M.
acknowledges support from NASA Headquarters under the
NASA Earth and Space Science Fellowship Program (grant

13

NNX15AQ62H). Simulations for this research were carried

out on the UCSC supercomputer Hyades, which is supported
by National Science Foundation (grant AST-1229745) and the
University of California, Santa Cruz.

REFERENCES
Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315
Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M.
2007, ApJ, 661, 1094
Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP,
125, 306
Fortney, J. J., & Hubbard, W. B. 2003, Icarus, 164, 228
Fortney, J. J., Ikoma, M., Nettelmann, N., Guillot, T., & Marley, M. S. 2011,
ApJ, 729, 32
Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504
French, M., Becker, A., Lorenzen, W., Nettelmann, N., Bethkenhagen, M.,
Wicht, J., & Redmer, R. 2012, ApJS, 202, 5
Graboske, Jr., H. C., Olness, R. J., Pollack, J. B., & Grossman, A. S. 1975,
ApJ, 199, 265
Guillot, T., Chabrier, G., Gautier, D., & Morel, P. 1995, ApJ, 450, 463
Guillot, T., Stevenson, D. J., Hubbard, W. B., & Saumon, D. 2004, The
interior of Jupiter, ed. F. Bagenal, T. E. Dowling, & W. B. McKinnon,
3557
Hanel, R., et al. 1981, Science, 212, 192
Hubbard, W. B. 1968, ApJ, 152, 745
. 1969, ApJ, 155, 333
. 1970, ApJ, 162, 687
Hubbard, W. B., & Militzer, B. 2016, ApJ, 820, 80
Hubbard, W. B., & Smoluchowski, R. 1973, Space Sci. Rev., 14, 599
Kato, S. 1966, PASJ, 18, 374
Kippenhahn, R., & Weigert, A. 1990, Stellar Structure and Evolution, 192
Kurokawa, H., & Inutsuka, S.-i. 2015, ApJ, 815, 78
Landau, L. D., & Lifshitz, E. M. 1959, Fluid mechanics
Leconte, J., & Chabrier, G. 2013, Nature Geoscience, 6, 347
Ledoux, P. 1947, ApJ, 105, 305
Lorenzen, W., Holst, B., & Redmer, R. 2009, Physical Review Letters, 102,
115701
. 2011, Phys. Rev. B, 84, 235109
Militzer, B., & Hubbard, W. B. 2013, ApJ, 774, 148
Militzer, B., Hubbard, W. B., Vorberger, J., Tamblyn, I., & Bonev, S. A.
2008, ApJ, 688, L45

Mirouh, G. M., Garaud, P., Stellmach, S., Traxler, A. L., & Wood, T. S.
2012, ApJ, 750, 61
Morales, M. A., Hamel, S., Caspersen, K., & Schwegler, E. 2013,
Phys. Rev. B, 87, 174105
Morales, M. A., Schwegler, E., Ceperley, D., Pierleoni, C., Hamel, S., &
Caspersen, K. 2009, Proceedings of the National Academy of Science,
106, 1324
Nettelmann, N., Fortney, J. J., Moore, K., & Mankovich, C. 2015, MNRAS,
447, 3422
Paxton, B., Bildsten, L., Dotter, A., Herwig, F., Lesaffre, P., & Timmes, F.
2011, ApJS, 192, 3
Paxton, B., et al. 2013, ApJS, 208, 4
. 2015, ApJS, 220, 15
Pearl, J. C., & Conrath, B. J. 1991, J. Geophys. Res., 96, 18
Pstow, R., Nettelmann, N., Lorenzen, W., & Redmer, R. 2016, Icarus, 267,
323
Rosenblum, E., Garaud, P., Traxler, A., & Stellmach, S. 2011, ApJ, 731, 66
Salpeter, E. E. 1973, ApJ, 181, L83
Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713
Saumon, D., & Guillot, T. 2004, ApJ, 609, 1170
Schwarzschild, M., & Hrm, R. 1958, ApJ, 128, 348
Seidelmann, P. K., et al. 2007, Celestial Mechanics and Dynamical
Astronomy, 98, 155
Stevenson, D. J. 1975, Phys. Rev. B, 12, 3999
. 1979, MNRAS, 187, 129
Stevenson, D. J., & Salpeter, E. E. 1977, ApJS, 35, 239
Straus, D. M., Ashcroft, N. W., & Beck, H. 1977, Phys. Rev. B, 15, 1914
Sweigart, A. V., & Gross, P. G. 1974, ApJ, 190, 101
Turner, J. S. 1974, Annual Review of Fluid Mechanics, 6, 37
von Zahn, U., Hunten, D. M., & Lehmacher, G. 1998, J. Geophys. Res., 103,
22815
Walin, G. 1964, Tellus, 16, 389
Wood, T. S., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 157