Models of Saturn's Interior Constructed with Accelerated Concentric Maclaurin Spheroid Method
MModels of Saturn’s Interior Constructed with Accelerated Concentric MaclaurinSpheroid Method
B. Militzer
Department of Earth and Planetary Science, University of California, Berkeley, CA 94720 andDepartment of Astronomy, University of California, Berkeley, CA 94720
S. Wahl
Department of Earth and Planetary Science, University of California, Berkeley, CA 94720
W. B. Hubbard
Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721
The
Cassini spacecraft’s Grand Finale orbits provided a unique opportunity to probe Saturn’sgravity field and interior structure. Doppler measurements [21] yielded unexpectedly large valuesfor the gravity harmonics J , J , and J that cannot be matched with planetary interior mod-els that assume uniform rotation. Instead we present a suite of models that assume the planet’sinterior 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 usingthe Concentric Maclaurin Spheroid (CMS) method. We present an acceleration technique for thismethod, which drastically reduces the computational cost, allows us to efficiently optimize modelparameters, map out allowed parameter regions with Monte Carlo sampling, and increases the pre-cision of the calculated J n gravity harmonics to match the error bars of the observations, whichwould be difficult without acceleration. Based on our models, Saturn is predicted to have a densecentral core of ∼ ± I. INTRODUCTION
Although Saturn’s deep interior was not a primary tar-get of the
Cassini spacecraft’s 13-year mission monitor-ing the Saturnian system, the final phase of the missionprovided unprecedentedly precise measurements of theplanet’s gravitational field [21]. This phase, from April 23to Sept. 15, 2017, culminated in 22 Grand Finale orbits,during which the
Cassini spacecraft dived between theplanet and its innermost ring. These measurements werecontemporaneous with the ongoing
Juno mission, whichis providing analogous measurements for Jupiter [8]. Asa result of both studies, the measured gravity fields arefar more precise than ever before, warranting a closerlook at the theory and numerical techniques linking theobserved gravity to the interior density structure of theplanet. Here we present models of Saturn’s interior struc-ture and interior rotation rate, matched to the
Cassini measurements, along with an acceleration technique forthe Concentric Maclaurin Spheroid (CMS) method [20]for calculating a self-consistent shape and gravity field.Prior to
Cassini’s
Grand Finale, the best determina-tion of Saturn’s gravity field was from earlier flyby mis-sions and from perturbations of the orbits of Saturn’s nat-ural satellites in combination with the orbit of
Cassini it-self [23]. However, this yielded significant measurementsof only the first three even zonal harmonics of the field, J , J and J . By contrast, X-band Doppler measure-ments during five of the 22 Grand Finale orbits produce a fit with significant determination of even zonal harmon-ics up to J , as well as odd zonal harmonics J and J [21].The distribution of mass within a planet depends onthe equation of state of hydrogen-helium mixtures at highpressures [38], as well as the radial distribution of heavierelements [48]. The interior density distribution influencesthe observed structure of the gravity field through devia-tions from spherical symmetry arising from rotation andtides. Thus, the measured field can place constraints, al-beit non-uniquely, on the internal structure of the planet.For the rapidly rotating Jovian planets, such terms areprimarily determined by the balance between centrifu-gal and gravitational forces. In the absence of internaldynamics, the density distribution and resulting gravityfield are axisymmetric and north-south symmetric, im-plying that only even zonal harmonics J n contribute tothe gravitational potential.If a planet in hydrostatic equilibrium rotates uniformlylike a solid body, the magnitudes of even zonal harmon-ics decay as | J n | ∼ q n rot , where q rot is the ratio of thecentrifugal and gravity accelerations at the equator. The J n of Jupiter measured by Juno spacecraft are broadlyconsistent with this relationship [8], meaning that it ispossible to find models with a uniform rotation rate thatmatch the observed J n , at least in the absence of otherconstraints, from the hydrogen-helium equation of stateand atmospheric composition. However, Fig. 1 illustrateshow the observed even moments J and higher for Saturn a r X i v : . [ a s t r o - ph . E P ] M a y deviate significantly from the expected relationship. Iesset al. [21] demonstrated that these observations cannotbe reproduced with models that assume uniform rota-tion, and that deep differential rotation [18] is requiredinstead. In this paper we expand upon the interpretationof Iess et al. [21] and introduce new analytical tools forhigh-precision gravity modeling. A. Differential Rotation
Over many years prior to and including the durationof the
Cassini mission, optical tracking of clouds has re-vealed large-scale zonal wind currents with respect to theaverage Saturn atmosphere, in particular a pronouncedeastward jet centered on the equator [11, 45]. However,prior to the gravity measurements discussed here, thedata were insufficient to constrain the depth of such zonalflows, and their effects were not considered in previousmodeling studies of Saturn’s interior [17, 43]. With theGrand Finale gravity data, it becomes possible to testa model in which the cloud-level zonal wind belts aremapped onto cylinders that extend to great depths. If thezonal-wind velocity profile continues to depths of manyscale heights, it will affect the observed gravity field intwo ways. First, it modifies the axisymmetric gravita-tional field, and thus changes the even J n from the val-ues expected for a uniformly rotating body with identicalinternal structure [18]. Second, to the extent that the ve-locity profile is not north-south symmetric, there arises acorresponding asymmetry in the gravity field, manifest-ing itself in non-zero odd J n [25]. The values of J and J reported by Iess et al. [21] thus exhibit the north-southasymmetric component of the differential rotation.There are currently two basic methods for incorporat- Degree, n | J n | × Differentialrotation
Saturn: Cassini measurementsSaturn: uniform rotation modelJupiter: Juno measurementsJupiter: uniform rotation model
FIG. 1. Comparison of the gravity harmonics measured forJupiter and Saturn with predictions from models assuminguniform rotation throughout the entire interior of both plan-ets. The deviations are small for Jupiter while substantial dis-crepancies emerge for Saturn. This illustrates that the effectsof differential rotation are much more important for Saturn. ing differential rotation into gravity models. The firstis to approximate the wind profile as rotation on cylin-ders, which can be described using potential theory [18]and can therefore be integrated directly into the poten-tial used in the CMS simulation [58]. This method hasthe benefit of being fully self-consistent; the dynamiccontribution to the potential modifies the shape of theequipotential surfaces, which feeds back into the calcu-lated gravitational field. The downside is that the windprofile must be constant on cylindrical surfaces and thuscannot decay inward, as would be expected due to in-teractions with the magnetic field as hydrogen becomesincreasingly more conductive with increasing pressure [3].For instance, winds at high latitude could not be in-cluded in this method, because they would correspondto cylinders extending all the way through the centerof the planet. Differential rotation on cylinders is alsonorth-south-symmetric by definition, so the odd J n areidentically zero and cannot be modeled. The models pre-sented in this paper are subject to these limitations.The second method starts with a gravity solution as-suming uniform rotation, using CMS or a similar method,and then uses the thermal wind equation [10, 25] orthe gravitational thermal wind equation [29] to calcu-late a correction to the density and gravitational mo-ments. While this introduces additional approximationsand does not produce a self-consistent solution for thegravitational field, it allows for more flexible wind fields,including cylinders of finite depth and flows with north-south asymmetries. Iess et al. [21] includes calculationsin which the observed J n are calculated with a decayingwind profile based on the observed cloud-level winds.Nevertheless, the models with differential rotation oncylinders that do not decay with depth are an impor-tant class of endmember models to consider for two rea-sons. They fit all even gravity moments measured bythe Cassini spacecraft and they are fully self-consistent,which means that predictions for the core mass, composi-tion of the envelope and rotation profile will be obtainedfrom just one theory.
B. Interior Model Background
Interior models of Saturn, like the ones presented here,have previously been fitted to gravity data from
Voyager [13, 14, 46] and pre-Grand Finale
Cassini data [17, 43].In all cases they take into account a reduction of heliummass fraction ( Y ) in the outer envelope arising from theimmiscibility and rainout of helium [49], although thereare some differences in the degree of rainout considered.The models differ primarily in the material equations ofstate used, whether the heavy element concentrations ( Z )are homogeneous or inhomogeneous between the innerand outer envelope, and whether they consider differen-tial rotation. The range of predicted core masses de-creased from ∼
10 – 25 to ∼ Galileo -era and pre-Grand Finale
Cassini gravity data [9], and some models consideringinhomogeneous Z had no central core at all [17].One persistent issue for modelling Saturn’s interior hasbeen the uncertainty of the planet’s deep rotation rate,due to the near-perfect alignment of the magnetic fielddipole with the rotation axis. Given this uncertainty, weconstructed ensembles of models for four published rota-tion periods: 10:32:45 h [16], 10:39:22 h [6], 10:45:45 h[15], and 10:47:06 h [12]. We also considered a very shortrotation period of 10:30:00 h in order to make the follow-ing calculation more robust. An independent constrainton the rotation are measurements of the planet’s degreeof flattening (oblateness) [32]. In Section III B, we usethis information to derive a new estimate for Saturn’sdeep rotation period that is fully consistent with our in-terior models, CMS method, and the Voyager oblatenessmeasurements.
II. METHODSA. Interior models
Molecular hydrogen (helium depleted) Metallic hydrogen (helium rich)Rock-ice core Helium rain lium A p p r o x i m a t e D e p t h o f W i n d s FIG. 2. Illustration of the four-layer models for Saturn’s in-terior that we constructed in this work. We assume an outermolecular and an inner metallic envelope, separated by a he-lium rain layer, with a dense core at the center of the planet.
Since planets cool by convection, models are typicallyconstructed under the assumption that most regions intheir interiors are adiabatic. However, novel ideas basedon double-diffusive convection have also been consid-ered [31, 42]. One example of non-adiabatic behavioroccurs at high pressure, where hydrogen and helium arepredicted to become immiscible because hydrogen turnsmetallic while helium remains an insulating fluid [49],leading to a region of helium rain. Following earlierwork [21, 55], we assume four-layer models with an outermolecular and an inner metallic envelope, separated by a helium rain layer, along with a dense core at the centerof the planet, as illustrated in Fig. 2. In both envelopelayers, an adiabat consistent with ab initio simulationsof hydrogen-helium mixtures [37, 38, 51] is determined.Each adiabat is characterized by an entropy, S , a heliummass fraction, Y , and a mass fraction of heavy elements, Z . We adopt the phase diagram for hydrogen-heliummixtures as derived by Morales et al. [40], and assumethat helium rain occurs wherever the P - T barotrope fallswithin the region of immiscibility in Fig. 3.We treat the helium rain layer as a smooth transitionfrom the parameters in the outer envelope ( S mol , Y mol , Z mol ) to inner envelope ( S met , Y met , Z met ) across a rangeof pressures P to P , defined by the intersections of theadiabat with the immiscibility curve. A summary of ourmodel parameters is given in Tab. IV. A collection ofrepresentative barotropes are shown in Fig. 3.Various core masses and radii are considered, but arenot independent, since the total mass of the core andenvelope must match that of Saturn. We first assumedfractional radius of 0.2 and later refined the core radiiby assuming either a terrestrial iron-silicate composition(0.325:0.675) or a solar iron-silicate-water ice composi-tion (0.1625:0.3375:0.5). We find the fractional core radiiof r C =0.188 and 0.231 respectively to be consistent withthese two compositions. We derived these core radii byadopting the additive volume rule for homogeneous mix-tures in combination with the equations of state for iron,MgSiO and water ice reported in Seager et al. [47] andWilson & Militzer [57] that relied on experimental dataand results from ab initio simulations. zo T e m p e r a t u r e ( K ) P r e ss u r e a t S a t u r n ' s c o r e - m a n t l e b o un d a r y H-He immiscibility region
Morales et al. (2013) Y o un g , f u ll y m i x e dg i a n t p l a n e t Helium rainstarts S a t u r n ' s e x t e r i o r i s e n t r o p e T ( a r ) = . K Range of interiorisentropes
FIG. 3. Temperature-pressure phase diagram of hydrogen-helium mixtures. The red lines show the interior models with P =10:39:22 h and r C =0.2 in relation to the shaded regionwhere the two fluids are predicted to become immiscible [39].The thin blue lines show various adiabats of an H-He mix-ture for a helium mass fraction of Y=0.245 [34]. The circlesmark the beginning and the end of the immiscibility regionsassumed in different models. For each set of model parameters, the CMS methodfinds a shape and gravitational field for the planet con-sistent with a prescribed rotation rate.The distribution of helium across the rain layer is rep-resented by a gradual gradient with depth between Y mol and Y met . Thus a value of Y mol , up to the solar heliumfraction Y =0.274 [34], is considered and a consistent Y met above the solar fraction is determined such that total,planet-wide helium mass fraction is conserved. The en-tropy of the outer envelope adiabat S mol is chosen to beconsistent with the observed temperature 142.7 K at 1bar [33]. B. Concentric Maclaurin Spheroid Method
The literature on the problem of the shape and gravi-tational potential of a liquid planet in hydrostatic equi-librium (also referred to as the theory of figures, TOF)extends back centuries Jeans [24]. Most geophysical im-plementations of TOF use a perturbation approach, byfinding the response, to various orders, to a small pertur-bation of the potential from spherical symmetry. For adiscussion of perturbation TOF, see Zharkov & Trubit-syn [59].Hubbard [19] developed a non-perturbative numeri-cal method, based on potential theory [50], for calcu-lating the self-consistent shape and gravitational fieldof a constant density, rotating fluid body to high pre-cision. This method was generalized to approximatea barotropic pressure-density relationship, discretizingthe interior into a series of concentric constant-density(Maclaurin) spheroids (CMS) by Hubbard [20]. Thespheroids comprise constant-potential level surfaces, de-formed in two dimensions for permanent rotation abouta fixed axis, and in three dimensions if a tidal potentialis included [54]. Thus, the surface of every spheroid isa surface of constant potential, density, pressure, tem-perature, and composition. The CMS method is non-perturbative and thus more general than methods thatapproximate the level surfaces as perturbed ellipsoids.The CMS method has been benchmarked against an in-dependent, non-perturbative numerical method [58].In this paper, we introduce an accelerated version ofthe CMS method, in which the shape of a subset ofspheroids is calculated explicitly, with the shape of mostspheroids obtained through interpolation of the radius.As we will show, this leads to a much more efficient algo-rithm for the same level of precision of the predicted grav-ity field. The acceleration technique enables us to con-struct ensembles of Saturn’s interior models with MonteCarlo sampling and to perform proof-of-principles CMScalculations with a large number of layers ( N L ) ∼ .Both would not have been feasible without accelerationof the method.As noted by Debras & Chabrier [5], while a model witha given number of spheroids generates an external grav-ity potential to a numerical precision of at least 10 − (much better than Juno or Cassini measurement preci-sion), the precision to which it approximates the smooth ρ ( P ) barotrope is limited by the number of layers. This leads to an N L -sensitivity of the generated gravity poten-tial that is larger than the uncertainty in the measuredpotential, as initially quantified by Wisdom & Hubbard[58]. The acceleration to the CMS method helps us rec-tify any uncertainty from discretization, allowing a muchsmoother discretization of the barotrope while the morecomputationally expensive part of the method is kept toa manageable number of layers. C. Self-consistent Shape and Gravity with CMS
The CMS technique, based on potential theory, allowsone to describe the interior of planets under the assump-tion of hydrostatic equilibrium. Baroclinic effects are ex-cluded from consideration, which implies that the tem-perature of a fluid parcel is only a function of its pressure, T ( P ). While this is well justified in the deep interior, itis more of an approximation at the 1 bar level when werelate the temperature of fluid parcels near the equatorwith those in the less irradiated polar regions. Under thisassumption, we combine T ( P ) with a realistic equationof state, ρ = ρ ( P, T ), of a mixture of hydrogen, helium,and a small amount of heavier elements in order to es-tablish a barotrope, a unique density-pressure relation ρ ( P ) = ρ ( P, T ( P )). This assumes knowledge of the com-position as a function of pressure.In hydrostatic equilibrium, the pressure, P , the massdensity, ρ , and the total potential, U , at any point in theplanet’s interior are related by ∇ P = ρ ∇ U. (1)The sign of the potentials is chosen such that forces aregiven by F = + ∇ U . In the co-rotating frame of theplanet, the total potential, U , includes contribution fromthe self-gravity, V , and the centrifugal term, Q , U = V + Q, (2)which we discuss in detail in the two following sections.For a planet with a uniform rotation rate, it is conve-nient to describe the relative strength of of the rotationalperturbation in terms of the parameter q rot = ω a GM , (3)where ω is the rotation rate, G is the universal gravita-tional constant, and M and a are the mass and equa-torial radius of the planet. Since CMS theory is non-perturbative, in principle the results are valid to all pow-ers of q rot .It follows that the pressure, density and potential canbe expressed in dimensionless, planetary units (PU): P pu ≡ a GM P ,ρ pu ≡ a M ρ , and U pu ≡ aGM U. (4)We label the N L spheroids with the indices i =0 , , , . . . , N L −
1, with i = 0 corresponding to the outer-most spheroid and i = N L − i can be described by a function r i ( µ ) where r i is the distance from the planet’s center and µ = cos( θ ) isa function of the polar angle, θ . We assume throughoutits interior, the planet is north-south symmetric, whichimplies, r i ( µ ) = r i ( − µ ).It is convenient to introduce a normalized shape func-tion, ζ i ( µ ) ≡ r i ( µ ) r i (0) ≤ r i (0) is equatorial radius of i th spheroid. ζ i ( µ )will approach unity for non-rotating planets. Further-more, we define λ i ≡ r i (0) /r (0) to be the ratio of theequatorial radius of the i th to the outermost spheroid.Note that r (0) ≡ a . These choices are illustrated inFig. 4. Axis of Rotation Equator
Spheroid i i =0 r i ( μ )= r i ( ) ζ i ( μ ) r i ( μ=0 )=a λ i ζ i ( μ=0 )=1 μ=+1μ=−1 r ( μ=0 )=a ζ ( μ=0 )=1 λ =1 FIG. 4. Illustration of the CMS method and variable defini-tions.
Hydrostatic equilibrium requires that the density in-creases monotonically with depth and thus with spheroidindex i . We can define δ i to be the density difference be-tween two adjacent spheroids, δ i = (cid:40) ρ i − ρ i − , i > ρ , i = 0 . (6)This parameterization of density has the added benefitof naturally handling discontinuities in ρ , as would beexpected for compositionally distinct layers.We represent the shape functions, ζ i ( µ ), on a gridof N µ points, µ m , such that ζ im ≡ ζ i ( µ m ). The CMSmethod refines the shape functions through an iterativeprocedure until the potential on every spheroid surfaceis constant and the equation of hydrostatic equilibriumis satisfied (Eq. 1). In the current implementation, wekeep equatorial radii of every spheroid fixed, r i (0) = λ i a , while the remaining spheroid points are adjusted until aself-consistent solution has been found.We start the iterations with all spheroids to be perfectspheres and thus initialize all normalized shape functionsto unity, ζ im = 1. A given set of spheroids defines a massdistribution and thus a gravity field. We can define afunction U i ( ζ, µ ) to calculate the total potential on thesurface of spheroid i . The spheroid shape has converged if U i ( ζ, µ ) is the same for all µ . However, at the beginningthere will always be significant deviations that we canencapsulate in a function, f im ( ζ im ) ≡ U i ( ζ im , µ m ) − U i (1 , , (7)that compares the potential at ζ im and µ m with that ofreference point on the equator of spheroid i . We com-pute the derivative f (cid:48) im ( ζ im ) = df im ( ζ im ) /dζ im analyti-cally and employ a single Newton step to derive an im-proved value for ζ im from ζ (new) im = ζ im − f im ( ζ im ) f (cid:48) im ( ζ im ) . (8)Once the points on all spheroids have been updated,we recalculate the zonal gravitational moments, J n , inorder to obtain an updated gravity field, U i . Assuminghydrostatic equilibrium (Eq. 1), we successively updatethe pressure on every spheroid P (new) i = P (new) i − + ρ i − ( U i − U i − ) . (9)starting from P that we keep fixed at 0.1 bar. Thisvalue is consistent with the observed gravity harmonicsthat were normalized to an equatorial radius of a = 60330km [21].Next we update the density of every spheroid, ρ (new) i = ρ ( ( P i +1 + P i ) / , (10)by evaluating the prescribed barotrope function, ρ ( P ),for the average of the pressure at the upper and the lowerboundaries of a particular spheroid.After every improvement of the spheroid shapes, ζ im ,an update step for the gravity harmonics, the potential,pressure, and spheroid densities follows. These two stepsare repeated until all of the moments, J n , have convergedsuch that the difference between successive iterations fallsbelow a specified tolerance. Occasionally, we find theconvergence of the algorithm to be slow if the shapesoscillate back and forth between two states. We detectsuch events and bypass them by inserting a regula falsi step.It is also necessary to have at least one free parameterfor a subset of the layers in order to obtain the correcttotal mass of the CMS model. In our implementationwe modify the mass of the central core to achieve thisbalance. D. Gravitational Potential
The gravitational potential at a vector coordinate, r ,due to an arbitrary mass distribution is given by V ( r ) = G (cid:90) d r (cid:48) ρ ( r (cid:48) ) | r − r (cid:48) | . (11)In the case of an axisymmetric mass distribution withthe center of mass at the origin, the potential can beexpanded in the following form [59], V ( r, µ ) = Gr ∞ (cid:88) n =0 P n ( µ ) (cid:90) dτ (cid:48) ρ ( r (cid:48) ) P n ( µ (cid:48) ) (cid:18) r (cid:48) r (cid:19) n (12)= GMr (cid:34) − ∞ (cid:88) n =1 ( a/r ) n J n P n ( µ ) (cid:35) . (13)where dτ (cid:48) = r (cid:48) dr (cid:48) dµ (cid:48) dφ (cid:48) . P n are the Legendre polyno-mials of order n . The gravity harmonics are given by J n = − πM a n +1 (cid:90) − dµ r max ( µ ) (cid:90) dr r n +2 P n ( µ ) ρ ( r, µ ) . (14) J represents the integral over all mass and has beennormalized to equal − J i,n , and the externalzonal harmonics, J (cid:48) i,n and J (cid:48)(cid:48) i,n , for every spheroid i andorder n . At the surface of the planet, the observable zonalharmonic is the sum of the moment from every spheroid.For convenience, the harmonics are normalized by theequatorial radius of the corresponding spheroid (cid:101) J i,n ≡ J i,n λ ni and (cid:101) J (cid:48) i,n ≡ J (cid:48) i,n λ ( n +1) i . (15)Following the derivation in Hubbard [20], we find thenormalized interior harmonics (cid:101) J i,n = − n + 3 2 πM δ i λ i +1 (cid:90) − dµ P n ( µ ) ζ i ( µ ) ( n +3) (16)and the exterior harmonics (cid:101) J (cid:48) i,n = − − n πM δ i λ i +1 (cid:90) − dµ P n ( µ ) ζ i ( µ ) (2 − n ) (17)with a special case for n = 2 (cid:101) J (cid:48) i,n = − πM δ i λ i +1 (cid:90) − dµ P n ( µ ) log( ζ i ) (18)and J (cid:48)(cid:48) i, = 2 πδ i a M , (19) where M is the total mass of the planet given by M = 2 π N − (cid:88) i =0 δ i λ i +1 (cid:90) − dµ ζ i ( µ ) (20)With this description of the planet’s self-gravity interms of J i,n , J (cid:48) i,n and J (cid:48)(cid:48) i,n , the expansion of of Eq. 12for a point on surface i yields V i ( ζ i , µ ) = − ζ i λ i N − (cid:88) j = i ∞ (cid:88) n =0 (cid:101) J j,n (cid:18) λ j λ i ζ i (cid:19) n P n ( µ ) (21)+ i − (cid:88) j =0 ∞ (cid:88) n =0 (cid:101) J (cid:48) j,n (cid:18) λ i ζ i λ j (cid:19) n +1 P n ( µ ) + i − (cid:88) j =0 J (cid:48)(cid:48) j, λ i ζ i . The gravitation potential on the equator of the outermostspheroid is given by V i =0 (1 ,
0) = − ∞ (cid:88) n =0 P n (0) J n , (22)where J n = N − (cid:88) i =0 λ ni (cid:101) J i,n (23)are the standard zonal gravity harmonics of the observ-able surface field in Eq. 14. In practical application ofthe CMS method, one finds that results converge rapidlywith increasing polynomial order, n . So we typically ter-minate the sum over n at 16 or 32. E. Centrifugal Potential
We assume potential theory throughout this work andwe are thus restricted to studying two cases: uniform ro-tation ( ω = constant) and differential rotation on cylin-ders where the angular frequency, ω ( l ), is solely a func-tion of the distance from the rotation axis, l . An illus-tration is shown in Fig. 5. Everywhere the centrifugalforce, (cid:126)F = lω (cid:126)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, Q = (cid:90) l dl (cid:48) l (cid:48) ω ( l (cid:48) ) (24)If ω is constant, one recovers the usual term Q ( l ) = l ω . It is not possible to give the cylinders a fi-nite depth, H , within potential theory. Calculationswith finite H can be performed with the thermal windequation [26–28] or the gravitational thermal wind equa-tion [29]. If one wanted to give the cylinders a finite depthor introduce any other z dependence, ω ( l, z ), one wouldinevitably introduce spurious force terms parallel to the z direction because the derivative ∂Q/∂z is no longer zero.This force would not be consistent with the assumptionthat the centrifugal force should be perpendicular to theaxis of rotation [50]. Therefore, the cylinders in our cal-culations penetrate through the equatorial plane of theplanet. As we will later see, this allows us to reproducethe observed winds in the equatorial regions but not thoseat higher latitudes, because they would involve very deepcylinders with too much mass.Most simply, one can represent the angular frequencyby an expansion in even powers of l , ω ( l ) = ω + c l + c l + c l + c l + . . . , (25)where ω is the rotation rate in the deep interior andthe expansion coefficients, c i , present the differentialpart. These coefficients need to be optimized jointly withthe parameters of our interior model in order to repro-duce the gravity coefficients that were measured by the Cassini spacecraft. While the expansion in Eq. 25 may beconvenient for analytical work, we found this functionalform to be impractical for numerical optimizations. Ifone changes one coefficient in the expansion, rotation ofall fluid parcels is affected. Changing the rotation ratein a small interval of l , requires changing several coeffi-cients in a coordinated fashion. Such inter-dependenciesare detrimental for the efficiency of any optimization al-gorithm. We therefore represent the angular frequency, ω ( l ), from l = 0 . . . l k , on which we adjust the frequency ω ( l k ). In this formulation, a change of ω ( l k ) will onlyaffect fluid parcels between l k − and l k +1 , which greatlysimplifies the optimization.We obtained good results with 11 and 21 knots. Fur-thermore, in our Monte Carlo (MC) calculations andsimplex optimizations, we observed that the angular fre-quency ω ( l k ) for radii interior to l k < . ω , presumably because the associated cylin-ders were so deep and involved too much mass. Based onthese observations, we exclude the ω ( l k ) values for small l from the optimization and set ω ( l k < .
7) = ω instead. F. Acceleration of the CMS method
Among numerical methods to solve partial differentialequations (PDE), one distinguishes between finite differ-ence and finite element techniques [41]. In the formerapproach, one approximates the derivatives in the PDEby computing differences between two adjacent points onthe integration domain. In the more sophisticated finiteelement approach, one also considers the properties ofthe interior of every integration interval. This typicallyenables one to derive a more accurate solution than ispossible with finite difference approaches, when the twomethods are compared for the same grid resolution.The acceleration of the CMS method, that we will nowintroduce, is comparable to switching from the finite dif-ference to a finite element approach. The goal is to re-
FIG. 5. Average of the rotation profiles in our suite of Saturninterior models that match the observed even gravity har-monics. It shows that differential rotation must be severalthousands of kilometers deep. Our models reproduces theEastward equatorial jet that rotates about 4% faster thanthe deep interior. The inset shows an illustration of the cylin-ders. The rotation frequencies inferred by tracking the cloudsin Saturn’s visible atmosphere [11, 45] are shown for compar-ison. duce the primary discretization error of the CMS meth-ods that arises from the approximation that the densitychanges in a step-wise fashion from one spheroid to thenext. The acceleration becomes possible because eachCMS iteration has two parts that have very differentcomputational costs. The expensive part (Eq. 8) involvesupdating the shape of every spheroid represented by thevariables ζ jm for a given gravity field. In the second,comparatively cheap step, one updates the interior andexterior gravity harmonics in Eqs. 16-19 for the currentspheroid geometry. As it turns out, the accuracy of thecomputed gravity harmonics depends sensitively on thenumber of spheroids, N L , which determines how preciselythe smooth density profile in the planet’s interior is ap-proximated by the step-wise representation of the nestedconstant-density spheroids.The core idea behind the acceleration is to only com-pute the spheroid shape explicitly at every n int layers.For the ( n int −
1) layers in between, we interpolate theshape functions ζ im as a function of λ i at constant µ m .This ζ im update is the most expensive part of the CMScalculation and scales like N L × N µ while the other partsof the calculation all scale like N L . Therefore, we evalu-ate the other parts of the calculation over the entire setof N L spheroids as before. The cost of the spline inter-polation is negligible compared to the explicit updates ofthe ζ im points according to Eq. 8. The inner and out-ermost spheroids are always updated explicitly to avoidextrapolations.Instead of updating ζ im for N L layers, we only needto update N L /n int layers [ ? ]. The reduction in compu-tational cost can be reinvested into increasing the totalnumber of layers. As we will show, the accuracy of anaccelerated CMS computation with a total layer numberof N acc L = N org L × n int will be much higher than that ofthe original calculation of N org L layers, while both havecomparable computational cost. The computation of allgravity harmonics is be performed with all N acc L layers,which significantly improves the accuracy compared withthe original calculations with N org L layers.In order to analyze the accuracy and the performanceof our acceleration technique, we constructed a represen-tative model for the interior structure of Saturn. Forthis analysis, we assume uniform rotation and performedcalculations for a variety of layer numbers with interpola-tion parameter, n int , ranging from 2 to 128. The resultsof the original method without acceleration are recoveredfor n int = 1. The resulting gravity coefficients that werecomputed with and without acceleration are comparedin Tabs. V and VI for different layer numbers. One findsthat all gravity harmonics converge smoothly as a func-tion of layer number, which allows one to extrapolate to N L → ∞ . We infer J n ( ∞ ) by employing the followingsemi-linear fit function:log | ∆ J n ( N L ) | ≡ log | J n ( N L ) − J n ( ∞ ) | = A − B log[ N L ](26)For every gravity coefficient, J n , we adjust the fit param-eter J n ( ∞ ) 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 ( ∞ ), areincluded in Tab. VI.Having access to extrapolated values, J n ( ∞ ), allowsus to study how the discretization error decays with in-creasing N L and to evaluate the effectiveness of the ac-celeration scheme. All curves in Fig. 6 show that thediscretization error decays quadratically as N − L . Thetop panel shows the behavior of the original method be-fore any acceleration was introduced. For J , one findsthat only 32 layers are needed for the discretization errorto be less than the error bar of the Cassini measure-ments because the uncertainty is comparatively large forthis gravity coefficient. Conversely J has been measuredwith a much higher precision and even CMS calculationswith 4096 layers are not sufficient to meet the accuracyof the measurements.The middle panel of Fig. 6 shows the discretization er-ror of accelerated CMS method with acceleration factor, n int = 16. The results show that calculations with 512explicit layers ( N tot L = 16384) are sufficiently accurate toreduce the discretization error of computed gravity coef-ficients below the uncertainty level of the Cassini mea-surements. This demonstrates that with the accelerationtechnique is very effective and enables us to match the E rr o r i n e v e n J n = | J n J n | Without acceleration: n int = 1 J J J J J Cassini J Cassini J Cassini J Cassini J Cassini J E rr o r i n e v e n J n = | J n J n | With acceleration: n int = 16 J J J J J Cassini J Cassini J Cassini J Cassini J Cassini J Number of explicitly treated CMS layers, N tot L / n int S u m o f e rr o r s i n e v e n J n = i = | J i J i | Without acceleration: n int =1With acceleration: n int =16With acceleration: n int =16with earlier gridWith acceleration: n int =32With acceleration: n int =64With acceleration: n int =128Sum of Cassini error bars J + J + J + J + J FIG. 6. Discretization error of the gravity harmonics cal-culated with the CMS method as a function of the numberof spheroids. The horizontal lines show the 1- σ uncertaintiesof the Cassini measurements of the even gravity harmonics, J n . The top panels show how the errors decay with increasingnumber of layers for calculations without acceleration. Themid panel displays results obtained with the accelerated CMSmethod where only one in n int = 16 layers are treated explic-itly. 512 explicit layers (total of 8192) are sufficient to reducedthe error in all calculated gravity harmonics below the uncer-tainties of the observations. Without the acceleration, wellover 4000 layers is required for this level of precision, as thetop panel shows. The bottom panel compares results derivedwith different acceleration factors, n int . For n int = 16, the ef-fects of two different λ discretization schemes are compared. accuracy of Juno and
Cassini measurements within theCMS framework.In the lower panel of Fig. 6, the discretization errorof different J n have been combined in order to compareresults for different acceleration factors, n int . The fig-ure confirms that an increase in n int leads to a signif-icant reduction of the discretization error when results TABLE I. Original and optimized λ grids for an interiormodel with 2049 sphoeroids and r C =0.231. The entire tablewill be published online.Spheroid index i Original λ i Optimized λ i . . . . . . . . . are compared for the same number layers that are treatedexplicitly, N totL /n int , which is also a measure of the com-putational cost.The lower panel of Fig. 6 also compares the discretiza-tion error that arises from two different λ grids. Thechoice of λ grid has an impact on how many spheroidsare needed to reach a certain level of accuracy. We showresults derived with an earlier λ grid from Wahl et al.[53], which was constructed by employing a denser meshof spheroids in the atmosphere and outer layers of theplanet where the density changes the most. We then de-veloped an alternate approach with the aim of construct-ing an optimal λ grid that further reduces the discretiza-tion error. This error arises from contrast in density be-tween two adjacent spheroids. To minimize this error,we construct a λ grid such that the relative differencein density is the same for all pairs of adjacent spheroidsthroughout the planet. This automatically places morelayers in the atmosphere, where the density changes mostrapidly. We construct our optimized λ grid by startingfrom a converged CMS calculation with our original grid,which provides us with series of ρ ( λ i ) points that we caninterpolate. We construct a geometric grid of ρ i valuesthat the spans the interval between the lowest and high-est density in our model while keeping ρ i +1 /ρ i constant.We derive our optimized λ i grid by solving ρ ( λ i ) = ρ i . InFig. 7 and Tab. I, we compare the original and optimized λ grids. During the optimization, more grid points areplaced in outer region of the planet where the densitychanges most rapidly. However, the inset of Fig. 7 showsthat the slopes of the two grid functions is very similarnear λ = 0. In this region the grid space should be afraction of the scale height of the atmosphere.In limit of N l → ∞ , CMS calculations with both λ grids will converge to identical results because the dis-cretization errors will gradually dimish in every part ofthe interior. However, an optimized λ grid may approachthis limit more rapidly. The lower panel of Fig. 6 showsthat our optimized λ grid reduces the discretization er-ror by a factor of 2.3 when compared to our original gridfor the same number of spheroids. For this reason, weemploy the optimized grid in all following calculations. i G r i d p o i n t i Original gridOptimized grid with( i +1 )/ ( i ) constant FIG. 7. Comparison of our original and optimized λ grids foran interior model with 2049 spheroids and r C =0.231. G. Planet models with polytrope index 1
Here we revisit standard planetary interior models thatapproximate the equation of state throughout the inte-rior by a polytropic equation of state, P ( ρ ) = Kρ /n with index n = 1. The constant K is adjusted so thatplanet’s total mass equals 1. Under these assumptions,potential and density are proportional and the planet’ssurface is given by P = ρ = 0. Wisdom & Hubbard[58] studied the properties of such planet models in greatdetail and compared the predictions from the consistentlevel curve (CLC) technique and from the CMS method.Here we present a comparison with our accelerated CMSapproach, which allows us to control density discretiza-tion error more carefully. We benchmark our resultsagainst Wisdom & Hubbard [58] using the identical valueof q rot =0.089195487.In Fig. 8, we show how discretization errors decay withincreasing number of spheroids. Overall the behavior issimilar to that of our more realistic Saturn interior modelin Fig. 6.We choose a acceleration factor of n int = 256 and per-formed a set of polytrope index 1 model calculations withincreasing precision. The number of explicitly treatedlayers, N L /n int were varied between 2 and 2 , whichbrought up the total number of layers to 131072 in ourlargest calculations, which is an increase of three ordersof magnitude compared to earlier CMS calculations. Weanalyze how our results improved with increasing layernumber and report the converged digits in Tab. II. Theagreement with the CLC predictions is excellent. All co-efficients J through J agree to 6, 7, or 8 significant dig-its, which is a better agreement than was reported in Wis-dom & Hubbard [58] where predictions from the CLC ap-proach and the non-accelerated CMS method were com-pared.0 Number of explicitly treated CMS layers, N tot L /n int -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 E rr o r i n e v e n J n = | J n − J ∞ n | ∆ J ∆ J ∆ J ∆ J ∆ J FIG. 8. Discretization error in the gravity harmonics of poly-trope index 1 planet models. The error of all gravity harmon-ics decays with increasing spheroid number, as we have seenfor the Saturn interior models in Fig. 6. All calculations forthis figure were performed with n int = 256.TABLE II. Gravity coefficients for the polytrope index 1planet models derived with the accelerated CMS (this work)and CLC [58] methods.Gravity coefficient CMS CLC10 × J × J − − × J × J − − × J × J − − × J × J − − × J × J − − H. Parameter Optimization
The primary goal of the model optimization is the gen-eration of Saturn interior models that reproduce the ob-served gravity harmonics. The agreement between mod-els and observations is typically expressed in some formof a χ function. Here we use, χ J = (cid:88) i =1 (cid:20) J observed2 i − J model2 i δJ observed2 i (cid:21) , (27)where δJ observed2 i are the 1- σ uncertainties in the obser-vations. Typically J is measured with much higher pre-cision than the higher order harmonics. To deal withthis imbalance, we find solutions that satisfy J observed2 i = J model2 i exactly by adjusting one model parameter like Z mol or Z met before χ J is evaluated. This optimization isperformed for converged CMS models that have reached hydrostatic balance and have matched the planet’s totalmass by adjusting the core mass.While Eq. 27 is certainly the most important optimiza-tion criterion, there are a number of other well motivatedconstraints to consider. For example, one would want toguide to the parameter optimization towards models withpressures P and P are close to the assumed immiscibil-ity curve in Fig. 3. From the assumed molecular andmetallic adiabats, we can infer the temperatures T and T that correspond to both pressures. For both pairs P - T and P - T , we find the closest points on the im-miscibility curve, P ∗ - T ∗ and P ∗ - T ∗ , that minimize thefollowing immiscibility penalty function, χ − He = (cid:88) i =1 C P (cid:12)(cid:12)(cid:12)(cid:12) P ∗ i − P i P i (cid:12)(cid:12)(cid:12)(cid:12) + C T (cid:12)(cid:12)(cid:12)(cid:12) T ∗ i − T i T i (cid:12)(cid:12)(cid:12)(cid:12) , (28)before we add the resulting minimum value to the total χ . C P and C T are weights that must be balanced withthose in other χ terms. We set C T /C P = 2. We chosenot to square the individual terms in Eq. 28 because,without an experimental confirmation of our immisci-bility curve, we do not want large deviations to enterquadratically. In Fig. 3, we shows some representativemodels to illustrate how much variation is in the P - T and P - T in our ensemble of models. Implicitly the χ − He term also introduces a penalty for metallic adia-bats that are too hot to be compatible with the assumedimmiscibility curve.Upon first introducing differential rotation into ourCMS models, we realized that a super-rotating equatorialjet improved the match to observed gravity harmonicsconsiderably. Furthermore, for l ≥ .
8, the inferred ro-tation profile was compatible with the wind speeds thatwere derived from tracking the clouds in Saturn’s atmo-sphere [11, 45]. From this point on, we favored modelsthat matched those observations by introducing the fol-lowing cloud penalty function, χ = C clouds (cid:88) k with l k > . | ω observed ( l k ) − ω model ( l k ) | , (29)where we sum over the knots, l k , in the outer region ofour rotation profile that often lead to good agreementwith the cloud tracking observations. C clouds plays therole of a weight.Finally we introduce one more penalty function, χ = C curvature (cid:88) all k [ ω (cid:48)(cid:48) model ( l k )] , (30)that favors smooth rotation profiles by penalizing largevalues in the second derivative of our rotational profile.We assume Y met ≥ Y mol and Z met ≥ Z mol , because weassume that helium rain can only lead to an enrichment ofthe metallic layer in helium and in heavy elements [56].We also constrain the helium abundance of the entireenvelope to match solar proportions.1 χ J FIG. 9. χ J deviations (see Eq. 27) between the calculatedand the observed gravity harmonics during the model opti-mization with the simplex algorithm. Only one model opti-mization (thick solid line) succeeded in converging to a statethat matched the spacecraft measurements well. We add Eqs. 27 through 30 to obtain one total χ function that we employ to optimize the model param-eters in Tab. IV. This turns out to be a very challeng-ing optimization problem, because many parameters arestrongly coupled and some optimization criteria are in-terdependent. We use the simplex algorithm [44] for theoptimization since it does not require any derivatives of χ with respect to the optimization parameters, whichare not available in analytical form. With this algorithm,it was very challenging to generate models that matchedobserved gravity data. In many cases, the algorithm getsstuck in a local minimum[ ? ]. Fig. 9 shows a couple of ex-amples of the χ J evolution during the simplex optimiza-tion. However, in 17 independent cases, the optimizationsucceeded and we were able to match the gravity har-monics within the uncertainties of the observations. Wesubsequently used these 17 solutions as starting pointsfor Markov chain Monte Carlo calculations in order tomap out the allowed parameter regions. We confirmedthat all 17 original solutions belong to the same param-eter region and one can go smoothly from one to theother. This provides strong evidence that the entire so-lution space is connected. I. Effects of an upper atmosphere
All CMS calculations presented so far, start from anoutermost spheroid with the pressure of 0.1 bar that wasanchored at the equatorial radius a . We had thereby ne-glected the effects of the tenuous upper atmosphere thatextends from the 0.1 bar level out into space. To studythe effects of this upper atmosphere quantitatively, weadded 64, 128, and 256 outer spheroids to CMS calcula-tions with 512, 1024, and 2048 layers, respectively. The TABLE III. Correction to the computed gravitional momentsdue to upper atmosphere that extends from 100 to 0.9 mbar. n × ∆ J n . × − − . × − . × − − . × −
10 +2 . × − − . × − − . × −
16 +4 . × − number of the additional spheroids was chosen such thatrange of pressure extended down to at least 1 mbar. Theoriginal outer spheroid is still associated with a pressureof 0.1 bar and remains anchored at the equatorial radius a . For all spheroids interior to this spheroid, we updatethe pressure according to Eq. 9 as we did before. How-ever, for all the additional exterior spheroids, we updatethe pressure with decreasing spheroid index according to, P (new) i = P (new) i +1 + ρ i ( U i +1 − U i ) . (31)For simplicity, we assume an isothermal upper atmo-sphere with a temperature set to the value at 0.1 bar.In all other respects the additional exterior spheroids aretreated in the same way as the interior spheroids. In prin-ciple, the temperature treatment could be made more re-alistic but, as we will show, the effects on the computedgravitational moments are negligible because there is solittle mass outside of the 0.1 bar level.Our extended CMS calculations converged to the samelevel of accuracy as they had previously. The pressure ofthe new outermost spheroid converged to 0.9 mbar. Thefractional mass outside of 0.1 bar level was found to beonly 7 . × − . This mass correction can also be in-terpreted as a change to the gravity coefficient J (seeEq. 14), which helps one to gauge the magnitude of thecorrection to the other gravity coefficients. In table III,we provide the differences in the gravitational momentsbetween our CMS calculation that included an extendedatmosphere to 0.9 mbar and our original calculations thatterminate at 0.1 bar. All values decay smoothly with in-creasing order n . For all the gravity coefficients that havebeen determined by the Cassini spacecraft, J through J , one finds that the correction due to the upper atmo-sphere is at least 18 times smaller than the uncertaintiesof the Cassini measurements [21]. For this reason, weconclude that our standard CMS calculations startingfrom the 0.1 bar level are sufficiently accurate for thisstudy.2
III. RESULTSA. Predictions for interior parameters
In Fig. 5, we plot the rotational profiles that haveemerged from our Monte Carlo calculations. Two promi-nent features are common to all models. There is a super-rotating equatorial jet in the equatorial region that ro-tates up to 4% faster than the deep interior. This be-havior is in agreement with the observed wind speedsfrom tracking the cloud motion in Saturn’s visible atmo-sphere [11, 45] and we have thus favored the sampling ofsuch models by introducing the term χ in Eq. 29.At a distance of approximately 50 000 km from the axisof rotation, our models require a sub-rotating region witha flow about 1% slower than in the deep interior. Thisfeature is not observed in the cloud motion at the surface,but is a common feature to all of our models that matchthe Cassini gravity harmonics. Both the super-rotatingequatoral jet, and the sub-rotating feature are presentregardless of the value we assume for the rotation periodof the deep anterior.In Fig. 11, we compare the predictions from ensem-bles of models that we generated with MC sampling fora range of core radii and rotation periods for the deepinterior. In panel (a), we plot the amount of heavy ele-ments in the envelope against the core mass. When onecompares models for the same core radius of r C = 0 . Z solar = 0 . Y mol values as low as 0.19are included. Y mol and Y met are tightly correlated in this figure because we assume the envelope overall has a solarhelium abundance.In Fig. 11d, we compare the heavy element abundancesin the molecular and metallic layers. Within the modelconstraint of Z met ≥ Z mol , a wide range of super-solarenrichments are predicted by our ensembles of models.There are plenty of models with Z met ≈ Z mol , which isin contrast to recent Jupiter models that required a dif-ferent amounts of enrichments in the two layers [55]. InFig. 11e and f, we plot the heavy element against thehelium abundances in the molecular and metallic layers,respectively. While both quantities are strongly corre-lated in the molecular layer, there appears to be muchmore flexibility in the metallic layer. One reason for thisis that a range of S met values are permitted in our mod-els while the entropy in the molecular layer is tied to thetemperature at the 1 bar level. In Fig. 11e, one can iden-tify a consistent trend for models with longer rotationperiods to predict larger values Y mol + Z mol and thus aslightly higher density for the molecular layer.The models with shorter rotation periods produce Y mol that are compatible with reanalyzed Voyager measure-ments of atmospheric helium, ∼ × solar [4], whilethe models with longer rotation rates require less de-pletion of helium in the outer envelope. Observationalconstraints on Z mol are uncertain; Fletcher et al. [7]observed atmospheric methane concentrations consistentwith ∼ × solar enrichment of carbon. There are no di-rect measurements of the abundance of atmospheric oxy-gen, the heavy element with the most significant con-tribution to the density and by consequence the grav-ity. Other heavy element ratios observations include bothmuch lower (N/H ∼ × solar) and higher S/H ∼ × so-lar), although these differences might reflect model de-pendence in determining the bulk elemental abundancefrom, or from measuring regions of the atmosphere thatour not well mixed [2].The models presented here predict values of both Z mol and Z met between 1–4 × solar for a uniform enrichmentof all heavy elements, which is lower than the observedenrichment in carbon. It is worth noting that in-situmeasurements of Jupiter’s atmosphere up to 22 bars bythe Galileo entry probe showed significant depletion inoxygen compared to carbon, but it is an outstandingquestion whether this accurately reflects the overall com-position of Jupiter’s molecular envelope. The heavy ele-ment content predicted by the models are also sensitiveto temperature of the adiabat. For Saturn, atmospherictemperature has never been measured in situ . So if S mol is higher than we expect, this tradeoff could account forhigher concentrations of heavy elements, without signifi-cantly affecting the other model predictions. B. Oblateness and Rotation Period
While the rotation period of Jupiter’s interior has beendetermined with high accuracy from magnetic field ob-3
Rotation period (s) O b l a t e n e ss , ( R e q − R p o l a r ) / R e q : : h : : h : : h : : h : : h : : h Occulationobservation,Lindal et al.(1985)CMS models,this workPreferred period
FIG. 10. Oblateness derived from CMS models with differentrotation periods compared to the radio occultation oblate-ness measurement by Lindal et al. [32], which determined anoblateness of 0.09822 ± ±
55 s. servations, the rotation period of Saturn’s deep interiorremains uncertain due to the remarkable alignment be-tween the dipole field and the axis of rotation. How-ever, the rotation period used in CMS calculations ofa planet significantly affects its shape. Saturn’s oblate-ness, ( R eq − R polar ) /R eq , has been measured with radiooccultation measurements of the Pioneer and
Voyager spacecraft [32]. Anderson and Schubert [1] constructedinterior models with uniform rotation that matched ob-served oblateness and pre-
Cassini gravity coefficients J , J , and J . They derived a rotation period of 10:32:44 h,which is significantly shorter than the system III periodof 10:39:22 h [6], as well as Cassini predictions of 10:45:45h and 10:47:06 h [12, 16].In Fig. 10, we compare models with differential rota-tion that we constructed for five different rotation periodsranging from 10:30:00 to 10:47:06 h. For all periods, itis possible to construct interior models with differentialrotation that match all even gravity coefficients. How-ever, the oblateness sensitively depends on rotation pe-riod that is assumed for the planet’s deep interior. InFig. 10, we compare the oblateness that derived from ourmodels with the radio occultation measurements by the
Pioneer and
Voyager spacecraft [32]. We find that rota-tion periods of 10:33:34 h ±
55 s are consistent with theseobservations, which represents a modest 1- σ increase overthe earlier determination of 10:32:44 h by Anderson and Schubert [1] who performed a similar analysis based oninterior models with uniform rotation. Thus, the 50 sdifference can primarily be attributed to effects of differ-ential rotation.Our determination of Saturn’s rotation rate isin remarkably good agreement with the value of10:33:38 h +112s − inferred from waves observed in Saturn’srings [36], even though the interior models for this anal-ysis were constructed without considering differential ro-tation.In Fig. 12, we compare predictions from models withour preferred rotation period of 10:33:34 h to those basedon the system III rotation period of 10:39:22 h. We stillpredict 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 rangeof Z met is considerably narrowed if the rotation periodis set to 10:33:34 h. Most Z met values now fall betweenvalues between 1.8 and 2.5 Z solar while in Fig. 11d, thesmaller and larger Z met values came from models withlonger and shorter rotation periods, now disfavored be-cause of the oblateness constraint. The Z mol values varybetween 1 and 3-fold Z solar as before.Finally, in Fig. 12c, Z mol and Y mol are now fairly tightlycorrelated when we assume a rotation period of 10:33:34h. These predictions can, in principle, be verified withremote observations or by an entry probe on a futuremissions. IV. CONCLUSIONS
We have presented an accelerated version of CMS thatallows us to construct planetary interior models withmany more layers than before and also enables construc-tion of ensembles of models using Monte Carlo methodsto efficiently optimize the parameters of individual mod-els. We have applied this accelerated CMS method toconstruct models for Saturn’s interiors with differentialrotation on cylinders, which permitted us to match theunexpectedly large values of the gravity harmonics J , J , and J that the Cassini spacecraft measured duringits Grand Finale orbits around Saturn. From our inte-rior models we infer that Saturn has a massive core of ∼ ±
55s for Saturn’s deep interior. [1] Anderson, J. D., & Schubert, G. 2007, Science (80-. ).,317, 1384[2] Atreya, S. K., Crida, A., Guillot, T., et al. 2019, in Sat-urn in the 21st Century, ed. K. H. Baines, F. M. Flasar, N. Krupp, & N. Stallard No. 1 (Cambdridge UniversityPress), 5–43[3] Cao, H., & Stevenson, D. J. 2017, J. Geophys. Res. Plan-ets, 122, 686 [4] Conrath, B. J., & Gautier, D. 2000, Icarus, 144, 124[5] Debras, F., & Chabrier, G. 2018, Astron. Astrophys.,609, doi:10.1051/0004-6361/201731682[6] Desch, M. D., & Kaiser, M. L. 1981, Geophys. Res. Lett.,8, 253[7] Fletcher, L. N., Orton, G. S., Teanby, N. A., Irwin, P. G.,& Bjoraker, G. L. 2009, Icarus, 199, 351. http://dx.doi.org/10.1016/j.icarus.2008.09.019 [8] Folkner, W. M., Iess, L., Anderson, J. D., et al. 2017,Geophys. Res. Lett., 44, 4694[9] Fortney, J. J., Helled, R., Nettelmann, N., et al. 2016,in Saturn in the 21st Century, ed. B. Kevin, M. Flasar,N. Krupp, & S. Thomas (Cambridge University Press),1–28. http://arxiv.org/abs/1609.06324 [10] Galanti, E., Kaspi, Y., & Tziperman, E. 2017, J. FluidMech., 810, 175[11] Garc´ıa-Melendo, E., P´erez-Hoyos, S., S´anchez-Lavega,A., & Hueso, R. 2011, Icarus, 215, 62[12] Giampieri, G., Dougherty, M. K., Smith, E. J., & Russell,C. T. 2006, Nature, 441, 62[13] Gudkova, T., & Zharkov, V. 1999, Planet. Space Sci.,47, 1201. http://linkinghub.elsevier.com/retrieve/pii/S0032063399000446 [14] Guillot, T. 1999, Planet. Space Sci., 47, 1183. http://linkinghub.elsevier.com/retrieve/pii/S0032063399000434 [15] Gurnett, D. A., Kurth, W. S., Hospodarsky, C. B., et al.2005, Science (80-. )., doi:10.1126/science.1105356[16] Helled, R., Galanti, E., & Kaspi, Y. 2015, Nature, 520,202. http://dx.doi.org/10.1038/nature14278%5Cn10.1038/nature14278 [17] Helled, R., & Guillot, T. 2013, Astrophys. J., 767, 113. http://stacks.iop.org/0004-637X/767/i=2/a=113?key=crossref.045d858be83734acdc0600277a318377 [18] Hubbard, W. 1982, Icarus, 52, 509. http://linkinghub.elsevier.com/retrieve/pii/0019103582900112 [19] Hubbard, W. B. 2012, Astrophys. J., 756, L15. http://stacks.iop.org/2041-8205/756/i=1/a=L15?key=crossref.34b95153bc3fdfb844cab51abbdf75d3 [20] —. 2013, Astrophys. J., 768, 43. http://stacks.iop.org/0004-637X/768/i=1/a=43?key=crossref.a31bd47c857111e805198695ba70780e [21] Iess, L., Militzer, B., Kaspi, Y., et al. 2019, Science,2965, eaat2965. [22] Jacobs, D. E. 1977, The State of the Art in NumericalAnalysis (Academic Press)[23] Jacobson, R. A., Antresian, P. G., Bordi, J. J., et al.2006, Astrophys. J., 132, 2520[24] Jeans, J. H. 2009, Problems of Cosmology and Stellar Dy-namics (Cambridge University Press). http://dx.doi.org/10.1017/CBO9780511694417 [25] Kaspi, Y. 2013, Geophys. Res. Lett., 40, 676[26] Kaspi, Y., & Galanti, E. 2016, Astrophys. J., 820, 91[27] Kaspi, Y., Guillot, T., Galanti, E., et al. 2017, Geophys.Res. Lett., 44, 5960[28] Kaspi, Y., Galanti, E., Hubbard, W., et al. 2018, Nature,555, doi:10.1038/nature25793[29] Kong, D., Liao, X., Zhang, K., & Schubert, G. 2013,Icarus, 226, 1425. http://linkinghub.elsevier.com/retrieve/pii/S0019103513003540 [30] Kuskov, O. L., & Kronrod, V. A. 2005, Icarus, 177, 550[31] Leconte, J., & Chabrier, G. 2013, Nat. Geosci., 6, 347. http://dx.doi.org/10.1038/ngeo1791 [32] Lindal, G. F., Sweetnam, D. N., & Eshelman, V. R. 1985,Astron. J., 90, 1136[33] Lindal, G. F., Wood, G. E., Levy, G. S., et al.1981, Journal of Geophysical Research: Space Physics,86, 8721. https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JA086iA10p08721 [34] Lodders, K. 2003, Astrophys. J., 591, 1220[35] —. 2010, in Astrophysics and Space Science Proceedings,ed. A. Goswami & B. E. Reddy (Berlin: Springer-Verlag),379–417[36] Mankovich, C., Marley, M. S., Fortney, J. J., &Movshovitz, N. 2019, The Astrophysical Journal, 871,1[37] Militzer, B. 2013, Phys. Rev. B, 87, 014202[38] Militzer, B., & Hubbard, W. B. 2013, Astrophys. J., 774,148[39] Morales, M. A., McMahon, J. M., Pierleonie, C., &Ceperley, D. M. 2013, Phys. Rev. Lett., 110, 065702[40] Morales, M. A., Pierleoni, C., Schwegler, E., & Ceperley,D. M. 2009, Proc. Nat. Acad. Sci., 106, 1324[41] Morton, K. W., & Mayers, D. F. 2005, Numerical solutionof partial differential equations: An introduction, 2ndedn. (Cambridge Univ. Press), arXiv:arXiv:1011.1669v3[42] Nettelmann, N., Fortney, J. J., Moore, K., &Mankovich, C. 2015, Mon. Not. R. Astron. Soc., 447,3422. http://mnras.oxfordjournals.org/cgi/doi/10.1093/mnras/stu2634 [43] Nettelmann, N., P¨ustow, R., & Redmer, R. 2013, Icarus,225, 548. http://linkinghub.elsevier.com/retrieve/pii/S0019103513001784 [44] Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flan-nery, B. P. 2001, Numerical Recipes in C++ (Cambridge,UK: Cambridge University Press)[45] Sanchez-Lavega, A., Rojas, J. F., & Sada, P. V. 2000,Icarus, 147, 405[46] Saumon, D., & Guillot, T. 2004, Astrophys. J., 609, 1170. http://stacks.iop.org/0004-637X/609/i=2/a=1170 [47] Seager, S., Kuchner, M. J., Hier-Majumder, C. A.,& Militzer, B. 2007, The Astrophysical Jour-nal, 669, 1279. http://adsabs.harvard.edu/cgi-bin/nph-data{_}query?bibcode=2007ApJ...669.1279S{&}link{_}type=ABSTRACT{%}5Cnpapers://0be24a46-325a-4116-a3c6-fd8a3b614472/Paper/p11167 [48] Soubiran, F., & Militzer, B. 2016, Astrophys. J., in press,http://arxiv.org/abs/1606.04162[49] Stevenson, D. J., & Salpeter, E. E. 1977, Astrophys JSuppl., 35, 239[50] Tassoul, J. L. 2015, Theory of Rotating Stars. (PSA-1), Princeton Series in Astrophysics (Princeton Uni-versity Press). https://books.google.com/books?id=nnJ9BgAAQBAJ [51] Vorberger, J., Tamblyn, I., Militzer, B., & Bonev, S.2007, Phys. Rev. B, 75, 024206[52] Wahl, S., Hubbard, W. B., & Militzer, B. 2016,Astrophys. J., 831, 14. http://stacks.iop.org/0004-637X/831/i=1/a=14?key=crossref.94e780e5342b88a796ca8fcc59534ff5 [53] —. 2017, Icarus, 282, 183[54] Wahl, S. M., Hubbard, W. B., & Militzer, B. 2017, Icarus,282, 183[55] Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017,Geophys. Res. Lett., 44, 4649 TABLE IV. Parameters and constraints in our Saturn models S mol Entropy H-He gas throughout the molecular layer.Constrained to match Saturns 1 bar temperature of 142.6 K [33] Y mol Helium mass fraction in the molecular layer. Constraint: Y mol ≤ Y solar = 0 . Z mol Mass fraction of heavy Z elements in the molecular layer. S met Entropy H-He gas throughout the metallic layer. S met ≥ S mol Y met Helium mass fraction in the metallic layer. Adjusted as function of Y mol to keep the overall composition of the envelope at Y solar . Z met Mass fraction of heavy Z elements in the metallic layer. Constraint: Z met ≥ Z mol . P Starting pressure of the helium rain layer (high pressure end of molecular layer) P Ending pressure of the helium rain layer (low pressure beginning of metallic layer)Core mass We assumed a compact core composed of heavy elementswith a sharp boundary to the metallic layer, with an equatorial radius, r C . ω ( l k ) Angular frequency for cylinder of radius, l k . Constraint: ω ( l k < .
7) = ω .TABLE V. Gravity coefficients predicted without acceleration scheme for different number of layers, N L . A representativeSaturn interior model with uniform rotation was used for this convergence analysis. N L J J J J J J J J
16 21058.747614 -1308.699233 119.572180 -13.513972 1.750222 -0.249135 0.037999 -0.00610832 17740.199991 -1047.173049 92.601757 -10.209559 1.294184 -0.180542 0.027002 -0.00425864 16789.968480 -978.984988 85.931473 -9.416858 1.186999 -0.164690 0.024500 -0.003843128 16552.657945 -962.339808 84.321088 -9.226705 1.161402 -0.160918 0.023907 -0.003745256 16493.667469 -958.249449 83.927467 -9.180331 1.155164 -0.159999 0.023762 -0.003721512 16478.115263 -957.156333 83.821482 -9.167812 1.153482 -0.159752 0.023724 -0.0037151024 16474.382140 -956.897102 83.796519 -9.164871 1.153087 -0.159694 0.023714 -0.0037132048 16473.437419 -956.831274 83.790168 -9.164122 1.152986 -0.159679 0.023712 -0.0037134096 16473.201201 -956.814816 83.788580 -9.163935 1.152961 -0.159675 0.023712 -0.0037138192 16473.142127 -956.810701 83.788183 -9.163888 1.152955 -0.159674 0.023711 -0.003713[56] Wilson, H. F., & Militzer, B. 2010, Phys. Rev. Lett., 104,121101[57] —. 2014, Astrophys. J., 973, 34[58] Wisdom, J., & Hubbard, W. B. 2016, Icarus, 267, 315. http://dx.doi.org/10.1016/j.icarus.2015.12.030 [59] Zharkov, V. N., & Trubitsyn, V. P. 1978, Physics of Plan-etary Interiors (Pachart), 388
ACKNOWLEDGMENTS
The authors acknowledge support from the NASA mis-sions
Cassini and
Juno . In part, the University of Cali-fornia supported this work through multicampus researchaward 00013725 for the Center for Frontiers in High En-ergy Density Science.6
TABLE VI. Gravity coefficients for the Saturn interior model in Tab. V predicted with acceleration factor n int = 16. The firstcolumn denotes the number of CMS layers that were treated explicitly and the second specifies the total layer number. Thelast row contains the extrapolated values for N tot L → ∞ . N tot L /n int N tot L J J J J J J J J ∞ ∞
15 16 17 18
Core mass (Earth masses) H e a vy e l e m e n t s i n e n v e l o p e ( E a r t h m a ss e s ) S o l a r e n v e l o p e × s o l a r e n v e l o p e × s o l a r e n v e l o p e × s o l a r e n v e l o p e (a) P = 10h30m00s, r C =0.200P = 10h32m44s, r C =0.200P = 10h45m45s, r C =0.200P = 10h47m06s, r C =0.200 P = 10h39m22s, r C =0.188P = 10h39m22s, r C =0.200P = 10h39m22s, r C =0.231 S met Y m e t + Z m e t (b)0.18 0.20 0.22 0.24 0.26 0.28 Y mol Y m e t Y mol =Y met =Y solar (c) 0.01 0.02 0.03 0.04 0.05 0.06 Z mol Z m e t (d)0.18 0.20 0.22 0.24 0.26 0.28 Y mol Z m o l Y solar ,Z solar Z mol =Z solar Z mol =2 × Z solar Z mol =3 × Z solar Z mol =4 × Z solar Z mol =5 × Z solar (e) 0.27 0.28 0.29 0.30 0.31 0.32 Y met Z m e t Z mol =Z solar Z mol =2 × Z solar Z mol =3 × Z solar Z mol =4 × Z solar Z mol = 5 × Z solar (f) FIG. 11. Comparison of parameters from seven sets of interior models including differential rotation. The assumed rotationperiods [6, 12, 15, 16] and fractional core radii are indicated by the color and symbol, as specified in the legend. r C = 0 . r C = 0 .
231 correspond to cores with iron-silicate and iron-silicate-ice compositions respectively. (a) The distribution ofheavy element mass between the core and envelope. (b) The variation of the mass fraction elements heavier than hydrogen withentropy in the inner, metallic envelope. (c) The variation of helium mass fraction in the molecular and metallic envelopes. (d)The variation of heavy element mass fraction in the molecular and metallic envelopes. (e) The tradeoff between heavy elementand helium mass fractions in the molecular envelope. (f) The tradeoff between heavy element and helium mass fractions in themetallic envelope. In panels (b),(c),(d),(e) and (f) solar values [35] are shown with a yellow star, corresponding to an assumedend-member case with no partitioning of helium through rain-out.
15 16 17 18
Core mass (Earth masses) H e a vy e l e m e n t s i n e n v e l o p e ( E a r t h m a ss e s ) S o l a r e n v e l o p e × s o l a r e n v e l o p e × s o l a r e n v e l o p e × s o l a r e n v e l o p e (a) P = 10h39m22s, r C =0.188P = 10h39m22s, r C =0.231 P = 10h33m34s, r C =0.188P = 10h33m34s, r C =0.231 Z mol Z m e t (b)0.18 0.20 0.22 0.24 0.26 0.28 Y mol Z m o l Y solar ,Z solar Z mol =Z solar Z mol =2 × Z solar Z mol =3 × Z solar Z mol =4 × Z solar Z mol =5 × Z solar (c) 0.27 0.28 0.29 0.30 0.31 0.32 Y met Z m e t Z mol =Z solar Z mol =2 × Z solar Z mol =3 × Z solar Z mol =4 × Z solar Z mol = 5 × Z solar (d) FIG. 12. Comparison of parameters from interior models with the preferred rotation rate of 10:33:34 h from this paper,compared to the System III rotation rate [6]. Core radii r C = 0 .
188 and r C = 0 ..