Models of Saturn's Interior Constructed with Accelerated Concentric Maclaurin Spheroid Method

TL;DR

This paper develops an accelerated Concentric Maclaurin Spheroid method to model Saturn's interior, matching gravity data and predicting a dense core and rotation period consistent with observations.

Contribution

It introduces an efficient acceleration technique for the CMS method, enabling precise modeling of Saturn's interior and gravity harmonics with reduced computational cost.

Findings

Saturn's core mass is 15-18 Earth masses.

Heavy elements in the envelope amount to 1.5-5 Earth masses.

Predicted rotation period is approximately 10:33:34 hours.

Abstract

The Cassini spacecraft's Grand Finale orbits provided a unique opportunity to probe Saturn's gravity field and interior structure. Doppler measurements yielded unexpectedly large values for the gravity harmonics J_6, J_8, and J_10 that cannot be matched with planetary interior models that assume uniform rotation. Instead we present a suite of models that assume the planet's interior rotates on cylinders, which allows us to match all the observed even gravity harmonics. For every interior model, the gravity field is calculated self-consistently with high precision using the Concentric Maclaurin Spheroid (CMS) method. We present an acceleration technique for this method, which drastically reduces the computational cost, allows us to efficiently optimize model parameters, map out allowed parameter regions with Monte Carlo sampling, and increases the precision of the calculated J_2n gravity harmonics to match the error bars of the observations, which would be difficult without acceleration. Based on our models, Saturn is predicted to have a dense central core of 15-18 Earth masses and an additional 1.5-5 Earth masses of heavy elements in the envelope. Finally, we vary the rotation period in the planet's deep interior and determine the resulting oblateness, which we compare with the value from radio occultation measurements by the Voyager spacecraft. We predict a rotation period of 10:33:34 h +- 55s, which is in agreement with recent estimates derived from ring seismology.

Figures

Tables

Table 1. Table 1: Original and optimized λ 𝜆 \lambda grids for an interior model with 2049 sphoeroids and r C subscript 𝑟 𝐶 r_{C} =0.231. The entire table will be published online.

Spheroid index $i$	Original ${\lambda }_i$	Optimized ${\lambda }_i$
0	1.0000000000	1.0000000000
1	0.9999958561	0.9999966866
2	0.9999912702	0.9999933632
3	0.9999861950	0.9999900322
$\dots$	$\dots$	$\dots$
2047	0.2316114866	0.2338981435
2048	0.2310000000	0.2310000000

Table 2. Table 2: Gravity coefficients for the polytrope index 1 planet models derived with the accelerated CMS (this work) and CLC  (Wisdom & Hubbard, 2016 ) methods.

Gravity coefficient	CMS	CLC
${10}^2\times J_2$	1.3988511	1.398851090
${10}^4\times J_4$	$-$5.318281	$-$5.318281001093
${10}^5\times J_6$	3.0118324	3.011832290534
${10}^6\times J_8$	$-$2.1321158	$-$2.132115710725
${10}^7\times J_{10}$	1.74067138	1.740671195866
${10}^8\times J_{12}$	$-$1.56821961	$-$1.568219505563
${10}^9\times J_{14}$	1.51809944	1.518099226841
${10}^{10}\times J_{16}$	$-$1.5519853	$-$1.551985081630
${10}^{11}\times J_{18}$	1.6559259	1.655925984019
${10}^{12}\times J_{20}$	$-$1.8285783	$-$1.828574676495

Table 3. Table 3: Correction to the computed gravitional moments due to upper atmosphere that extends from 100 to 0.9 mbar.

$n$	${10}^6\times \Delta J_n$
2	$+6.3\times {10}^{-3}$
4	$-1.1\times {10}^{-3}$
6	$+2.5\times {10}^{-4}$
8	$-7.6\times {10}^{-5}$
10	$+2.2\times {10}^{-5}$
12	$-1.7\times {10}^{-6}$
14	$-2.3\times {10}^{-6}$
16	$+4.3\times {10}^{-7}$

Table 4. Table 4: Parameters and constraints in our Saturn models

$S_{\mathrm{mol}}$	Entropy H-He gas throughout the molecular layer.
	Constrained to match Saturn’s 1 bar temperature of 142.6 K (Lindal et al., 1981)
$Y_{\mathrm{mol}}$	Helium mass fraction in the molecular layer. Constraint: $Y_{\mathrm{mol}}\le Y_{\mathrm{solar}}=0.2741$ (Lodders, 2003).
$Z_{\mathrm{mol}}$	Mass fraction of heavy Z elements in the molecular layer.
$S_{\mathrm{met}}$	Entropy H-He gas throughout the metallic layer. $S_{\mathrm{met}}\ge S_{\mathrm{mol}}$
$Y_{\mathrm{met}}$	Helium mass fraction in the metallic layer. Adjusted as function of $Y_{\mathrm{mol}}$
	to keep the overall composition of the envelope at $Y_{\mathrm{solar}}$.
$Z_{\mathrm{met}}$	Mass fraction of heavy Z elements in the metallic layer. Constraint: $Z_{\mathrm{met}}\ge Z_{\mathrm{mol}}$.
$P_1$	Starting pressure of the helium rain layer (high pressure end of molecular layer)
$P_2$	Ending pressure of the helium rain layer (low pressure beginning of metallic layer)
Core mass	We assumed a compact core composed of heavy elements
	with a sharp boundary to the metallic layer, with an equatorial radius, $r_C$.
$\omega \left(l_k\right)$	Angular frequency for cylinder of radius, $l_k$. Constraint: .

Table 5. Table 5: Gravity coefficients predicted without acceleration scheme for different number of layers, N L subscript 𝑁 𝐿 N_{L} . A representative Saturn interior model with uniform rotation was used for this convergence analysis.

$N_L$	$J_2$	$J_4$	$J_6$	$J_8$	$J_{10}$	$J_{12}$	$J_{14}$	$J_{16}$
16	21058.747614	-1308.699233	119.572180	-13.513972	1.750222	-0.249135	0.037999	-0.006108
32	17740.199991	-1047.173049	92.601757	-10.209559	1.294184	-0.180542	0.027002	-0.004258
64	16789.968480	-978.984988	85.931473	-9.416858	1.186999	-0.164690	0.024500	-0.003843
128	16552.657945	-962.339808	84.321088	-9.226705	1.161402	-0.160918	0.023907	-0.003745
256	16493.667469	-958.249449	83.927467	-9.180331	1.155164	-0.159999	0.023762	-0.003721
512	16478.115263	-957.156333	83.821482	-9.167812	1.153482	-0.159752	0.023724	-0.003715
1024	16474.382140	-956.897102	83.796519	-9.164871	1.153087	-0.159694	0.023714	-0.003713
2048	16473.437419	-956.831274	83.790168	-9.164122	1.152986	-0.159679	0.023712	-0.003713
4096	16473.201201	-956.814816	83.788580	-9.163935	1.152961	-0.159675	0.023712	-0.003713
8192	16473.142127	-956.810701	83.788183	-9.163888	1.152955	-0.159674	0.023711	-0.003713

Table 6. Table 6: Gravity coefficients for the Saturn interior model in Tab.  5 predicted with acceleration factor n int = 16 subscript 𝑛 int 16 n_{\rm int}=16 . The first column denotes the number of CMS layers that were treated explicitly and the second specifies the total layer number. The last row contains the extrapolated values for N L tot → ∞ → superscript subscript 𝑁 𝐿 tot N_{L}^{\rm tot}\to\infty .

$N_L^{\mathrm{tot}}/n_{\mathrm{int}}$	$N_L^{\mathrm{tot}}$	$J_2$	$J_4$	$J_6$	$J_8$	$J_{10}$	$J_{12}$	$J_{14}$	$J_{16}$
8	128	16538.897354	-961.501994	84.237944	-9.215607	1.159703	-0.160646	0.023863	-0.003738
16	256	16498.268862	-958.176106	83.922106	-9.180301	1.155204	-0.160008	0.023764	-0.003721
32	512	16478.086394	-957.150893	83.821249	-9.167769	1.153475	-0.159751	0.023723	-0.003715
64	1024	16474.446344	-956.898550	83.796655	-9.164888	1.153089	-0.159694	0.023714	-0.003713
128	2048	16473.448194	-956.831549	83.790193	-9.164125	1.152986	-0.159679	0.023712	-0.003713
256	4096	16473.202955	-956.814864	83.788584	-9.163936	1.152961	-0.159675	0.023712	-0.003713
512	8192	16473.142447	-956.810710	83.788183	-9.163888	1.152955	-0.159674	0.023711	-0.003713
1024	16384	16473.127416	-956.809674	83.788084	-9.163877	1.152953	-0.159674	0.023711	-0.003713
2048	32768	16473.123665	-956.809415	83.788059	-9.163874	1.152953	-0.159674	0.023711	-0.003713
$\infty$	$\infty$	16473.122342	-956.809322	83.788050	-9.163873	1.152952	-0.159674	0.023711	-0.003713

Equations

$$\nabla P=\rho \nabla U.$$

$$U=V+Q,$$

$$q_{\mathrm{rot}}=\frac{{\omega }^2{a}^3}{GM},$$

$$P_{\mathrm{pu}}\equiv$$

$${\rho }_{\mathrm{pu}}\equiv$$

$$U_{\mathrm{pu}}\equiv$$

$${\zeta }_i(\mu )\equiv \frac{r_i(\mu )}{r_i(0)}\le 1$$

$${\delta }_i=\{\begin{array}{ll}{\rho }_i-{\rho }_{i-1},& i>0\\ {\rho }_0,& i=0.\end{array}$$

$$f_{im}({\zeta }_{im})\equiv U_i({\zeta }_{im},{\mu }_m)-U_i(1,0),$$

$${\zeta }_{im}^{(\mathrm{new})}={\zeta }_{im}-\frac{f_{im}({\zeta }_{im})}{f_{im}^{\prime }({\zeta }_{im})}.$$

$$P_i^{(\mathrm{new})}=P_{i-1}^{(\mathrm{new})}+{\rho }_{i-1}(U_i-U_{i-1}).$$

$${\rho }_i^{(\mathrm{new})}=\rho ((P_{i+1}+P_i)/2),$$

$$V(\mathbf{r})=G\int d^3{\mathbf{r}}^{\prime }\frac{\rho ({\mathbf{r}}^{\prime })}{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}.$$

$$V(r,\mu )$$

$$=$$

$$J_n=-\frac{2\pi }{M{a}^n}\underset{-1}{\overset{+1}{\int }}d\mu \underset{0}{\overset{r_{\max }(\mu )}{\int }}dr{r}^{n+2}{P}_n(\mu )\rho (r,\mu ).$$

$${\tilde{J}}_{i,n}\equiv \frac{J_{i,n}}{{\lambda }_i^n}\mathrm{and}{\tilde{J}}_{i,n}^{\prime }\equiv J_{i,n}^{\prime }{\lambda }_i^{(n+1)}.$$

$${\tilde{J}}_{i,n}=-\frac{1}{n+3}\frac{2\pi }{M}{\delta }_i{\lambda }_i^3\underset{-1}{\overset{+1}{\int }}d\mu P_n(\mu ){\zeta }_i(\mu {)}^{(n+3)}$$

$${\tilde{J}}_{i,n}^{\prime }=-\frac{1}{2-n}\frac{2\pi }{M}{\delta }_i{\lambda }_i^3\underset{-1}{\overset{+1}{\int }}d\mu P_n(\mu ){\zeta }_i(\mu {)}^{(2-n)}$$

$${\tilde{J}}_{i,n}^{\prime }=-\frac{2\pi }{M}{\delta }_i{\lambda }_i^3\underset{-1}{\overset{+1}{\int }}d\mu P_n(\mu )\log ({\zeta }_i)$$

$$J_{i,0}^{\prime \prime }=\frac{2\pi {\delta }_i{a}^3}{3M},$$

$$M=\frac{2\pi }{3}\sum _{i=0}^{N-1}{\delta }_i{\lambda }_i^3\underset{-1}{\overset{+1}{\int }}d\mu {\zeta }_i(\mu {)}^3$$

$$V_i({\zeta }_i,\mu )$$

$$V_{i=0}(1,0)=-\sum _{n=0}^{\infty }{P}_n(0)J_n,$$

$$J_n=\sum _{i=0}^{N-1}{\lambda }_i^n{\tilde{J}}_{i,n}$$

$$Q={\int }_0^{l}d{l}^{\prime }{l}^{\prime }\omega (l^{\prime }{)}^2$$

$$\omega (l)={\omega }_0+c_2{l}^2+c_4{l}^4+c_6{l}^6+c_8{l}^8+\dots ,$$

$$\log \left|\Delta J_n(N_L)\right|\equiv \log \left|J_n(N_L)-J_n(\infty )\right|=A-B\log [N_L]$$

$${\chi }_J^2=\sum _{i=1}^5{\left[\frac{J_{2i}^{\mathrm{observed}}-J_{2i}^{\mathrm{model}}}{\delta J_{2i}^{\mathrm{observed}}}\right]}^2,$$

$${\chi }_{\mathrm{H}-\mathrm{He}}^2=\sum _{i=1}^2{C}_P\left|\frac{P_i^{\ast }-P_i}{P_i}\right|+C_T\left|\frac{T_i^{\ast }-T_i}{T_i}\right|,$$

$${\chi }_{\mathrm{clouds}}^2=C_{\mathrm{clouds}}\sum _{k\mathrm{with}{l}_k>0.8}\left|{\omega }_{\mathrm{observed}}(l_k)-{\omega }_{\mathrm{model}}(l_k)\right|,$$

$${\chi }_{\mathrm{curvature}}^2=C_{\mathrm{curvature}}\sum _{\mathrm{all}k}{\left[{\omega }_{\mathrm{model}}^{\prime \prime }(l_k)\right]}^2,$$

Full text

I Introduction

Although Saturn’s deep interior was not a primary target of the Cassini spacecraft’s 13-year mission monitoring the Saturnian system, the final phase of the mission provided unprecedentedly precise measurements of the planet’s gravitational field (Iess et al., 2019). This phase, from April 23 to Sept. 15, 2017, culminated in 22 Grand Finale orbits, during which the Cassini spacecraft dived between the planet and its innermost ring. These measurements were contemporaneous with the ongoing Juno mission, which is providing analogous measurements for Jupiter (Folkner et al., 2017). As a result of both studies, the measured gravity fields are far more precise than ever before, warranting a closer look at the theory and numerical techniques linking the observed gravity to the interior density structure of the planet. Here we present models of Saturn’s interior structure and interior rotation rate, matched to the Cassini measurements, along with an acceleration technique for the Concentric Maclaurin Spheroid (CMS) method (Hubbard, 2013) for calculating a self-consistent shape and gravity field.

Prior to Cassini’s Grand Finale, the best determination of Saturn’s gravity field was from earlier flyby missions and from perturbations of the orbits of Saturn’s natural satellites in combination with the orbit of Cassini itself (Jacobson et al., 2006). However, this yielded significant measurements of only the first three even zonal harmonics of the field, $J_{2}$, $J_{4}$ and $J_{6}$. By contrast, X-band Doppler measurements during five of the 22 Grand Finale orbits produce a fit with significant determination of even zonal harmonics up to $J_{12}$, as well as odd zonal harmonics $J_{3}$ and $J_{5}$ (Iess et al., 2019).

The distribution of mass within a planet depends on the equation of state of hydrogen-helium mixtures at high pressures (Militzer & Hubbard, 2013), as well as the radial distribution of heavier elements (Soubiran & Militzer, 2016). The interior density distribution influences the observed structure of the gravity field through deviations from spherical symmetry arising from rotation and tides. Thus, the measured field can place constraints, albeit non-uniquely, on the internal structure of the planet. For the rapidly rotating Jovian planets, such terms are primarily determined by the balance between centrifugal and gravitational forces. In the absence of internal dynamics, the density distribution and resulting gravity field are axisymmetric and north-south symmetric, implying that only even zonal harmonics $J_{2n}$ contribute to the gravitational potential.

If a planet in hydrostatic equilibrium rotates uniformly like a solid body, the magnitudes of even zonal harmonics decay as $|J_{2n}|\sim q_{\rm rot}^{n}$, where $q_{\rm rot}$ is the ratio of the centrifugal and gravity accelerations at the equator. The $J_{2n}$ of Jupiter measured by Juno spacecraft are broadly consistent with this relationship (Folkner et al., 2017), meaning that it is possible to find models with a uniform rotation rate that match the observed $J_{2n}$, at least in the absence of other constraints, from the hydrogen-helium equation of state and atmospheric composition. However, Fig. 1 illustrates how the observed even moments $J_{8}$ and higher for Saturn deviate significantly from the expected relationship. Iess et al. (2019) demonstrated that these observations cannot be reproduced with models that assume uniform rotation, and that deep differential rotation (Hubbard, 1982) is required instead. In this paper we expand upon the interpretation of Iess et al. (2019) and introduce new analytical tools for high-precision gravity modeling.

I.1 Differential Rotation

Over many years prior to and including the duration of the Cassini mission, optical tracking of clouds has revealed large-scale zonal wind currents with respect to the average Saturn atmosphere, in particular a pronounced eastward jet centered on the equator (Sanchez-Lavega et al., 2000; García-Melendo et al., 2011). However, prior to the gravity measurements discussed here, the data were insufficient to constrain the depth of such zonal flows, and their effects were not considered in previous modeling studies of Saturn’s interior (Helled & Guillot, 2013; Nettelmann et al., 2013). With the Grand Finale gravity data, it becomes possible to test a model in which the cloud-level zonal wind belts are mapped onto cylinders that extend to great depths. If the zonal-wind velocity profile continues to depths of many scale heights, it will affect the observed gravity field in two ways. First, it modifies the axisymmetric gravitational field, and thus changes the even $J_{2n}$ from the values expected for a uniformly rotating body with identical internal structure (Hubbard, 1982). Second, to the extent that the velocity profile is not north-south symmetric, there arises a corresponding asymmetry in the gravity field, manifesting itself in non-zero odd $J_{n}$ (Kaspi, 2013). The values of $J_{3}$ and $J_{5}$ reported by Iess et al. (2019) thus exhibit the north-south asymmetric component of the differential rotation.

There are currently two basic methods for incorporating differential rotation into gravity models. The first is to approximate the wind profile as rotation on cylinders, which can be described using potential theory (Hubbard, 1982) and can therefore be integrated directly into the potential used in the CMS simulation (Wisdom & Hubbard, 2016). This method has the benefit of being fully self-consistent; the dynamic contribution to the potential modifies the shape of the equipotential surfaces, which feeds back into the calculated gravitational field. The downside is that the wind profile must be constant on cylindrical surfaces and thus cannot decay inward, as would be expected due to interactions with the magnetic field as hydrogen becomes increasingly more conductive with increasing pressure (Cao & Stevenson, 2017). For instance, winds at high latitude could not be included in this method, because they would correspond to cylinders extending all the way through the center of the planet. Differential rotation on cylinders is also north-south-symmetric by definition, so the odd $J_{n}$ are identically zero and cannot be modeled. The models presented in this paper are subject to these limitations.

The second method starts with a gravity solution assuming uniform rotation, using CMS or a similar method, and then uses the thermal wind equation (Kaspi, 2013; Galanti et al., 2017) or the gravitational thermal wind equation (Kong et al., 2013) to calculate a correction to the density and gravitational moments. While this introduces additional approximations and does not produce a self-consistent solution for the gravitational field, it allows for more flexible wind fields, including cylinders of finite depth and flows with north-south asymmetries. Iess et al. (2019) includes calculations in which the observed $J_{n}$ are calculated with a decaying wind profile based on the observed cloud-level winds.

Nevertheless, the models with differential rotation on cylinders that do not decay with depth are an important class of endmember models to consider for two reasons. They fit all even gravity moments measured by the Cassini spacecraft and they are fully self-consistent, which means that predictions for the core mass, composition of the envelope and rotation profile will be obtained from just one theory.

I.2 Interior Model Background

Interior models of Saturn, like the ones presented here, have previously been fitted to gravity data from Voyager (Gudkova & Zharkov, 1999; Guillot, 1999; Saumon & Guillot, 2004) and pre-Grand Finale Cassini data (Helled & Guillot, 2013; Nettelmann et al., 2013). In all cases they take into account a reduction of helium mass fraction ($Y$) in the outer envelope arising from the immiscibility and rainout of helium (Stevenson & Salpeter, 1977), although there are some differences in the degree of rainout considered. The models differ primarily in the material equations of state used, whether the heavy element concentrations ($Z$) are homogeneous or inhomogeneous between the inner and outer envelope, and whether they consider differential rotation. The range of predicted core masses decreased from $\sim$10 – 25 to $\sim$5 – 20 Earth Masses when models were fitted to Galileo-era and pre-Grand Finale Cassini gravity data (Fortney et al., 2016), and some models considering inhomogeneous $Z$ had no central core at all (Helled & Guillot, 2013).

One persistent issue for modelling Saturn’s interior has been the uncertainty of the planet’s deep rotation rate, due to the near-perfect alignment of the magnetic field dipole with the rotation axis. Given this uncertainty, we constructed ensembles of models for four published rotation periods: 10:32:45 h (Helled et al., 2015), 10:39:22 h (Desch & Kaiser, 1981), 10:45:45 h (Gurnett et al., 2005), and 10:47:06 h (Giampieri et al., 2006). We also considered a very short rotation period of 10:30:00 h in order to make the following calculation more robust. An independent constraint on the rotation are measurements of the planet’s degree of flattening (oblateness) (Lindal et al., 1985). In Section III.2, we use this information to derive a new estimate for Saturn’s deep rotation period that is fully consistent with our interior models, CMS method, and the Voyager oblateness measurements.

II Methods

II.1 Interior models

Since planets cool by convection, models are typically constructed under the assumption that most regions in their interiors are adiabatic. However, novel ideas based on double-diffusive convection have also been considered (Leconte & Chabrier, 2013; Nettelmann et al., 2015). One example of non-adiabatic behavior occurs at high pressure, where hydrogen and helium are predicted to become immiscible because hydrogen turns metallic while helium remains an insulating fluid (Stevenson & Salpeter, 1977), leading to a region of helium rain. Following earlier work (Wahl et al., 2017c; Iess et al., 2019), we assume four-layer models with an outer molecular and an inner metallic envelope, separated by a helium rain layer, along with a dense core at the center of the planet, as illustrated in Fig. 2. In both envelope layers, an adiabat consistent with ab initio simulations of hydrogen-helium mixtures (Vorberger et al., 2007; Militzer, 2013; Militzer & Hubbard, 2013) is determined. Each adiabat is characterized by an entropy, $S$, a helium mass fraction, $Y$, and a mass fraction of heavy elements, $Z$. We adopt the phase diagram for hydrogen-helium mixtures as derived by Morales et al. (2009), and assume that helium rain occurs wherever the $P$-$T$ barotrope falls within the region of immiscibility in Fig. 3.

We treat the helium rain layer as a smooth transition from the parameters in the outer envelope ($S_{\rm mol}$, $Y_{\rm mol}$, $Z_{\rm mol}$) to inner envelope ($S_{\rm met}$, $Y_{\rm met}$, $Z_{\rm met}$) across a range of pressures $P_{1}$ to $P_{2}$, defined by the intersections of the adiabat with the immiscibility curve. A summary of our model parameters is given in Tab. 4. A collection of representative barotropes are shown in Fig. 3.

Various core masses and radii are considered, but are not independent, since the total mass of the core and envelope must match that of Saturn. We first assumed fractional radius of 0.2 and later refined the core radii by assuming either a terrestrial iron-silicate composition (0.325:0.675) or a solar iron-silicate-water ice composition (0.1625:0.3375:0.5). We find the fractional core radii of $r_{C}$=0.188 and 0.231 respectively to be consistent with these two compositions. We derived these core radii by adopting the additive volume rule for homogeneous mixtures in combination with the equations of state for iron, MgSiO3 and water ice reported in Seager et al. (2007) and Wilson & Militzer (2014) that relied on experimental data and results from ab initio simulations.

For each set of model parameters, the CMS method finds a shape and gravitational field for the planet consistent with a prescribed rotation rate.

The distribution of helium across the rain layer is represented by a gradual gradient with depth between $Y_{\rm mol}$ and $Y_{\rm met}$. Thus a value of $Y_{\rm mol}$, up to the solar helium fraction $Y$=0.274 (Lodders, 2003), is considered and a consistent $Y_{\rm met}$ above the solar fraction is determined such that total, planet-wide helium mass fraction is conserved. The entropy of the outer envelope adiabat $S_{\rm mol}$ is chosen to be consistent with the observed temperature 142.7 K at 1 bar  (Lindal et al., 1981).

II.2 Concentric Maclaurin Spheroid Method

The literature on the problem of the shape and gravitational potential of a liquid planet in hydrostatic equilibrium (also referred to as the theory of figures, TOF) extends back centuries Jeans (2009). Most geophysical implementations of TOF use a perturbation approach, by finding the response, to various orders, to a small perturbation of the potential from spherical symmetry. For a discussion of perturbation TOF, see Zharkov & Trubitsyn (1978).

Hubbard (2012) developed a non-perturbative numerical method, based on potential theory (Tassoul, 2015), for calculating the self-consistent shape and gravitational field of a constant density, rotating fluid body to high precision. This method was generalized to approximate a barotropic pressure-density relationship, discretizing the interior into a series of concentric constant-density (Maclaurin) spheroids (CMS) by Hubbard (2013). The spheroids comprise constant-potential level surfaces, deformed in two dimensions for permanent rotation about a fixed axis, and in three dimensions if a tidal potential is included (Wahl et al., 2017b). Thus, the surface of every spheroid is a surface of constant potential, density, pressure, temperature, and composition. The CMS method is non-perturbative and thus more general than methods that approximate the level surfaces as perturbed ellipsoids. The CMS method has been benchmarked against an independent, non-perturbative numerical method (Wisdom & Hubbard, 2016).

In this paper, we introduce an accelerated version of the CMS method, in which the shape of a subset of spheroids is calculated explicitly, with the shape of most spheroids obtained through interpolation of the radius. As we will show, this leads to a much more efficient algorithm for the same level of precision of the predicted gravity field. The acceleration technique enables us to construct ensembles of Saturn’s interior models with Monte Carlo sampling and to perform proof-of-principles CMS calculations with a large number of layers ($N_{L}$) $\sim 10^{5}$. Both would not have been feasible without acceleration of the method.

As noted by Debras & Chabrier (2018), while a model with a given number of spheroids generates an external gravity potential to a numerical precision of at least $10^{-12}$ (much better than Juno or Cassini measurement precision), the precision to which it approximates the smooth $\rho(P)$ barotrope is limited by the number of layers. This leads to an $N_{L}$-sensitivity of the generated gravity potential that is larger than the uncertainty in the measured potential, as initially quantified by Wisdom & Hubbard (2016). The acceleration to the CMS method helps us rectify any uncertainty from discretization, allowing a much smoother discretization of the barotrope while the more computationally expensive part of the method is kept to a manageable number of layers.

II.3 Self-consistent Shape and Gravity with CMS

The CMS technique, based on potential theory, allows one to describe the interior of planets under the assumption of hydrostatic equilibrium. Baroclinic effects are excluded from consideration, which implies that the temperature of a fluid parcel is only a function of its pressure, $T(P)$. While this is well justified in the deep interior, it is more of an approximation at the 1 bar level when we relate the temperature of fluid parcels near the equator with those in the less irradiated polar regions. Under this assumption, we combine $T(P)$ with a realistic equation of state, $\rho=\rho(P,T)$, of a mixture of hydrogen, helium, and a small amount of heavier elements in order to establish a barotrope, a unique density-pressure relation $\rho(P)=\rho(P,T(P))$. This assumes knowledge of the composition as a function of pressure.

In hydrostatic equilibrium, the pressure, $P$, the mass density, $\rho$, and the total potential, $U$, at any point in the planet’s interior are related by

[TABLE]

The sign of the potentials is chosen such that forces are given by $F=+\nabla U$. In the co-rotating frame of the planet, the total potential, $U$, includes contribution from the self-gravity, $V$, and the centrifugal term, $Q$,

[TABLE]

which we discuss in detail in the two following sections.

For a planet with a uniform rotation rate, it is convenient to describe the relative strength of of the rotational perturbation in terms of the parameter

[TABLE]

where $\omega$ is the rotation rate, $G$ is the universal gravitational constant, and $M$ and $a$ are the mass and equatorial radius of the planet. Since CMS theory is non-perturbative, in principle the results are valid to all powers of $q_{\rm rot}$.

It follows that the pressure, density and potential can be expressed in dimensionless, planetary units (PU):

[TABLE]

We label the $N_{L}$ spheroids with the indices $i=0,1,2,\dots,N_{L}-1$, with $i=0$ corresponding to the outermost spheroid and $i=N_{L}-1$ corresponding to the innermost spheroid. All models presented here are symmetric with respect to the axis of rotation. We neglect any non-axisymetric contributions to potential, such as tidal perturbation by a satellite (Wahl et al., 2017a, 2016). So the shape of every spheroid $i$ can be described by a function $r_{i}(\mu)$ where $r_{i}$ is the distance from the planet’s center and $\mu=\cos(\theta)$ is a function of the polar angle, $\theta$. We assume throughout its interior, the planet is north-south symmetric, which implies, $r_{i}(\mu)=r_{i}(-\mu)$.

It is convenient to introduce a normalized shape function,

[TABLE]

where $r_{i}(0)$ is equatorial radius of $i$th spheroid. $\zeta_{i}(\mu)$ will approach unity for non-rotating planets. Furthermore, we define $\lambda_{i}\equiv r_{i}(0)/r_{0}(0)$ to be the ratio of the equatorial radius of the $i$th to the outermost spheroid. Note that $r_{0}(0)\equiv a$. These choices are illustrated in Fig. 4.

Hydrostatic equilibrium requires that the density increases monotonically with depth and thus with spheroid index $i$. We can define $\delta_{i}$ to be the density difference between two adjacent spheroids,

[TABLE]

This parameterization of density has the added benefit of naturally handling discontinuities in $\rho$, as would be expected for compositionally distinct layers.

We represent the shape functions, $\zeta_{i}(\mu)$, on a grid of $N_{\mu}$ points, $\mu_{m}$, such that $\zeta_{im}\equiv\zeta_{i}(\mu_{m})$. The CMS method refines the shape functions through an iterative procedure until the potential on every spheroid surface is constant and the equation of hydrostatic equilibrium is satisfied (Eq. 1). In the current implementation, we keep equatorial radii of every spheroid fixed, $r_{i}(0)=\lambda_{i}a$, while the remaining spheroid points are adjusted until a self-consistent solution has been found.

We start the iterations with all spheroids to be perfect spheres and thus initialize all normalized shape functions to unity, $\zeta_{im}=1$. A given set of spheroids defines a mass distribution and thus a gravity field. We can define a function $U_{i}(\zeta,\mu)$ to calculate the total potential on the surface of spheroid $i$. The spheroid shape has converged if $U_{i}(\zeta,\mu)$ is the same for all $\mu$. However, at the beginning there will always be significant deviations that we can encapsulate in a function,

[TABLE]

that compares the potential at $\zeta_{im}$ and $\mu_{m}$ with that of reference point on the equator of spheroid $i$. We compute the derivative $f^{\prime}_{im}(\zeta_{im})=df_{im}(\zeta_{im})/d\zeta_{im}$ analytically and employ a single Newton step to derive an improved value for $\zeta_{im}$ from

[TABLE]

Once the points on all spheroids have been updated, we recalculate the zonal gravitational moments, $J_{n}$, in order to obtain an updated gravity field, $U_{i}$. Assuming hydrostatic equilibrium (Eq. 1), we successively update the pressure on every spheroid

[TABLE]

starting from $P_{0}$ that we keep fixed at 0.1 bar. This value is consistent with the observed gravity harmonics that were normalized to an equatorial radius of $a=60330$ km (Iess et al., 2019).

Next we update the density of every spheroid,

[TABLE]

by evaluating the prescribed barotrope function, $\rho(P)$, for the average of the pressure at the upper and the lower boundaries of a particular spheroid.

After every improvement of the spheroid shapes, $\zeta_{im}$, an update step for the gravity harmonics, the potential, pressure, and spheroid densities follows. These two steps are repeated until all of the moments, $J_{n}$, have converged such that the difference between successive iterations falls below a specified tolerance. Occasionally, we find the convergence of the algorithm to be slow if the shapes oscillate back and forth between two states. We detect such events and bypass them by inserting a regula falsi step.

It is also necessary to have at least one free parameter for a subset of the layers in order to obtain the correct total mass of the CMS model. In our implementation we modify the mass of the central core to achieve this balance.

II.4 Gravitational Potential

The gravitational potential at a vector coordinate, $\mathbf{r}$, due to an arbitrary mass distribution is given by

[TABLE]

In the case of an axisymmetric mass distribution with the center of mass at the origin, the potential can be expanded in the following form (Zharkov & Trubitsyn, 1978),

[TABLE]

where $d\tau^{\prime}=r^{\prime 2}dr^{\prime}\,d\mu^{\prime}\,d\phi^{\prime}$. $P_{n}$ are the Legendre polynomials of order $n$. The gravity harmonics are given by

[TABLE]

$J_{0}$ represents the integral over all mass and has been normalized to equal $-1$ by convention.

Following Hubbard (2013), the self-gravity contribution to the potential is found by expanding Eq. (12) in terms of the interior zonal harmonics, $J_{i,n}$, and the external zonal harmonics, $J_{i,n}^{\prime}$ and $J_{i,n}^{\prime\prime}$, for every spheroid $i$ and order $n$. At the surface of the planet, the observable zonal harmonic is the sum of the moment from every spheroid.

For convenience, the harmonics are normalized by the equatorial radius of the corresponding spheroid

[TABLE]

Following the derivation in Hubbard (2013), we find the normalized interior harmonics

[TABLE]

and the exterior harmonics

[TABLE]

with a special case for $n=2$

[TABLE]

and

[TABLE]

where $M$ is the total mass of the planet given by

[TABLE]

With this description of the planet’s self-gravity in terms of $J_{i,n}$, $J_{i,n}^{\prime}$ and $J_{i,n}^{\prime\prime}$, the expansion of of Eq. 12 for a point on surface $i$ yields

[TABLE]

The gravitation potential on the equator of the outermost spheroid is given by

[TABLE]

where

[TABLE]

are the standard zonal gravity harmonics of the observable surface field in Eq. 14. In practical application of the CMS method, one finds that results converge rapidly with increasing polynomial order, $n$. So we typically terminate the sum over $n$ at 16 or 32.

II.5 Centrifugal Potential

We assume potential theory throughout this work and we are thus restricted to studying two cases: uniform rotation ($\omega$ = constant) and differential rotation on cylinders where the angular frequency, $\omega(l)$, is solely a function of the distance from the rotation axis, $l$. An illustration is shown in Fig. 5. Everywhere the centrifugal force, $\vec{F}=l\omega^{2}\vec{e}_{l}$, is perpendicular to the axis of rotation, which we assume to be the $z$ axis. In potential theory, this force is represented by the centrifugal potential,

[TABLE]

If $\omega$ is constant, one recovers the usual term $Q(l)=\frac{1}{2}l^{2}\omega^{2}$. It is not possible to give the cylinders a finite depth, $H$, within potential theory. Calculations with finite $H$ can be performed with the thermal wind equation (Kaspi & Galanti, 2016; Kaspi et al., 2017, 2018) or the gravitational thermal wind equation (Kong et al., 2013). If one wanted to give the cylinders a finite depth or introduce any other $z$ dependence, $\omega(l,z)$, one would inevitably introduce spurious force terms parallel to the $z$ direction because the derivative $\partial Q/\partial z$ is no longer zero. This force would not be consistent with the assumption that the centrifugal force should be perpendicular to the axis of rotation (Tassoul, 2015). Therefore, the cylinders in our calculations penetrate through the equatorial plane of the planet. As we will later see, this allows us to reproduce the observed winds in the equatorial regions but not those at higher latitudes, because they would involve very deep cylinders with too much mass.

Most simply, one can represent the angular frequency by an expansion in even powers of $l$,

[TABLE]

where $\omega_{0}$ is the rotation rate in the deep interior and the expansion coefficients, $c_{2i}$, present the differential part. These coefficients need to be optimized jointly with the parameters of our interior model in order to reproduce the gravity coefficients that were measured by the Cassini spacecraft. While the expansion in Eq. 25 may be convenient for analytical work, we found this functional form to be impractical for numerical optimizations. If one changes one coefficient in the expansion, rotation of all fluid parcels is affected. Changing the rotation rate in a small interval of $l$, requires changing several coefficients in a coordinated fashion. Such inter-dependencies are detrimental for the efficiency of any optimization algorithm. We therefore represent the angular frequency, $\omega(l)$, from $l=0\ldots 1$ by a spline function with a fixed number of knots, $l_{k}$, on which we adjust the frequency $\omega(l_{k})$. In this formulation, a change of $\omega(l_{k})$ will only affect fluid parcels between $l_{k-1}$ and $l_{k+1}$, which greatly simplifies the optimization.

We obtained good results with 11 and 21 knots. Furthermore, in our Monte Carlo (MC) calculations and simplex optimizations, we observed that the angular frequency $\omega(l_{k})$ for radii interior to $l_{k}<0.7$ never deviated from $\omega_{0}$, presumably because the associated cylinders were so deep and involved too much mass. Based on these observations, we exclude the $\omega(l_{k})$ values for small $l$ from the optimization and set $\omega(l_{k}<0.7)=\omega_{0}$ instead.

II.6 Acceleration of the CMS method

Among numerical methods to solve partial differential equations (PDE), one distinguishes between finite difference and finite element techniques (Morton & Mayers, 2005). In the former approach, one approximates the derivatives in the PDE by computing differences between two adjacent points on the integration domain. In the more sophisticated finite element approach, one also considers the properties of the interior of every integration interval. This typically enables one to derive a more accurate solution than is possible with finite difference approaches, when the two methods are compared for the same grid resolution.

The acceleration of the CMS method, that we will now introduce, is comparable to switching from the finite difference to a finite element approach. The goal is to reduce the primary discretization error of the CMS methods that arises from the approximation that the density changes in a step-wise fashion from one spheroid to the next. The acceleration becomes possible because each CMS iteration has two parts that have very different computational costs. The expensive part (Eq. 8) involves updating the shape of every spheroid represented by the variables $\zeta_{jm}$ for a given gravity field. In the second, comparatively cheap step, one updates the interior and exterior gravity harmonics in Eqs. 16-19 for the current spheroid geometry. As it turns out, the accuracy of the computed gravity harmonics depends sensitively on the number of spheroids, $N_{L}$, which determines how precisely the smooth density profile in the planet’s interior is approximated by the step-wise representation of the nested constant-density spheroids.

The core idea behind the acceleration is to only compute the spheroid shape explicitly at every $n_{\rm int}$ layers. For the $(n_{\rm int}-1)$ layers in between, we interpolate the shape functions $\zeta_{im}$ as a function of $\lambda_{i}$ at constant $\mu_{m}$. This $\zeta_{im}$ update is the most expensive part of the CMS calculation and scales like $N_{L}\times N_{\mu}$ while the other parts of the calculation all scale like $N_{L}$. Therefore, we evaluate the other parts of the calculation over the entire set of $N_{L}$ spheroids as before. The cost of the spline interpolation is negligible compared to the explicit updates of the $\zeta_{im}$ points according to Eq. 8. The inner and outermost spheroids are always updated explicitly to avoid extrapolations.

Instead of updating $\zeta_{im}$ for $N_{L}$ layers, we only need to update $N_{L}/n_{\rm int}$ layers 111To keep the following analysis simple, we write $N_{L}/n_{\rm int}$ for the number of layers that we treat explicitly while it is in fact $(N_{L}-1)/n_{\rm int}+1$.. The reduction in computational cost can be reinvested into increasing the total number of layers. As we will show, the accuracy of an accelerated CMS computation with a total layer number of $N_{L}^{\rm acc}=N_{L}^{\rm org}\times n_{\rm int}$ will be much higher than that of the original calculation of $N_{L}^{\rm org}$ layers, while both have comparable computational cost. The computation of all gravity harmonics is be performed with all $N_{L}^{\rm acc}$ layers, which significantly improves the accuracy compared with the original calculations with $N_{L}^{\rm org}$ layers.

In order to analyze the accuracy and the performance of our acceleration technique, we constructed a representative model for the interior structure of Saturn. For this analysis, we assume uniform rotation and performed calculations for a variety of layer numbers with interpolation parameter, $n_{\rm int}$, ranging from 2 to 128. The results of the original method without acceleration are recovered for $n_{\rm int}=1$. The resulting gravity coefficients that were computed with and without acceleration are compared in Tabs. 5 and 6 for different layer numbers. One finds that all gravity harmonics converge smoothly as a function of layer number, which allows one to extrapolate to $N_{L}\to\infty$. We infer $J_{n}(\infty)$ by employing the following semi-linear fit function:

[TABLE]

For every gravity coefficient, $J_{n}$, we adjust the fit parameter $J_{n}(\infty)$ and derive the linear fit coefficients $A$ and $B$ until we have obtained the best possible match to the $J_{n}(N_{L})$ data set. The extrapolated values, $J_{n}(\infty)$, are included in Tab. 6.

Having access to extrapolated values, $J_{n}(\infty)$, allows us to study how the discretization error decays with increasing $N_{L}$ and to evaluate the effectiveness of the acceleration scheme. All curves in Fig. 6 show that the discretization error decays quadratically as $N_{L}^{-2}$. The top panel shows the behavior of the original method before any acceleration was introduced. For $J_{12}$, one finds that only 32 layers are needed for the discretization error to be less than the error bar of the Cassini measurements because the uncertainty is comparatively large for this gravity coefficient. Conversely $J_{2}$ has been measured with a much higher precision and even CMS calculations with 4096 layers are not sufficient to meet the accuracy of the measurements.

The middle panel of Fig. 6 shows the discretization error of accelerated CMS method with acceleration factor, $n_{int}=16$. The results show that calculations with 512 explicit layers ($N_{L}^{\rm tot}=16384$) are sufficiently accurate to reduce the discretization error of computed gravity coefficients below the uncertainty level of the Cassini measurements. This demonstrates that with the acceleration technique is very effective and enables us to match the accuracy of Juno and Cassini measurements within the CMS framework.

In the lower panel of Fig. 6, the discretization error of different $J_{n}$ have been combined in order to compare results for different acceleration factors, $n_{\rm int}$. The figure confirms that an increase in $n_{\rm int}$ leads to a significant reduction of the discretization error when results are compared for the same number layers that are treated explicitly, $N_{L}^{tot}/n_{\rm int}$, which is also a measure of the computational cost.

The lower panel of Fig. 6 also compares the discretization error that arises from two different $\lambda$ grids. The choice of $\lambda$ grid has an impact on how many spheroids are needed to reach a certain level of accuracy. We show results derived with an earlier $\lambda$ grid from Wahl et al. (2017a), which was constructed by employing a denser mesh of spheroids in the atmosphere and outer layers of the planet where the density changes the most. We then developed an alternate approach with the aim of constructing an optimal $\lambda$ grid that further reduces the discretization error. This error arises from contrast in density between two adjacent spheroids. To minimize this error, we construct a $\lambda$ grid such that the relative difference in density is the same for all pairs of adjacent spheroids throughout the planet. This automatically places more layers in the atmosphere, where the density changes most rapidly. We construct our optimized $\lambda$ grid by starting from a converged CMS calculation with our original grid, which provides us with series of $\rho(\lambda_{i})$ points that we can interpolate. We construct a geometric grid of $\rho_{i}$ values that the spans the interval between the lowest and highest density in our model while keeping $\rho_{i+1}/\rho_{i}$ constant. We derive our optimized $\lambda_{i}$ grid by solving $\rho(\lambda_{i})=\rho_{i}$. In Fig. 7 and Tab. 1, we compare the original and optimized $\lambda$ grids. During the optimization, more grid points are placed in outer region of the planet where the density changes most rapidly. However, the inset of Fig. 7 shows that the slopes of the two grid functions is very similar near $\lambda=0$. In this region the grid space should be a fraction of the scale height of the atmosphere.

In limit of $N_{l}\to\infty$, CMS calculations with both $\lambda$ grids will converge to identical results because the discretization errors will gradually dimish in every part of the interior. However, an optimized $\lambda$ grid may approach this limit more rapidly. The lower panel of Fig. 6 shows that our optimized $\lambda$ grid reduces the discretization error by a factor of 2.3 when compared to our original grid for the same number of spheroids. For this reason, we employ the optimized grid in all following calculations.

II.7 Planet models with polytrope index 1

Here we revisit standard planetary interior models that approximate the equation of state throughout the interior by a polytropic equation of state, $P(\rho)=K\rho^{1+1/n}$ with index $n=1$. The constant $K$ is adjusted so that planet’s total mass equals 1. Under these assumptions, potential and density are proportional and the planet’s surface is given by $P=\rho=0$. Wisdom & Hubbard (2016) studied the properties of such planet models in great detail and compared the predictions from the consistent level curve (CLC) technique and from the CMS method. Here we present a comparison with our accelerated CMS approach, which allows us to control density discretization error more carefully. We benchmark our results against Wisdom & Hubbard (2016) using the identical value of $q_{\rm rot}$=0.089195487.

In Fig. 8, we show how discretization errors decay with increasing number of spheroids. Overall the behavior is similar to that of our more realistic Saturn interior model in Fig. 6.

We choose a acceleration factor of $n_{\rm int}=256$ and performed a set of polytrope index 1 model calculations with increasing precision. The number of explicitly treated layers, $N_{L}/n_{\rm int}$ were varied between $2^{2}$ and $2^{9}$, which brought up the total number of layers to 131072 in our largest calculations, which is an increase of three orders of magnitude compared to earlier CMS calculations. We analyze how our results improved with increasing layer number and report the converged digits in Tab. 2. The agreement with the CLC predictions is excellent. All coefficients $J_{2}$ through $J_{20}$ agree to 6, 7, or 8 significant digits, which is a better agreement than was reported in Wisdom & Hubbard (2016) where predictions from the CLC approach and the non-accelerated CMS method were compared.

II.8 Parameter Optimization

The primary goal of the model optimization is the generation of Saturn interior models that reproduce the observed gravity harmonics. The agreement between models and observations is typically expressed in some form of a $\chi^{2}$ function. Here we use,

[TABLE]

where $\delta J_{2i}^{\rm observed}$ are the 1-$\sigma$ uncertainties in the observations. Typically $J_{2}$ is measured with much higher precision than the higher order harmonics. To deal with this imbalance, we find solutions that satisfy $J_{2i}^{\rm observed}=J_{2i}^{\rm model}$ exactly by adjusting one model parameter like $Z_{\rm mol}$ or $Z_{\rm met}$ before $\chi^{2}_{J}$ is evaluated. This optimization is performed for converged CMS models that have reached hydrostatic balance and have matched the planet’s total mass by adjusting the core mass.

While Eq. 27 is certainly the most important optimization criterion, there are a number of other well motivated constraints to consider. For example, one would want to guide to the parameter optimization towards models with pressures $P_{1}$ and $P_{2}$ are close to the assumed immiscibility curve in Fig. 3. From the assumed molecular and metallic adiabats, we can infer the temperatures $T_{1}$ and $T_{2}$ that correspond to both pressures. For both pairs $P_{1}$-$T_{1}$ and $P_{2}$-$T_{2}$, we find the closest points on the immiscibility curve, $P^{*}_{1}$-$T^{*}_{1}$ and $P^{*}_{2}$-$T^{*}_{2}$, that minimize the following immiscibility penalty function,

[TABLE]

before we add the resulting minimum value to the total $\chi^{2}$. $C_{P}$ and $C_{T}$ are weights that must be balanced with those in other $\chi^{2}$ terms. We set $C_{T}/C_{P}=2$. We chose not to square the individual terms in Eq. 28 because, without an experimental confirmation of our immiscibility curve, we do not want large deviations to enter quadratically. In Fig. 3, we shows some representative models to illustrate how much variation is in the $P_{1}$-$T_{1}$ and $P_{2}$-$T_{2}$ in our ensemble of models. Implicitly the $\chi^{2}_{\rm H-He}$ term also introduces a penalty for metallic adiabats that are too hot to be compatible with the assumed immiscibility curve.

Upon first introducing differential rotation into our CMS models, we realized that a super-rotating equatorial jet improved the match to observed gravity harmonics considerably. Furthermore, for $l\geq 0.8$, the inferred rotation profile was compatible with the wind speeds that were derived from tracking the clouds in Saturn’s atmosphere (Sanchez-Lavega et al., 2000; García-Melendo et al., 2011). From this point on, we favored models that matched those observations by introducing the following cloud penalty function,

[TABLE]

where we sum over the knots, $l_{k}$, in the outer region of our rotation profile that often lead to good agreement with the cloud tracking observations. $C_{\rm clouds}$ plays the role of a weight.

Finally we introduce one more penalty function,

[TABLE]

that favors smooth rotation profiles by penalizing large values in the second derivative of our rotational profile.

We assume $Y_{\rm met}\geq Y_{\rm mol}$ and $Z_{\rm met}\geq Z_{\rm mol}$, because we assume that helium rain can only lead to an enrichment of the metallic layer in helium and in heavy elements (Wilson & Militzer, 2010). We also constrain the helium abundance of the entire envelope to match solar proportions.

We add Eqs. 27 through 30 to obtain one total $\chi^{2}$ function that we employ to optimize the model parameters in Tab. 4. This turns out to be a very challenging optimization problem, because many parameters are strongly coupled and some optimization criteria are interdependent. We use the simplex algorithm (Press et al., 2001) for the optimization since it does not require any derivatives of $\chi^{2}$ with respect to the optimization parameters, which are not available in analytical form. With this algorithm, it was very challenging to generate models that matched observed gravity data. In many cases, the algorithm gets stuck in a local minimum222We computed the derivative numerically and tested the BFGS (Jacobs, 1977) optimization algorithm but this did not lead to an algorithm that is more efficient overall because of the cost of computing the derivative with finite differences.. Fig. 9 shows a couple of examples of the $\chi_{J}^{2}$ evolution during the simplex optimization. However, in 17 independent cases, the optimization succeeded and we were able to match the gravity harmonics within the uncertainties of the observations. We subsequently used these 17 solutions as starting points for Markov chain Monte Carlo calculations in order to map out the allowed parameter regions. We confirmed that all 17 original solutions belong to the same parameter region and one can go smoothly from one to the other. This provides strong evidence that the entire solution space is connected.

II.9 Effects of an upper atmosphere

All CMS calculations presented so far, start from an outermost spheroid with the pressure of 0.1 bar that was anchored at the equatorial radius $a$. We had thereby neglected the effects of the tenuous upper atmosphere that extends from the 0.1 bar level out into space. To study the effects of this upper atmosphere quantitatively, we added 64, 128, and 256 outer spheroids to CMS calculations with 512, 1024, and 2048 layers, respectively. The number of the additional spheroids was chosen such that range of pressure extended down to at least 1 mbar. The original outer spheroid is still associated with a pressure of 0.1 bar and remains anchored at the equatorial radius $a$. For all spheroids interior to this spheroid, we update the pressure according to Eq. 9 as we did before. However, for all the additional exterior spheroids, we update the pressure with decreasing spheroid index according to,

[TABLE]

For simplicity, we assume an isothermal upper atmosphere with a temperature set to the value at 0.1 bar. In all other respects the additional exterior spheroids are treated in the same way as the interior spheroids. In principle, the temperature treatment could be made more realistic but, as we will show, the effects on the computed gravitational moments are negligible because there is so little mass outside of the 0.1 bar level.

Our extended CMS calculations converged to the same level of accuracy as they had previously. The pressure of the new outermost spheroid converged to 0.9 mbar. The fractional mass outside of 0.1 bar level was found to be only $7.5\times 10^{-8}$. This mass correction can also be interpreted as a change to the gravity coefficient $J_{0}$ (see Eq. 14), which helps one to gauge the magnitude of the correction to the other gravity coefficients. In table 3, we provide the differences in the gravitational moments between our CMS calculation that included an extended atmosphere to 0.9 mbar and our original calculations that terminate at 0.1 bar. All values decay smoothly with increasing order $n$. For all the gravity coefficients that have been determined by the Cassini spacecraft, $J_{2}$ through $J_{10}$, one finds that the correction due to the upper atmosphere is at least 18 times smaller than the uncertainties of the Cassini measurements (Iess et al., 2019). For this reason, we conclude that our standard CMS calculations starting from the 0.1 bar level are sufficiently accurate for this study.

III Results

III.1 Predictions for interior parameters

In Fig. 5, we plot the rotational profiles that have emerged from our Monte Carlo calculations. Two prominent features are common to all models. There is a super-rotating equatorial jet in the equatorial region that rotates up to 4% faster than the deep interior. This behavior is in agreement with the observed wind speeds from tracking the cloud motion in Saturn’s visible atmosphere (Sanchez-Lavega et al., 2000; García-Melendo et al., 2011) and we have thus favored the sampling of such models by introducing the term $\chi^{2}_{\rm clouds}$ in Eq. 29. At a distance of approximately 50 000 km from the axis of rotation, our models require a sub-rotating region with a flow about 1% slower than in the deep interior. This feature is not observed in the cloud motion at the surface, but is a common feature to all of our models that match the Cassini gravity harmonics. Both the super-rotating equatoral jet, and the sub-rotating feature are present regardless of the value we assume for the rotation period of the deep anterior.

In Fig. 11, we compare the predictions from ensembles of models that we generated with MC sampling for a range of core radii and rotation periods for the deep interior. In panel (a), we plot the amount of heavy elements in the envelope against the core mass. When one compares models for the same core radius of $r_{C}=0.2$, a simple trend emerges. With increasing rotation period, the amount of heavy elements in the atmosphere decreases from approximately 4-fold to 1.2-fold the solar value ($Z_{\rm solar}=0.0153$ Lodders (2010)) while the core mass increases from approximately 15.3 to 16.9 Earth masses. Larger variations in the predicted core masses are seen when the fractional core radius is varied between 0.188 (rocky composition) and 0.231 (rock-ice core in Callisto’s proportion (Kuskov & Kronrod, 2005). A smaller core radius leads to a smaller core masses because the H-He mixture that surrounds that core is exposed to higher pressure, which increases its density and lets it mimic the behavior of the dense core to a larger extent than this is the case in models with larger core radii. Therefore, the uncertainty of the core composition is the primary reason why the predicted core masses vary between 15 and 18 Earth masses.

In Fig. 11b, we plot the combined enrichment in helium and heavy elements in the metallic layer against the entropy in this layer. Since a higher entropy implies a higher temperature and thus a slightly lower density, the enrichment rises with increasing entropy. We find the models with a very long rotation period of 10:45:45 h and 10:47:06 h are confined to a very narrow region of available parameter space predicting the lowest enrichment and the highest entropy for the metallic layer. The long period models appear similarly confined in Fig. 11c where they predict almost no helium rain had occurred while models with shorter rotation periods predict various amount of helium rain. $Y_{\rm mol}$ values as low as 0.19 are included. $Y_{\rm mol}$ and $Y_{\rm met}$ are tightly correlated in this figure because we assume the envelope overall has a solar helium abundance.

In Fig. 11d, we compare the heavy element abundances in the molecular and metallic layers. Within the model constraint of $Z_{\rm met}\geq Z_{\rm mol}$, a wide range of super-solar enrichments are predicted by our ensembles of models. There are plenty of models with $Z_{\rm met}\approx Z_{\rm mol}$, which is in contrast to recent Jupiter models that required a different amounts of enrichments in the two layers (Wahl et al., 2017c). In Fig. 11e and f, we plot the heavy element against the helium abundances in the molecular and metallic layers, respectively. While both quantities are strongly correlated in the molecular layer, there appears to be much more flexibility in the metallic layer. One reason for this is that a range of $S_{\rm met}$ values are permitted in our models while the entropy in the molecular layer is tied to the temperature at the 1 bar level. In Fig. 11e, one can identify a consistent trend for models with longer rotation periods to predict larger values $Y_{\rm mol}+Z_{\rm mol}$ and thus a slightly higher density for the molecular layer.

The models with shorter rotation periods produce $Y_{\rm mol}$ that are compatible with reanalyzed Voyager measurements of atmospheric helium, $\sim$0.6–0.8$\times$solar (Conrath & Gautier, 2000), while the models with longer rotation rates require less depletion of helium in the outer envelope. Observational constraints on $Z_{\rm mol}$ are uncertain; Fletcher et al. (2009) observed atmospheric methane concentrations consistent with $\sim$9$\times$solar enrichment of carbon. There are no direct measurements of the abundance of atmospheric oxygen, the heavy element with the most significant contribution to the density and by consequence the gravity. Other heavy element ratios observations include both much lower (N/H $\sim$3$\times$ solar) and higher S/H $\sim$13$\times$ solar), although these differences might reflect model dependence in determining the bulk elemental abundance from, or from measuring regions of the atmosphere that our not well mixed (Atreya et al., 2019).

The models presented here predict values of both $Z_{\rm mol}$ and $Z_{met}$ between 1–4$\times$solar for a uniform enrichment of all heavy elements, which is lower than the observed enrichment in carbon. It is worth noting that in-situ measurements of Jupiter’s atmosphere up to 22 bars by the Galileo entry probe showed significant depletion in oxygen compared to carbon, but it is an outstanding question whether this accurately reflects the overall composition of Jupiter’s molecular envelope. The heavy element content predicted by the models are also sensitive to temperature of the adiabat. For Saturn, atmospheric temperature has never been measured in situ. So if $S_{\rm mol}$ is higher than we expect, this tradeoff could account for higher concentrations of heavy elements, without significantly affecting the other model predictions.

III.2 Oblateness and Rotation Period

While the rotation period of Jupiter’s interior has been determined with high accuracy from magnetic field observations, the rotation period of Saturn’s deep interior remains uncertain due to the remarkable alignment between the dipole field and the axis of rotation. However, the rotation period used in CMS calculations of a planet significantly affects its shape. Saturn’s oblateness, $(R_{eq}-R_{polar})/R_{eq}$, has been measured with radio occultation measurements of the Pioneer and Voyager spacecraft (Lindal et al., 1985). Anderson and Schubert (Anderson & Schubert, 2007) constructed interior models with uniform rotation that matched observed oblateness and pre-Cassini gravity coefficients $J_{2}$, $J_{4}$, and $J_{6}$. They derived a rotation period of 10:32:44 h, which is significantly shorter than the system III period of 10:39:22 h (Desch & Kaiser, 1981), as well as Cassini predictions of 10:45:45 h and 10:47:06 h (Giampieri et al., 2006; Helled et al., 2015).

In Fig. 10, we compare models with differential rotation that we constructed for five different rotation periods ranging from 10:30:00 to 10:47:06 h. For all periods, it is possible to construct interior models with differential rotation that match all even gravity coefficients. However, the oblateness sensitively depends on rotation period that is assumed for the planet’s deep interior. In Fig. 10, we compare the oblateness that derived from our models with the radio occultation measurements by the Pioneer and Voyager spacecraft (Lindal et al., 1985). We find that rotation periods of 10:33:34 h $\pm$ 55 s are consistent with these observations, which represents a modest 1-$\sigma$ increase over the earlier determination of 10:32:44 h by Anderson and Schubert (Anderson & Schubert, 2007) who performed a similar analysis based on interior models with uniform rotation. Thus, the 50 s difference can primarily be attributed to effects of differential rotation.

Our determination of Saturn’s rotation rate is in remarkably good agreement with the value of 10:33:38 h${}^{+112\rm s}_{~{}-89\rm s}$ inferred from waves observed in Saturn’s rings (Mankovich et al., 2019), even though the interior models for this analysis were constructed without considering differential rotation.

In Fig. 12, we compare predictions from models with our preferred rotation period of 10:33:34 h to those based on the system III rotation period of 10:39:22 h. We still predict a core mass range from 15 to 18 Earth masses, primarily set by the uncertainty in the core composition. When we compare Figs.  11d and  12b, we find the range of $Z_{\rm met}$ is considerably narrowed if the rotation period is set to 10:33:34 h. Most $Z_{\rm met}$ values now fall between values between 1.8 and 2.5 $Z_{\rm solar}$ while in Fig. 11d, the smaller and larger $Z_{\rm met}$ values came from models with longer and shorter rotation periods, now disfavored because of the oblateness constraint. The $Z_{\rm mol}$ values vary between 1 and 3-fold $Z_{\rm solar}$ as before.

Finally, in Fig. 12c, $Z_{\rm mol}$ and $Y_{\rm mol}$ are now fairly tightly correlated when we assume a rotation period of 10:33:34 h. These predictions can, in principle, be verified with remote observations or by an entry probe on a future missions.

IV Conclusions

We have presented an accelerated version of CMS that allows us to construct planetary interior models with many more layers than before and also enables construction of ensembles of models using Monte Carlo methods to efficiently optimize the parameters of individual models. We have applied this accelerated CMS method to construct models for Saturn’s interiors with differential rotation on cylinders, which permitted us to match the unexpectedly large values of the gravity harmonics $J_{6}$, $J_{8}$, and $J_{10}$ that the Cassini spacecraft measured during its Grand Finale orbits around Saturn. From our interior models we infer that Saturn has a massive core of $\sim$15–18 Earth masses and there are additional heavy elements worth 1.5–5 Earth masses distributed throughout its envelope. In our models, we have also varied the rotation period of Saturn’s deep interior and studied the effects on Saturn’s oblateness. By matching occultation measurements of spacecraft we predict a rotation period of $10:33:34$h $\pm$ 55s for Saturn’s deep interior.