UC Office of the President
Recent Work
Title
Saturn’sProbable Interior: An Exploration of Saturn’s Potential Interior Density Structures
Permalink
https://escholarship.org/uc/item/0vn7g57f
Journal
The Astrophysical Journal, 891(2)
ISSN
1538-4357
Authors
Movshovitz, Naor
Fortney, Jonathan J
Mankovich, Chris
et al.
Publication Date
2020-03-11
DOI
10.3847/1538-4357/ab71ff
Peer reviewed
eScholarship.org Powered by the California Digital Library
University of California
Saturn’s Probable Interior: An Exploration of Saturn’s Potential Interior Density
Structures
Naor Movshovitz
1
, Jonathan J. Fortney
1
, Chris Mankovich
1,4
, Daniel Thorngren
2,5
, and Ravit Helled
3
1
Department of Astronomy and Astrophysics, University of California, Santa Cruz, CA, USA; nmovshov@ucsc.edu
2
Department of Physics, University of California, Santa Cruz, CA, USA
3
Institute for Computational Science, Center for Theoretical Astrophysics & Cosmology, University of Zurich, Switzerland
4
Division of Geological and Planetary Sciences, Caltech, Pasadena, CA, US
5
Institute for research on exoplanets, Universite ‘de Montreal, Montreal, Canada
Received 2019 July 23; revised 2020 January 25; accepted 2020 January 29; published 2020 March 11
Abstract
The gravity field of a giant planet is typically our best window into its interior structure and composition. Through
comparison of a model planet’s calculated gravitational potential with the observed potential, inferences can be
made about interior quantities, including possible composition and the existence of a core. Necessarily, a host of
assumptions go into such calculations, making every inference about a giant planet’s structure strongly model
dependent. In this work, we present a more general picture by setting Saturn ’s gravity field, as measured during the
Cassini Grand Finale, as a likelihood function driving a Markov Chain Monte Carlo exploration of the possible
interior density profiles. The result is a posterior distribution of the interior structure that is not tied to assumed
composition, thermal state, or material equations of state. Constraints on interior structure derived in this Bayesian
framework are necessarily less informative, but are also less biased and more general. These empirical and
probabilistic constraints on the density structure are our main data product, which we archive for continued
analysis. We find that the outer half of Saturn’s radius is relatively well constrained, and we interpret our findings
as suggesting a significant metal enrichment, in line with atmospheric abundances from remote sensing. As
expected, the inner half of Saturn ’s radius is less well constrained by gravity, but we generally find solutions that
include a significant density enhancement, which can be interpreted as a core, although this core is often lower in
density and larger in radial extent than typically found by standard models. This is consistent with a dilute core
and/or composition gradients.
Unified Astronomy Thesaurus concepts: Saturn (1426); Planetary interior (1248); Planetary structure (1256); Planet
formation (1241)
1. Introduction
1.1. The Gravity Field as a Probe on the Interior
There are a number of fundamental questions that we would
like to understand about giant planets. Do they have a heavy-
element core? If so, what is its mass? Is it distinct from
the overlying H/He envelope, or partially mixed into it? Is the
H/He envelope enriched in heavy elements compared to the
Sun? Is the envelope fully convective and well mixed?
Unfortunately, the vast mass of a giant planet is completely
hidden from view, so that we must use indirect methods to try
to answer these questions. Most of our knowledge about the
interiors of giant planets comes from interpreting their gravity
fields, as recently reviewed for Saturn by Fortney et al. (2018).
Because the planets are fluid and rapidly rotating, they assume
an oblate shape and their gravitational potential differs from
that of a spherically symmetric body of the same mass. The
external gravitational potential
V
e
is a function of the colatitude
θ and distance r from the center of the planet, and is typically
written as an expansion in powers of
R
r
,
eq
where
R
e
q
is the
equatorial radius of the planet:
() ( ) ( ()
⎛
⎝
⎜
⎞
⎠
⎟
å
qq=- -
=
¥
Vr
GM
r
RrJP,1 cos.1
n
n
nne
1
eq
In Equation (1), P
n
are Legendre polynomials of degree n. The
coefficients J
n
(“the Js”) are measurable for many solar system
bodies by fitting a multiparameter orbit model to Doppler
residuals of spacecraft on close approach. For fluid planets in
hydrostatic equilibrium, where azimuthal and north–south
symmetry hold, only even-degree coefficients are nonzero.
When an interior model for a planet is created, the J
n
values
are calculated as integrals of the interior density over the
planetary volume:
()() ( ) ()
ò
rqt=- ¢ ¢ ¢ ¢rMR J r P dcos . 2
n
n
n
n
eq
These model J
n
can then be compared to measured ones. As is
well known, the different Js sample the density at different
depths (with J
2
probing deepest) but with significant overlap,
and with most of the weighting over the planet’s outer half in
radius. This point is illustrated in Figure 1.
The gravity field is a nonunique feature of the interior mass
distribution. In other words, different mass distributions can
lead to identical gravity signals. This complicates the process of
making inferences about the interior structure based only on the
external gravity field. In principle, one should explore a wide
range of possible interior structures, of possible
()
r
¢r
in
Equation (2), to see the full range of solutions that fit the
gravity field. Initially, researchers had focused on finding a
single, best-fit solution subject to a host of assumptions, chosen
for computational convenience and not necessarily following
reality. More recently, there have been efforts to explore an
expanded range of interior structures, usually by making
alternative assumptions about the prototypical planet.
The Astrophysical Journal, 891:109 (19pp), 2020 March 10 https://doi.org/10.3847/1538-4357/ab71ff
© 2020. The American Astronomical Society. All rights reserved.
1
The main contribution of the present work is the introduction
of a different approach to the task of inferring interior structure
from gravity and the application of this approach to Saturn. The
result is a suite of interior structure models of Saturn computed
with fewer assumptions and therefore showing a fuller range of
structures consistent with observation. We describe our method
in detail and compare it with previous work of similar spirit in
Sections 1.2.4 and 2.6.
1.2. Common Assumptions in Planetary Interior Models
There are typically at least three significant assumptions or
choices that modelers make when constructing interior models
of giant planets, thereby implicitly constraining the possible
inferences from these models.
1.2.1. The Planets Have Three Layers
Perhaps the most constraining assumption is the prototypical
picture of three layers, each well mixed enough to be considered
homogeneous. For Jupiter and Saturn, these are a helium-poor
outer envelope, a helium-rich inner envelope, and a heavy-
element, usually constant-density core. Investigators also adjust
the abundance of heavy elements in the He-rich and He-poor
layers, with little physical motivation other than it seems to
facilitate finding an acceptable match to the gravity field.
While a core–envelope structure is certainly a plausible one,
and indeed rooted in well-studied planet formation theories, the
assumption of compositionally homogeneous layers may well
be a significantly limiting oversimplification.
1.2.2. The Interior Pressure–Temperature Profile Is Isentropic
A typical assumption of interior modeling is that the
pressure–temperature profile is isentropic, lying on a single
(P, T) adiabat that is continued from a measured or inferred
temperature at 1bar. This second assumption is likely to be
true over some of the interior, but there are good reasons to
doubt that this holds throughout the interior.
Jupiter and Saturn have an atmospheric He depletion
compared to the Sun, and it has long been suggested that this
is due to He phase separation from liquid metallic hydrogen in
the deep interiors (Stevenson & Salpeter 1977; Fortney &
Hubbard 2003). There is likely a region with a He abundance
gradient starting between 1 and 2Mbar in both planets
(Nettelmann et al. 2015; Mankovich et al. 2016).Inmodels
that attempt to interpret the gravity field, if such a layer is
included at all, it is by interpolation between the outer and inner
homogeneous layers (e.g., Wahl et al. 2017; Militzer et al. 2019),
but this interpolation is unlikely to capture fully the effects of
composition gradients. Composition gradients can inhibit large-
scale convection (Ledoux 1947), implying that heat is
transported via layered convection or radiation/conduction. This
leads to higher internal temperatures that in return allow higher
heavy-element enrichment at a given density–pressure. Indeed,
nonadiabatic structures have been recently suggested for all outer
planets in the solar system (e.g., Leconte & Chabrier 2012;
Vazan et al. 2016; Podolak et al. 2019).
1.2.3. The Inferred Composition Relies on Equation of State
Calculations
As the field progresses, the equations of state (EOSs) used
for modeling giant planet interiors become a better representa-
tion of reality. Nevertheless, the EOSs for all relevant materials
and mixtures are not perfectly known. Simulations from first
principles of hydrogen, helium, and their mixtures over the
conditions relevant for giant planets have been carried out and
partially validated against experimental data (Nettelmann et al.
2008; Militzer & Hubbard 2013). These EOSs for hydrogen
show good agreement with data up to
~1.5 Mbar
(e.g., Militzer
et al. 2016). However, the pressure at the bottom of Saturn’s
H/He envelope is about 10Mbar and for Jupiter it is about
40Mbar, well beyond the realm of experiment. Recent
structure models used EOSs for hydrogen and helium based
on density functional theory (DFT) simulations (Nettelmann
et al. 2008; Miguel et al. 2016; Militzer et al. 2016). Until
recently, different EOSs led to different inferred compositions
for Jupiter due to different approaches to calculating the
entropy. Today, there is good agreement between state-of-the-
art EOSs (Nettelmann 2017), but it should be kept in mind that
DFT also suffers from approximations (Mazzola et al. 2018),
and there remains an uncertainty of
~2%
in the hydrogen EOS,
which increases significantly when it comes to predicting
hydrogen–helium demixing (Morales et al. 2009).
The heavy elements must also be represented by an EOS
(typically for water or silicates), which introduces another
source of uncertainty. Therefore, the range of possible
compositions and internal structures from such interior models
Figure 1. Contribution functions of the gravitational harmonics J
2
(blue solid), J
4
(red dashed), and J
6
(yellow dotted) for a typical, three-layer Saturn model. The
contribution “density” (
()µJrdr
) is plotted in the left panel and the cumulative contribution in the right panel. The horizontal line intersects the curves at a depth
where the corresponding J reaches 90% of its final value.
2
The Astrophysical Journal, 891:109 (19pp), 2020 March 10 Movshovitz et al.
cannot be taken to be the true range of allowed values, even if
the parameter space of possible EOSs, H/He/Z mixing ratios,
and outer/inner envelope transition pressures were thoroughly
explored.
1.2.4. Appreciating the Complexities
The drawbacks of the assumptions discussed above have
long been known and the reality that giant planets are surely
more complicated than the traditional modeling framework
allows for is generally accepted (e.g., Stevenson 1985). More
recent investigations are attempting to allow for a more
complex structure. Interior composition gradients due to
remnants of formation ( Leconte & Chabrier 2012; Helled &
Stevenson 2017), core dredge-up (Militzer et al. 2016),
convective mixing of primordial composition gradients (Vazan
et al. 2016, 2018), and He sedimentation (Nettelmann et al.
2015; Mankovich et al. 2016) have been considered and were
found to lead to different structures. Additionally, some
investigators have begun using what may be referred to as
“empirical” models. In this context, an empirical model is one
that is focused on the more direct connection between gravity
and density (e.g., Helled et al. 2009) or gravity and equilibrium
shape (Helled et al. 2015), without invoking the compositional
and thermodynamical origin of these structures.
The work we present here is in the spirit of empirical models.
We explore systematically, in a Bayesian inference framework,
the possible density profiles of Saturn. We limit our assumptions
as much as possible, in order to find the widest range of interior
structures with their probability distribution based on their
gravitational potential matching the observed field.
2. Composition-independent Interior Density Calculation
The premise of removing some assumptions and deriving
composition-free interior density profiles (sometimes referred
to as empirical models) is simple and in fact has been pursued
in previous works (e.g., Marley et al. 1995; Podolak et al. 2000;
Helled et al. 2009). (We discuss similarities and differences
with these works in Section 2.6 below.) The only information
that is needed to calculate a gravity field is the density
everywhere inside the planet,
()
r
r
, and so this is the only
quantity we will directly vary. In fact, hydrostatic equilibrium
produces level surfaces—closed surfaces of constant density,
pressure, and potential—and therefore a one-dimensional
description of the mass distribution is sufficient: we can use
() ()
r
r=r s
, where s is the volumetric mean radius of the
unique level surface of density ρ.
All other properties of the planet will be inferences, rather
than input parameters. Because there is unavoidable uncertainty
associated with the measurement of the gravity field (and also
with its theoretical calculation from interior models; see
Section 2.4), this means that there must be a continuous
distribution of possible density profiles that fit the gravity
solution, and we must base our inferences on the entire
distribution. In practice, because we can only ever consider a
finite number of solutions, this means that we must base our
inferences on a random sample from this unknown distribution
of allowed solutions.
In this section, we describe the process of obtaining this
random sample, as applied to Saturn. For the most part, the
same process would apply equally well to the other giant
planets. We mention in places modifications that may be
needed if the same method is to be applied to Jupiter, Uranus,
or Neptune.
2.1. Overview
Formally, the distribution we are after is the posterior
probability
(∣ )r
Jp
, the probability that the planet’s interior
density follows
()
r
r= s
given that the gravity coefficients
were measured as
J
. This consists of several subtasks. First,
we must find a suitable parameterization of
(
)r
s
. This
parameterization should be able to represent all the physically
reasonable
(
)r
s
curves without undue loss of generality, but
this is not particularly difficult. It is also necessary that the
range and behavior of the numeric values of all parameters are
such that they can be efficiently sampled, e.g., with a Markov
Chain Monte Carlo (MCMC ) algorithm. This is easier said than
done, and the best parameterization may be different for
different planets.
For Saturn, we find that a piecewise-quadratic function of
density as a function of normalized radius works best for the
bulk of the planet, with a quartic (degree 4) polynomial
required to represent the uppermost region (for
P 2GP
a
).
We describe this parameterization in detail in Section 2.2. Note
that this is one place where modifications might be needed
before applying the same procedure to Jupiter or the ice giants.
To drive the sampling algorithm, we need a way to evaluate
the relative likelihood of two model planets, and we do this by
comparing how well they match Saturn’s observed mass and
gravity field. The details of this calculation are given in
Section 2.3.
The likelihood calculation requires that we know the
equilibrium shape and gravity field of a given density profile
to sufficient accuracy. Note that in Equation (2) the integrand is
known but the integration bounds are unknown. We need to
first determine the planet’s equilibrium shape. The shape is
determined by a balance between the centrifugal acceleration of
the rotating planet and the gravitational acceleration. This is
therefore a circular problem, requiring an iterative calculation
to converge to a self-consistent solution.
We use an implementation of fourth-order Theory of Figures
(ToF) using the coefficients given in Nettelmann ( 2017) and
employ optimization techniques that allow us to solve the
hydrostatic equilibrium state to desired precision very quickly.
The details are given below in the Section 2.4.
The emphasis on speed is necessary, as the next subtask is to
employ a suitable MCMC algorithm to draw a large sample
of possible
r
. There is no generally agreed-upon method of
predicting the number of sampling steps required for
convergence.
6
By experimentation, we find that our Saturn
parameterization requires tens of thousands of steps to become
independent of its seed state and has a long autocorrelation
time, requiring a large number of steps following convergence
to obtain the desired effective sample size. Producing a valid
sample required the computation of about 10 million model
planets in total. We give the details of our sampling method
and convergence tests in Section 2.5 and Appendix C.
The last step is calculating some derived physical quantities
of interest, based on the obtained
r
sample. Given the gravity
field, the pressure on each level surface can be computed from
the hydrostatic equilibrium equation. And with knowledge of
the pressure and density at each level, we may begin to estimate
6
Or even of being sure that convergence was reached.
3
The Astrophysical Journal, 891:109 (19pp), 2020 March 10 Movshovitz et al.
other quantities of interest, e.g.,the helium fraction, the heavy-
element content, etc. These quantities are not determined
directly by the gravity field but can be inferred, with additional
assumptions. We discuss the results of this analysis, as applied
to Saturn, in Section 3.
2.2. Parameterization of
(
)r
s
Our goal is to sample from a space of
(
)r
s
curves that is as
general as possible, making a minimum of assumptions about
(
)r
s
while still restricting the sample to physically meaningful
density profiles and, importantly, keeping the number of free
parameters small, for sampling efficiency. These competing
requirements are not easy to satisfy, and it may be that the best
parameterization depends on the planet being studied as well as
on the available sampling algorithms and computing resources.
When looking for a good parameterization of
(
)r
s
, we were
guided by previously published work on Saturn’s interior.
Traditionally derived models are less general than we would
like, but they are physically sound. Examining them exposes
the major features expected of a
(
)r
s
curve representing
Saturn’s interior. Figure 2 shows the density profiles of several
Saturn models recently published by Mankovich et al. (2019,
hereafter M19). These models assume a three-layer structure
for Saturn along the lines of what was considered by
Nettelmann et al. (2013). They consist of a homogeneous
outer envelope with helium mass fraction
=
Y
Y
1
and water
mass fraction
=ZZ
1
, a homogeneous inner envelope with
=
Y
Y
2
and Z= Z
2
, and finally a central core with Z=1.
These models assume an additive-volume mixture of hydrogen,
helium, and water as described by the Saumon et al. (1995) and
French et al. (2009) EOSs, and are assumed to have adiabatic
temperature profiles throughout the envelope with an iso-
thermal core.
The general feature is a monotonic and piecewise-smooth
function in three segments. This is not surprising, as these
models were all derived with the assumption of three layers of
homogeneous composition, commonly thought of as an upper
envelope, lower envelope, and core. While we do not wish to
make such a strong assumption, we find it necessary to make
the much weaker assumption that
(
)r
s
is a monotonic,
piecewise-smooth function, with no more than (but possibly
fewer than!) two density discontinuities. Further, between
discontinuities, the density appears to follow very smooth
curves, suggesting that it may be well approximated by
a quadratic function of s/R
m
for each segment, where
=
R
58, 232 km
m
is Saturn’s volumetric mean 1 bar radius
(Lindal et al. 1985). By experimentation, we find no advantage
in using higher-order polynomials to approximate any of the
main segments.
This piecewise-continuous model should not be confused
with the traditional three-layer one. The assumption of density
being piecewise continuous is much less strict than that of
composition being piecewise constant, even if they lead to
visually similar plots. Nevertheless, it would be even better to
allow more discontinuities or, better yet, a variable number of
them. While this may seem like a relatively straightforward
generalization, it would in fact greatly increase the computa-
tional cost of sampling the parameter space. To understand
why, consider that each additional discontinuity in
(
)r
s
not
only introduces four additional parameters (the three para-
meters required to describe the quadratic plus the location of
the additional break point ), these parameters will also be highly
correlated with the rest. As it turns out, this correlation is
already evident with just two discontinuities. Informally, each
of the two density “jumps” can substitute for the other in the
large subset of models where only a single pronounced
discontinuity appears. This evident “redundancy” is by no
means proof that there cannot be more than two sharp density
jumps in Saturn’s interior. But it helps us accept, at least
temporarily, a compromise between maximum generality and
minimum CPU hours.
When we examine more closely the very top of the density
curves in Figure 2 we find that the uppermost part of the
envelope (by radius,
=
r
rR0.94
a m
) does not follow the
same quadratic as the rest of the upper envelope. Instead, it is
more similar to a quartic polynomial. This is demonstrated
visually for one density profile in the inset of Figure 2 and in
more detail in Appendix A. In this low-pressure region, the
physical models are based on well-tested EOSs of H and He,
and the assumption of an adiabatic gradient is appropriate, so
we would be well advised to constrain our profiles to make use
of this information. In Appendix A we explain how we derive a
one-parameter family of quartic functions that keeps us
grounded to realistic density values in the region above r
a
,
while still allowing variation by letting the value of
()
r
r==sr
a
a
be sampled.
It is important to note that to date, all EOS-based models of
Saturn find solutions consistent with the measured gravity field
that predict a concentration of heavy elements in the envelope
of at most a few times the protosolar value (e.g., Helled &
Guillot 2013; Nettelmann et al. 2013; Militzer et al. 2019),
while Saturn’s atmospheric spectra indicate a higher value,
perhaps as high as 10 times the protosolar metallicity (Atreya
et al. 2016). In principle, atmospheric enrichment might not
represent the bulk composition of the outer envelope, as was
Figure 2. Three representative Saturn density profiles from M19. These
profiles were derived using the standard, three-layer assumption, and thus
represent only a subset of possible profiles. On the other hand, they are known
to be in strict agreement (by construction) with theoretical EOSs throughout the
interior. The inset shows a zoomed-in view of the top part of the envelope.
The red solid line is the same curve as in the full-scale figure; the black dotted
line is a quadratic fit, a good approximation of the upper envelope overall, and
the blue dashed line is a quartic fit to the segment
>sR 0.9
4
m
, a much better
fit there (Appendix A).
4
The Astrophysical Journal, 891:109 (19pp), 2020 March 10 Movshovitz et al.
