## Cassini's last look at Saturn's rings

During the final stages of the Cassini mission, the spacecraft flew between the planet and its rings, providing a new view on this spectacular system (see the Perspective by Ida). Setting the scene, Spilker reviews the numerous discoveries made using Cassini during the 13 years it spent orbiting Saturn. Iess *et al.* measured the gravitational pull on Cassini, separating the contributions from the planet and the rings. This allowed them to determine the interior structure of Saturn and the mass of its rings. Buratti *et al.* present observations of five small moons located in and around the rings. The moons each have distinctive shapes and compositions, owing to accretion of ring material. Tiscareno *et al.* observed the rings directly at close range, finding complex features sculpted by the gravitational interactions between moons and ring particles. Together, these results show that Saturn's rings are substantially younger than the planet itself and constrain models of their origin.

*Science*, this issue p. 1046, p. eaat2965, p. eaat2349, p. eaau1017; see also p. 1028

## Structured Abstract

### INTRODUCTION

The interior structure of Saturn, the depth of its winds, and the mass and age of its rings have been open questions in planetary science. The mass distribution inside a fluid and rapidly rotating planet like Saturn is largely driven by the ratio between centrifugal and gravitational forces. In static conditions, the planet should rotate uniformly and its gravity field should be axially and hemispherically symmetric and thus described by even zonal harmonics. However, optical tracking of clouds in the atmospheres of gas giants indicates that their rotation is not uniform. If the velocity field seen at cloud level extends deep into the interior, then there is a redistribution of mass that modifies the gravity field. Gravity measurements therefore constrain both Saturn’s deep interior and the depth of its winds. They also provide the mass of its rings, which dynamical and compositional dating methods have shown is related to their age.

### RATIONALE

In its final 22 orbits, the Cassini spacecraft passed between the rings and the atmosphere of Saturn. In six of these orbits, while the spacecraft was in free fall under the combined attraction of Saturn and its rings, radio tracking from an antenna on Earth was established to measure the radial velocity of the spacecraft with accuracies between 0.023 and 0.088 mm·s^{−1} at a 30 s integration time. Fitting these data with a dynamical model of the forces acting on the spacecraft allowed us to determine the separate signatures from each zonal harmonic and the ring mass.

### RESULTS

The estimated quantities were a zonal gravity model for Saturn and the gravity from a uniformly distributed mass spanning the A, B, and C rings. Additional small accelerations of unknown origin (possibly related to a time-variable gravity field) were required to obtain a good fit of the data. The determination of Saturn’s static zonal coefficients and ring mass does not depend on the assumed source of these small effects. We found that the measured values of *J*_{6}, *J*_{8}, and *J*_{10} are so large that they cannot be matched with interior models relying on uniform rotation and plausible compositions, but they are in agreement with interior models that assume deep differential rotation, extending from the equator into the interior up to distances of 0.7 to 0.8 Saturn radii from the spin axis. The equatorial outer layers of Saturn must rotate at an angular velocity that is 4% faster than that of the deep interior, whereas regions at higher latitudes must rotate 1 to 2% slower, regardless of the assumed rotation period. A thermal-wind approach shows that flow profiles that are similar in general character to the observed one yield solutions within the uncertainty of the gravity measurements when extended to a depth of ~9000 km, confirming that the flows are very deep and likely extend down to the levels where magnetic dissipation occurs.

We also found that the mass of the rings is 0.41 times that of the Saturnian moon Mimas, which is at the lower end of the expected values.

### CONCLUSION

Saturn’s gravity field measured by Cassini implies a strong and deep differential rotation, extending to a depth of ~9000 km. This differs from Jupiter, where winds are shallower (~3000 km). The gravity measurements are consistent with a mass of Saturn’s core of 15 to 18 Earth masses.

The low value of the ring mass suggests a scenario where the present rings of Saturn are young, probably just 10 million to 100 million years old, to be consistent with their pristine icy composition. Nevertheless, the rings may have evolved substantially since their formation and were perhaps once more massive than they are today. Models for a young ring system invoke the chance capture and tidal disruption of a comet or an icy outer Solar System body, suggesting that catastrophic events continued to occur in the Solar System long after its formation 4.6 billion years ago.

## Abstract

The interior structure of Saturn, the depth of its winds, and the mass and age of its rings constrain its formation and evolution. In the final phase of the Cassini mission, the spacecraft dived between the planet and its innermost ring, at altitudes of 2600 to 3900 kilometers above the cloud tops. During six of these crossings, a radio link with Earth was monitored to determine the gravitational field of the planet and the mass of its rings. We find that Saturn’s gravity deviates from theoretical expectations and requires differential rotation of the atmosphere extending to a depth of at least 9000 kilometers. The total mass of the rings is (1.54 ± 0.49) × 10^{19} kilograms (0.41 ± 0.13 times that of the moon Mimas), indicating that the rings may have formed 10^{7} to 10^{8} years ago.

The mass distribution inside a fluid and rapidly rotating planet like Saturn is largely driven by the ratio between centrifugal and gravitational forces. In the absence of internal dynamics, axial and hemispherical symmetry is expected, implying that when the gravitational potential is decomposed into spherical harmonics, only the even zonal harmonics appear (zonal harmonics are longitude-independent). Assuming hydrostatic equilibrium, interior models of gas giant planets indicate that the zonal coefficients *J*_{2}* _{k}* can be approximated by

*J*

_{2}

*=*

_{k}*a*, where

_{k}q^{k}*q*is the ratio of the centrifugal and gravitational acceleration at the equator (~0.16 for Saturn),

*k*is a positive integer, and the coefficients

*a*depend on the density profile inside the planet (

_{k}*1*).

Optical tracking of clouds indicates that dynamical phenomena operate on Saturn and Jupiter. The measured zonal (west-east) wind velocity field suggests a state of differential rotation (DR), whereby the angular velocity at any location depends on its distance from the axis of rotation and the depth along this axis (*2*, *3*). If the velocity field seen at the top of the atmosphere (conventionally defined as the level of 1 bar pressure) continues into the interior, then internal dynamics are expected to affect the gravitational field in two ways. First, the equipotential surfaces are perturbed symmetrically, redistributing mass in such a way that the even zonal coefficients deviate from the relation *J*_{2}* _{k}* ~

*q*(

^{k}*2*). Second, any north-south asymmetry in the velocity field leads to nonzero values of the odd zonal harmonics (

*4*). These theoretical expectations have been confirmed by the Juno mission at Jupiter (

*5*–

*7*), where gravity measurements showed that zonal winds are 2000 to 3000 km deep and suggest that the heavy-element core is diffuse (

*8*).

Gravity measurements at Saturn can be used to determine the mass of its rings, which is related to their age, by dynamical and compositional dating methods (*9*–*11*). Before the Grand Finale phase of the mission, the pericenter of Cassini’s orbit was always outside Saturn’s A ring, so that the gravitational effects of the rings could not be separated from those of the oblateness of the planet. During the Grand Finale, Cassini flew between the planet and the rings. This geometry effectively breaks the degeneracy between the even zonal field and the mass of the rings, providing a direct, dynamical estimate of the ring mass.

## Cassini gravity measurements

We determined Saturn’s gravitational field by reconstruction of Cassini’s trajectory during the Grand Finale using a coherent microwave link between Earth tracking stations and the spacecraft. Range-rate measurements were obtained from the Doppler shift of a carrier signal sent from the ground at 7.2 GHz (X band) and retransmitted back to Earth by Cassini’s onboard transponder at 8.4 GHz. An auxiliary downlink at Ka band (32.5 GHz) was also recorded.

In April 2017, Cassini was inserted into a series of inclined, highly eccentric orbits, grazing Saturn’s cloud tops at each pericenter (table S1). The orbit nodes were chosen such that the angle between the orbit normal and the direction to Earth is close to 90°. This edge-on condition provides the maximum projection of the spacecraft velocity along the line of sight, making it optimal for range-rate measurements and gravity estimations.

Of the 22 Grand Finale orbits (designated Rev 271 through Rev 293), six were selected for gravity measurements. Five orbits (Revs 273, 274, 278, 280, and 284) provided useful data (data from Rev 275 were lost because of a station configuration error). These orbits were selected to minimize neutral particle drag and maximize the spacecraft viewing period around closest approach (C/A).

We produced an orbital solution based on 5 data arcs, using Doppler observations with count times of 30 s, spanning a period of 24 to 36 hours enclosing each closest approach. Tracking data were acquired by the antennas of the three complexes (Goldstone, USA; Madrid, Spain; and Canberra, Australia) of NASA’s Deep Space Network (DSN), and two deep space antennas from the European Space Agency’s ESTRACK network, located in the southern hemisphere (Malargüe, Argentina, and New Norcia, Australia).

Two-way Doppler measurements at X band make up 93% of the dataset (*12*), with the addition of a few three-way passes to fill gaps during ground station handovers. Available X-band range observables are also included. All data with either uplink or downlink elevation below 15° were discarded to avoid systematic measurement errors due to imperfect calibration of path delays in Earth’s troposphere. Calibrations of the tropospheric and ionospheric path delays were provided by the DSN on the basis of pressure, humidity, and Global Positioning System data. Noise from solar plasma, which can strongly affect X-band radio links, was low owing to the large solar elongation angle (>142° on all six gravity orbits) (*13*). The data quality was statistically equivalent in X- and Ka-band data, with a root mean square Doppler noise between 0.020 and 0.088 mm·s^{−1} at a 30-s integration time (fig. S1 and table S1). For comparison, the Doppler signals due to the weakest measurable harmonics (*J*_{3}, *J*_{10}) and Saturn’s ring are 40 to 200 times larger than the average Doppler noise (fig. S2).

The X-band Doppler data were favored in our analysis because of the higher signal-to-noise ratio during the unavoidable ring occultation periods. The radio link was maintained during the occultations, except for a blockage of ~10 min on each orbit when the signal crossed the opaque core of the optically thick B ring (region B3). A second, longer blockage period due to ring occultations occurred on the outbound leg of the gravity orbits. Diffraction and near-forward scattering of the radio signal by ring particles caused an increase in the Doppler noise at X band by a factor of 2 to 3 during the occultation periods. Ka band is more sensitive to this effect and suffered repeated signal losses.

## Dynamical model

Our orbital fitting is based on the dynamical model previously adopted and tested in the determination of the gravity fields of Titan (*14*) and Enceladus (*15*), implemented in the Jet Propulsion Laboratory (JPL) navigation code MONTE (*16*). The model was extended to include Saturn’s gravitational parameter *GM* (where *G* is the gravitational constant and *M* is its mass), the zonal harmonic coefficients *J*_{2} to *J*_{20}, the tesseral (longitude-dependent) field of degree 2 (to account for possible nonprincipal axis rotation), and the mass of the rings. The truncation of the zonal field was set to twice the degree of the highest gravity harmonic whose central value is above the uncertainty (degree 10). The estimate of the even zonal harmonics used in our interpretation (degree ≤ 10) is insensitive to changes in the truncation. Although this model was adequate for fitting Juno gravity data at Jupiter (*5*), obtained in an orbital configuration similar to Cassini’s Grand Finale, it could not reduce Cassini’s Doppler residuals to the noise level. Instead, small stochastic accelerations were added to the model (see below) to produce signature-free residuals (fig. S1).

Although all the rings are included in the dynamical model, only rings A, B, and C can produce an acceleration potentially detectable by Cassini. The rings are assumed to be coplanar with Saturn’s equator, and each is assumed to have a constant surface mass density (the data are insensitive to radial variations of the density). Space- and ground-based measurements of ring occultations provide the determination of Saturn’s spin axis position and precession rate (*17*). Although Doppler data are sensitive to the orientation of Saturn’s equatorial plane, we adopt the prior determination because it is about ten times more accurate than that obtained from our orbital fitting.

Accelerations that are known to be large enough to produce noticeable signatures on range-rate data were accounted for. These include the point-mass gravitational accelerations from Saturn and its satellites (including the ring moons), computed from the JPL planetary and satellite ephemerides DE430, SAT389, and SAT393 (*18*), the acceleration from the Sun, the planets and satellites of the Solar System, and Saturn’s tidal response to its satellites (*12*).

The acquired Doppler data were combined in a multi-arc, weighted least-squares estimation filter. In the multi-arc approach, the entire time span of the observations is decomposed into shorter intervals, and a distinction is made between global and local estimated parameters. Global parameters are common to all arcs and estimated using all available observables. These include Saturn’s *GM*, *J*_{2} through *J*_{20}, *C*_{21}, *C*_{22}, *S*_{21}, *S*_{22}, and the masses of the A, B, and C rings. Although the estimates of the masses of the three rings are highly correlated, their sum is well determined. The ring masses were initially set to the current best estimates of their values from ring occultation data (*19*–*21*), with an a priori uncertainty of 100% (see table S2). The mass of the B ring had a large uncertainty of 10 Mimas masses (Saturn’s moon Mimas has a *GM* of 2.5026 km^{3}·s^{−2}). The sum of the masses of Saturn and its rings was constrained to be equal to the value (and uncertainty) estimated from satellite ephemerides (*18*).

Local parameters are those that belong only to a single arc and a different value was estimated for each arc. These included the Cassini position and velocity when the gravity observations began in each of the five orbits, at least 12 hours before transit at pericenter. The a priori uncertainties for position and velocity were set at 100 km and 1 m·s^{−1}, respectively.

In the estimation filter, parameters whose uncertainties are large enough to contribute to the final covariance matrix have been included as consider parameters. These include the Love number *k*_{22} (determining Saturn’s tidal response), Saturn’s pole direction, accelerations due to thermal emission from Cassini’s radioisotope thermoelectric generators (RTGs), and solar radiation pressure. The nominal value and a priori uncertainty of the Love number are taken from previous Cassini constraints (*22*), as were the pole position and precession rate (*17*). The anisotropic thermal acceleration from the RTGs was determined to 5% or better during the Cassini mission by the navigation team. The relative uncertainty associated with solar radiation pressure acceleration was set to 20%.

## Gravity determination

Our deterministic model, based on the geophysical expectations for the gravity field of a gas giant like Saturn, can adequately fit the Doppler data if each arc is analyzed separately. However, the same model cannot jointly fit all passes in a combined, multi-arc, gravity and orbital solution. The signatures in range-rate residuals are as large as 0.2 mm·s^{−1} over time scales of 20 to 60 min, corresponding to radial accelerations of the order of 10^{−7} m·s^{−2}. This value is an underestimate of the real unmodeled accelerations, as a large fraction of them are aliased in those associated with estimated parameters. The unmodeled accelerations acting on Cassini must be compensated for, to avoid biases in the estimates of the gravity harmonics and ring mass.

The missing accelerations could be due to longitudinally varying density anomalies resulting from wind dynamics or convection in Saturn’s deep interior (*23*). For a rocky planet, the corresponding gravitational field would be static and described by tesseral harmonics. However, this approach is not immediately applicable to a fluid planet like Saturn, because vortices move longitudinally with different speeds depending on latitude. The resulting gravity disturbances (caused by the nonzonal wind dynamics) would not be static in any reference frame for the 60-day duration of the Cassini gravity measurements. However, if the density anomalies are deep-seated, below the level of the zonal winds and in the region of uniform rotation, then a static tesseral field associated with convection might be maintained over time scales much longer than the duration of the experiment.

Although a high-degree tesseral field can fit the data, the required degree increases with the number of passes analyzed and depends on the assumed rotation period. Several measurements of Saturn’s rotation period have been made using different techniques. In our analysis we use four values—10 hours 32 min 45 s (*24*), 10 hours 39 min 22 s (also called System III) (*25*), 10 hours 45 min 45 s (*26*), and 10 hours 47 min 6 s (*27*)—leading to a required tesseral field of degree and order ranging from 8 to 12 (fig. S4). However, the determination of the tesseral field amplitudes is nonunique, because of the uncertainty in the underlying rotation rate of the deep interior (*12*).

Another potential solution is to assume a time-varying field, generated by acoustic oscillations in the planet (*28*). Fundamental mode oscillations (f-mode, with radial order equal to zero) have been detected by analysis of density waves in the C ring (*29*). These waves were also detected in Cassini ring stellar occultation data and are identified with sectoral modes which have an azimuthal order 2 ≤ *m* ≤ 10 (sectoral harmonics depend only on longitude). However, their absolute amplitudes are not well determined. In principle, an infinite number of modes can exist within Saturn, only a subset of which can drive waves in the rings. The inclusion of such normal modes in the dynamical model is possible, but the choice of the relevant modes is not unique. We found that different combinations of modes can fit the data but that it is impossible to identify unambiguously the dominant ones with the limited data available.

A more general approach, often adopted in orbit determinations in the presence of unknown or poorly modeled dynamics, is to assume empirical, random accelerations unrelated to any specific physical model. In the absence of geophysical evidence to pinpoint the root cause of the residual accelerations, we used this approach to obtain the baseline solution shown in Table 1. Our goal in this analysis was to disentangle the estimates of the zonal harmonics and the ring mass from the unmodeled accelerations. As an additional check on the robustness of the solution, we compared the baseline solution with solutions obtained using a tesseral field and multiple normal modes (*12*).

These random accelerations are assumed to act on Cassini’s orbit for a limited time span centered on pericenter and estimated in the orbital frame as local parameters. On each arc, these accelerations can act for equal-duration intervals at a constant value. The duration and the a priori uncertainty of the constant accelerations are adjustable parameters. In principle, it is desirable to optimize the number and duration of the intervals to avoid a degradation of the solution accuracy because of the over-parametrization of the problem. The a priori uncertainty must also be minimized to avoid aliasing part of the zonal gravity signal into the piecewise-constant accelerations. The minimum value of the a priori uncertainty (corresponding approximately to the largest allowed value of the random accelerations) represents the order of magnitude of the dynamical model’s incompleteness.

We found that random accelerations with an a priori uncertainty of 4 × 10^{−7} m·s^{−2}, acting for 10-min intervals within ±1 hour from pericenter, are able to remove the excess signatures in the Doppler residuals. The profile of the estimated random accelerations is reported in fig. S3.

The reference multi-arc solution is reported in Table 1 and Fig. 1. Before Cassini’s Grand Finale, the determination of Saturn’s even zonal coefficients *J*_{2}, *J*_{4}, and *J*_{6} was obtained from perturbations of the orbits of the moons, as well as from direct perturbations on Cassini itself (*18*). The data acquired during the Grand Finale orbits are consistent with these previous estimates; they also provide determination of the higher-degree gravity harmonics (*J*_{8} and *J*_{10}) and the odd harmonic coefficients *J*_{3} and *J*_{5}, the only odd harmonics whose values are larger than the associated uncertainties.

The data effectively constrain the sum of the masses of the A, B, and C rings, but the individual masses are poorly determined because of their large mutual correlations. However, the masses of the more-transparent A and C rings have been estimated from density waves seen in stellar occultations (*20*, *21*). We therefore estimated the masses of the A and C rings subject to those constraints (a priori values and uncertainties). Our final estimates of the A and C ring masses are essentially the same as the a priori values (see table S2). The uncertainty in the total ring mass (Table 1) is obtained from the correlation submatrix of the A, B, and C ring masses. We find a total ring *GM* value of 1.02 ± 0.41 km^{3}·s^{−2}, equivalent to 0.41 ± 0.13 Mimas masses (1σ uncertainties) (table S2).

The alternative tesseral field and normal modes approaches were compared with the piecewise-constant acceleration case to check the stability of our solution (*12*). The three solutions are in agreement with each other, as shown by the error ellipses in fig. S4 for pairs of gravity harmonics and in table S2 for the ring masses. We conclude that the estimate of the zonal harmonics is robust and largely independent of the dynamical model. This is especially true for the ring mass, which is consistent in all models and insensitive to the assumptions.

## Saturn interior models with uniform rotation

Although gravity measurements provide constraints on the interior of gas giants, every model using gravity data unavoidably suffers from uncertainties. We tested whether interior models based on uniform rotation can explain the measured gravity harmonics. We developed a suite of interior models based on the common assumption that the fluid in Saturn’s interior rotates uniformly like a solid body. All our models have a molecular and a metallic layer, each represented by an adiabat (a constant entropy curve in the temperature-pressure diagram) consistent with an equation of state (EOS) for hydrogen-helium mixtures determined from first-principles simulations (*30*, *31*) and characterized by an entropy, *S*, a helium mass fraction, *Y*, and a mass fraction of heavy elements, *Z*. We assume that helium rain occurs in Saturn’s interior wherever hydrogen and helium become immiscible, because hydrogen becomes metallic but helium does not (*32*). We thus introduce a helium rain layer that starts and ends at pressure-temperature conditions that are compatible with the hydrogen-helium immiscibility zone. We treat the helium rain layer as a smooth transition of the parameters across a range of pressures *P*_{1} to *P*_{2}, defined by the intersections of the adiabat with the immiscibility curve (*32*). These adiabatic profiles are thus described by six parameters (*S*_{mol}, *Y*_{mol}, *Z*_{mol}, *S*_{met}, *Y*_{met}, *Z*_{met}), where the subscript “mol” denotes the outer molecular envelope and “met” the inner metallic envelope. The models match the planet’s total mass and include dense cores of various sizes. Initially, we assumed a fractional radius of the core as *r*_{c} = 0.2 (*33*) and later refined the core radii by assuming a terrestrial iron-silicate composition (0.325:0.675) as well as a solar iron-silicate-water ice composition (0.1625:0.3375:0.5), which is consistent with Callisto’s interior (*34*). For these two compositions, we derived fractional core radii of *r*_{c} = 0.188 and 0.231, respectively, assuming hydrostatic equilibrium inside the core and previously published equations of state (*35*).

For each set of parameters, we construct a model for Saturn’s interior by performing a calculation with the concentric MacLaurin spheroid (CMS) method (*36*, *37*) to find a self-consistent shape and gravitational field for the planet. The simulation suite includes models with the four uniform rotation periods for Saturn’s interior observations (*24*–*27*). A range of interior distributions of helium were adopted, with a gradual gradient at a depth consistent with the phase diagram (*32*), and the planet-wide mass fraction of helium matching the solar fraction (*38*). A wide range of *Y*_{mol} was considered, bounded by the solar helium fraction *Y* = 0.274 and the lowest prediction *Y* = 0.055 from infrared measurements of Saturn’s atmosphere by Voyager (*39*). The range in entropy of the deep interior (*S*_{met}) is between the value in the outer envelope *S*_{mol} = 6.84 and the maximum value consistent with the predicted hydrogen-helium phase separation conditions, *S* ~ 7.20. Values of *Z*_{mol} up to five times solar fraction and fractional core radii up to 0.3 were considered. In each model, the core mass and deep heavy-element fraction (*Z*_{met}) were tuned to match the observed *J*_{2}, and an iterative approach was used to identify the narrower range of model parameters matching the observed *J*_{4}.

The resulting gravity harmonics are compared in Table 2 with the Cassini observations. Although a subset of the uniform rotation models is found to simultaneously match the observed *J*_{2} and *J*_{4}, the magnitudes of the observed values for *J*_{6}, *J*_{8}, and *J*_{10} are much larger than those predicted by the uniform rotation models (Fig. 1 and section 6 of the supplementary materials). Despite the wide range of physical parameters and interior structures considered, the tight correlation between *J*_{4} and the higher-order harmonics *J*_{6} to *J*_{10} precludes any possibility of finding a model that simultaneously matches all the observed harmonics. This result is unlike the gravity measurements of Jupiter from the Juno spacecraft. In this case, static models can be constructed, for which all even harmonics differ by less than 0.1 × 10^{−6} from the measurements (*8*), whereas Table 2 shows that the deviations for Saturn are much larger.

To test the robustness of this result, we added additional flexibility to the interior density distribution. Starting from a reasonable physical model, we introduced arbitrary density jumps of up to ±4% at four points in pressure that were chosen at random. We fit the models to the observed *J*_{2} value but report the full range of the higher *J*_{2}* _{k}* value without minimizing the discrepancy from the observations. These additional density modifications increased the range of all

*J*

_{2}

*values in our ensemble of models (Table 2), but a sizable discrepancy from the observed*

_{k}*J*

_{8}and

*J*

_{10}values remained. Figure 1, fig. S6, and Table 2 show the large and systematic deviations between the observations and models that assume uniform rotation. This demonstrates that specific assumptions regarding the EOS and conditions of helium rain cannot explain the inability of the models to match the higher-order harmonics.

The inability to reproduce the unusually large values of *J*_{6}, *J*_{8}, and *J*_{10} leads us to conclude that models with uniform rotation are ruled out by the gravity data. Instead, the observations require us to introduce DR (*2*) into our models for Saturn’s interior. We study its effects with two different techniques. First, we introduce DR on cylinders into the CMS method, which allows us to construct consistent interior models but requires constant rotation rate on cylinders that penetrate all the way through the planet. Second, we use the thermal-wind method (*3*) to derive the difference in the gravity signature between a uniformly and a differentially rotating body. The latter approach has the advantage that the DR profiles can have a finite depth and are not required to be north-south symmetric.

## Differential rotation on cylinders with the CMS method

In the CMS with DR approach, we simultaneously optimize interior parameters and DR profiles, ω(*l*), where ω is angular frequency and *l* is the cylindrical radius from the axis of rotation. We derive a centrifugal potential, *2*, *40*, *41*). Because this approach is based on potential theory, ω(*l*) can only depend on *l* and may not decay with depth along the cylinders. The benefit of the combined CMS + DR approach is that it represents a fully self-consistent velocity profile, interior density distribution, and gravitational field, rather than treating the winds as a correction to a model with uniform rotation.

We find that the assumption of DR can account for the unusually large magnitude of coefficients *J*_{6} to *J*_{10}. By assuming Saturn’s equatorial region rotates ~4% faster than the deep interior, agreement between models and data for all coefficients *J*_{2} to *J*_{10} can be achieved (Table 2). This assumption is also compatible with the equatorial part of the wind profile obtained from optical observations (*42*), showing strong eastward flow near the equator. We thus favored models that were in general agreement with the observed winds when we subsequently constructed an ensemble of interior models. For the planet’s deep interior, we assumed the same four rotation periods as above (*24*–*27*).

Starting from a parameter set chosen at random, we performed ~10^{4} independent model optimizations with the simplex algorithm. The 11 best models were used as starting points for subsequent Markov chain Monte Carlo (MCMC) calculations to explore the allowed parameter space more carefully. Models that matched the observed gravitational moments *J*_{2} to *J*_{10} are compared in Figs. 2 and 3. Figure 2 shows the rotational velocity profiles as a function of *l* and compares these to the corresponding north-south averaged surface winds derived from optical tracking of clouds (*42*, *43*).

From the best-fitting MCMC models, we selected quantities to elucidate the allowed region of parameter space (Fig. 3). Models that assume core masses between 15.0 and 18.2 Earth masses best match the gravity measurements (Fig. 3A). This range of core masses is broadly consistent with predictions of earlier models (*44*–*47*) built on previous gravity measurements and different assumptions for the EOS or structure of the deep interior. The predicted core masses increase if a smaller core radius is assumed, because the H-He mixture surrounding the core is exposed to higher pressures and becomes denser. The core masses also increase when a longer rotation period is assumed. As the total planet mass is constrained, the core mass anticorrelates with the fraction of heavy elements in the envelope (*Z*_{mol}). In comparison to the core (which we assume here consists entirely of elements heavier than helium), the mass of heavy elements in the envelope is relatively small, with a permitted range of 1.3 to 4.8 Earth masses. By contrast, analogous interior models of Jupiter predict a core mass slightly lower than Saturn’s, but with greater heavy-element enrichment in the deep envelope (*8*). Figure 3B shows the correlation between *Y*_{met} *+* *Z*_{met} and entropy of the metallic region *S*_{met.} An increase in entropy implies a higher temperature and a lower density. This density reduction can then be compensated by higher fractions of helium and heavy elements.

The helium fraction in the metallic layer and in the molecular layer are linearly related (Fig. 3C), as we require all models to have an overall solar helium abundance in the envelope. All models assumed that some helium rain has occurred, depleting the molecular layer and enriching the metallic layer in helium. We find that almost no helium sequestration is predicted for models that assume the longer rotation periods of 10 hours 45 min 45 s or 10 hours 47 min 6 s. This effectively removes one degree of freedom for the interior models and further constrains the range of other model parameters. For example, when we performed MCMC simulations with variable core radii for a period of 10 hours 47 min 6 s, no models with core radius larger than 0.21 and core masses larger than 17.2 Earth masses were obtained. Figure 3D shows that the heavy-element fraction in the molecular envelope may vary between 1 and 3.5 times the solar value. The heavy *Z* fraction is anticorrelated with the helium fraction, which implies that there is only limited capacity for elements heavier than hydrogen. Longer rotation periods imply that there is a greater capacity for such elements.

All models that match the even gravity harmonics (*J*_{2} to *J*_{10}) show a rapid decrease in angular velocity from *l* = 1.0 to *l* ~ 0.9 (Fig. 2). For all rotation periods under consideration, our CMS + DR models also require a minimum region near *l* ~ 0.83 where the fluid rotates slower than in the deep interior. This radius corresponds to a depth of ~10,000 km from the surface at the equator. By construction, the DR profiles converge to the assumed rotation rate for small *l*. The CMS + DR models are unable to reproduce the wind profile for small *l* (high latitudes) because we assume DR on cylinders.

## Differential rotation with finite depth flows

An alternative approach, the thermal-wind method (*3*), can incorporate wind profiles with a finite depth. In this approach, the flow is not limited to full cylindrical symmetry and the flows are not required to be equal in the Northern and Southern hemispheres. This allows the model to account for both the even and odd gravity harmonics; the latter reflect a north-south asymmetry in the gravity field (*4*). It has been shown that for the barotropic case (flows limited to cylindrical symmetry), the thermal-wind and the CMS method with DR produce consistent results (*41*, *48*). We apply the thermal-wind method using an adjoint inverse model (*49*), where the measured gravity harmonics in Table 1 are used to identify the flow structure that best fits the data. Although it has these advantages, the thermal-wind method can only account for the dynamical part of the gravity spectrum (attributable to the flows), and therefore relies on having prior knowledge of a reference interior model for the even gravity harmonics. The predicted gravity values thus depend on the reference model (Table 2).

For Jupiter, the extension of the cloud-level flow into the interior results in a gravity spectrum that matches the observations (*6*) and the background density profile. The case of Saturn is more complicated. We find that regardless of the vertical structure of the zonal flow, with the background density profiles corresponding to the 11 models discussed in the previous section, extending the observed cloud-level flow (*43*) into the interior results in dynamical gravity harmonics (Δ*J _{n}*) that are too small by at least a factor of two. For example, the gravity signature of the observed surface winds extended along cylinders to large depths gives roughly Δ

*J*

_{8}= −1.5 × 10

^{−6}and Δ

*J*

_{10}= 1 × 10

^{−6}(

*4*,

*50*). The measured

*J*

_{8}and

*J*

_{10}are −14.6 × 10

^{−6}and 4.7 × 10

^{−6}, respectively, whereas the uniform rotation contribution is, at most, −9 × 10

^{−6}and 1 × 10

^{−6}(Table 2). This implies that the observed cloud-level flow cannot account for the difference between the measurements and the uniform rotation contribution regardless of how it extends into the interior. As seen in Fig. 2, the only way to match the data is to alter the latitudinal profile of the zonal wind, by assuming that, because of other processes, the bulk of the sub-cloud-level flow does not resemble that observed at cloud level.

Therefore, to match the Cassini gravity measurements, another degree of freedom is added to the thermal-wind gravity inversion, whereby the meridional profile of the flow is allowed to vary from the observed winds in addition to the flow depth. We then optimize for the depth and meridional structure of the zonal flow so that the calculated gravity harmonics match the measured values. To compare the results to those of the hemispherically symmetric CMS method, and because, unlike Jupiter (*6*), the odd harmonics are small and contribute little to the flow structure, we use only the even gravity harmonics, focusing on those which are most strongly affected by the flows, namely *J*_{8} and *J*_{10.} Because of the uncertainty resulting from the uniform rotation reference model, and the nonuniqueness of the problem, there is little merit in optimizing for the vertical profile of the flow, as can be done for Jupiter (*6*), so the vertical decay is approximated by a radial hyperbolic tangent profile with a width of 500 km and is then optimized for its depth.

The reconstructed meridional profile of the zonal flow at the planet’s surface is shown in Fig. 4A, achieved with a vertical profile depth of 9363 ± 357 km. With these parameters, we are able to match the measured values of *J*_{8} and *J*_{10}, once the uniform rotation contribution in Table 2 is subtracted. Both gravity harmonics are within the uncertainties of the measured values (Table 3). The Δ*J* we optimize for are calculated from the difference between the measurements and the 11 preferred models run without DR and using the Voyager rotation period. Although the reconstructed zonal wind profile varies from the observed one, it still retains its main characteristics (Fig. 4A). The uncertainty of the wind depth is obtained from the variance in the 11 solutions of the internal models, the variance in the internal models and thermal-wind calculation for different rotation rates between 10 hours 32 min and 10 hours 47 min, and using different internal models, by statistically varying the depth, and finally searching for solutions that satisfy the condition that both harmonics are within the uncertainties. In all cases, the reported depth is the inflection point of the hyperbolic tangent vertical profile. The exact meridional structure of the zonal winds is not uniquely determined, meaning that other profiles can give results within the measurement uncertainties.

Regardless of the exact meridional profile, all vertical flow profiles are constrained to contain a very deep flow of ~9000 km. This depth, corresponding to 15% of Saturn’s radius, matches the predictions from magnetohydrodynamic theories, suggesting that the flow should extend down to the levels of magnetic dissipation (*51*, *52*). It has been suggested that a higher-order expansion to the vorticity balance, beyond a thermal wind, is necessary to reproduce the dynamical gravity harmonics (*53*). To compare those solutions [known as the thermal-gravity equation, (*53*)], we calculate the gravity harmonics using the thermal-gravity equation for a similar wind profile. The comparison (Table 3) shows that the solutions from both methods are within the observational uncertainties, confirming that the thermal-wind solution is the leading-order vorticity balance (*48*). The thermal-gravity solution was obtained by solving the full integro-differential equation on a sphere (*54*).

The CMS solutions with DR on full cylinders show a flow signature only for distances from the axis of rotation *l* ≥ 0.6 (Fig. 2). Using the thermal-wind inversion model, we can examine the nature of the optimized flow obtained under a similar restriction. In Fig. 4B, we show that taking the same approach with the thermal-wind model, restricting the flow to be equatorward only for *l* ≥ 0.6 (latitudes ≤ 60°), gives similar results. The reconstructed flow in this case is not as close to the observed flow as in the unrestricted case (Fig. 4A) but still gives results within 1σ of the deviation uncertainty (Table 3). We also find a large flow depth of 8832 ± 295 km in this solution. The deviations between the two flow models in Fig. 4 illustrate the nonuniqueness in the resulting flows when gravity harmonics Δ*J*_{8} and Δ*J*_{10} are modeled. However, all our thermal-wind flow solutions show strong westward flow near latitude 35°, which is in agreement with the westward flowing region calculated with the CMS model in Fig. 2. These results were calculated for a rotation period of 10 hours 39 min 22 s; because the high-degree harmonics are weakly affected by the rotation rate (*55*), we expect consistent results if a different period is assumed.

We find therefore that with the thermal-wind approach, extending the exact observed flows into the interior in a simple manner does not lead to an exact match of measured gravity values. However, flow profiles that are similar in general character to the observed ones (Fig. 4) yield solutions within the deviation uncertainty when extended to a depth of ~9000 km, indicating that the flows are very deep and likely extend down to the levels of magnetic dissipation.

The agreement between the CMS and thermal-wind solutions, based on two substantially different approaches, indicate that the interpretation of the Cassini gravity data is robust. Regardless of the exact interior model or vertical decay profile chosen, the flows must extend deep below the surface, to ~15% of the planet’s radius.

## Mass and age of Saturn’s rings

We now consider the mass of the main rings (i.e., the A, B, and C rings) and how this quantity can constrain theoretical models for the age and origin of the rings. (The masses of the diffuse D, E, F, and G rings are expected to be negligible.) Before the Voyager flybys in 1980 and 1981, it was known that the cross section–weighted average ring particle radius was at least 10 cm, as required to account for the rings’ high radar reflectivity (*56*). Voyager results provided estimates of the rings’ mass, on the basis of local surface mass densities determined from the wavelengths of 13 density waves in the A and B rings and an optical depth profile derived from a stellar occultation (*57*). The estimated total ring mass was 2.8 × 10^{19} kg, or 0.75 Mimas masses. However, it was subsequently argued that substantially more mass might be hidden in the opaque parts of the B ring (*58*).

Cassini observations of density waves in the rings have led to additional local surface density estimates of the A (*20*), C (*29*), and B rings (*19*). The B-ring density is higher than the others, with a mean value of ~600 kg·m^{−2}, but not higher by as much as its large optical depth would suggest. Combining these estimates, we find a likely total ring mass of *GM* = 1.01 km^{3}·s^{−2}, or 0.40 Mimas masses. However, the presence of self-gravity wakes in the A and B rings has led to suggestions that substantial amounts of material may be present but does not contribute to density wave estimates. Numerical simulations of these structures (*58*) suggest wake surface densities as high as 5000 kg·m^{−2} and an upper limit to the total ring mass of 9.7 × 10^{19} kg, or ~2.5 Mimas masses.

Our estimate of *GM* = 0.58 ± 0.48 km^{3}·s^{−2} for the B ring, when combined with the previously determined *GM* products of the A (0.38 km^{3}·s^{−2}) and C (0.06 km^{3}·s^{−2}) rings, leads to a total ring *GM* of 1.02 ± 0.43 km^{3}·s^{−2}, equivalent to 0.41 ± 0.13 Mimas masses (table S2). This estimate is somewhat smaller than the Voyager value but consistent with those inferred from Cassini observations of density waves.

The ring mass has implications for the age of the rings (*59*). Traditional age estimates of Saturn’s ring-satellite system fall into three categories: dynamical estimates based on the rate of recession of the small satellites from the rings owing to gravitational torques and the back-reaction on the A ring (*10*, *60*); structural evolution time scales based on the evolution of unconfined edges, such as the inner edges of the A and B rings (*11*); and compositional time scales based on the assumption that the rings were formed as pure water ice and have subsequently been steadily darkened by the infall of interplanetary debris (*9*). Most of these time scale estimates depend, directly or indirectly, on the mass of the rings; both dynamical and compositional considerations suggest that low-mass rings are likely to be young, and both approaches yield evolutionary ages ~10^{8} years for the A and B rings, assuming the Voyager measurements of ring masses and interplanetary impact fluxes (*59*). Low-mass rings therefore pose a challenge for models in which the ring system is assumed to have formed in the early Solar System.

The Cassini results have sharpened these arguments. Analysis of the main rings’ non-icy fraction derived from their passive microwave emission, using the pre-Cassini interplanetary flux estimates and a minimum-mass B ring, led to estimates of the ages of the A and B rings of 80 to 150 million years (Ma) and 30 to 100 Ma, respectively (*61*). Our revised mass estimate for the B ring would increase the latter estimate by ~25%. Cassini dust measurements have also provided a refined estimate of the interplanetary dust flux at the rings (*62*), indicating that the interplanetary dust impact flux on the rings is higher by almost a factor of 10 compared to the Voyager estimates, because of an improved understanding of the gravitational focusing by Saturn. These new dust fluxes would reduce the age estimates for the A and B rings.

As the rings evolve viscously, they may spread radially beyond the Roche limit and therefore gradually lose mass to form new generations of small satellites, which would then move rapidly away from the rings (*63*). In this scenario, the rings would have been more massive in the past, with proportionally longer evolutionary time scales. Current models of viscous evolution (*64*) predict that such a ring would approach an asymptotic mass of ~1.5 × 10^{19} kg, or 0.40 Mimas masses, after 5 billion years, close to our estimate of 0.41 ± 0.13 Mimas masses. Taken at face value, and subject to the limitations of the simulations (which ignore the gravitational torques from the newly formed moons), this suggests that the rings may be older and formed with more mass than they have today. However, the simulations also show that a massive, primordial ring initially loses mass very quickly before settling into a long period of progressively slower evolution (*64*). An old high-mass ring would thus quickly evolve to a mass not much greater than that measured today and would therefore still be subjected to rapid ballistic and compositional evolution.

On balance, we favor a scenario whereby the present rings of Saturn are relatively young, at least compared with the planet itself, although they may have evolved substantially in the past 10^{7} to 10^{8} years and were perhaps once more massive than they are today. Our data do not indicate how the ring system formed within such a recent period. Models of that process invoke the chance capture and tidal disruption of a comet or Centaur (*65*, *66*). Alternatively, rings may arise from the catastrophic disruption of an earlier population of midsized icy satellites ~100 Ma ago (*67*), although whether the resulting debris can migrate into the ring region before it reaccretes into new satellites is uncertain (*68*). Regardless of how the rings formed, the ring mass we derived from Cassini gravity data indicates a recent origin for Saturn’s ring system, which we consider a fitting way to end Cassini’s mission.

## Supplementary Materials

science.sciencemag.org/content/364/6445/eaat2965/suppl/DC1

Materials and Methods

Figs. S1 to S6

Tables S1 and S2

This is an article distributed under the terms of the Science Journals Default License.

## References and Notes

**Acknowledgments:**We thank J. W. Armstrong, J. Cuzzi, M. Di Benedetto, S. Finocchiaro, J. Fortney, R. French, L. Gomez-Casajus, M. Gregnanin, M. Hedman, R. A. Jacobson, E. Marouf, N. Movshovitz, and L. Spilker for their suggestions and contributions. We are grateful to the Cassini Radio Science operations team members (E. Barbinis, D. Kahan, J. Gao, C. Lee) for supporting the experiments described in this paper. We also thank the DSN and ESA personnel involved in the acquisition of the data. L.I. thanks the late Bruno Bertotti for his mentorship and contributions to the Cassini mission.

**Funding:**Support for this work was provided by the Italian Space Agency (L.I., D.D., P.R., M.J.M., P.T., and M.Z.), NASA’s CDAP program (B.M., W.H., and S.W.), the Cassini Project (P.N. and A.A.), and the Israeli Space Agency (Y.K. and E.G.). All authors acknowledge support from the Cassini Project.

**Author contributions:**L.I. led the experiment and supervised the data analysis. D.D. and P.R. carried out the gravity data analysis. B.M., W.H., and S.W. contributed interior models and the differential rotation with CMS method. Y.K. and E.G. contributed the models on differential rotation with finite depth flows. P.N. contributed the ring mass interpretation. A.A. coordinated the experiment operations and data acquisition. M.J.M., P.T., and M.Z. helped with the data analysis.

**Competing interests:**The authors declare no competing interests.

**Data and materials availability:**The Cassini data used in this paper is available through the NASA Planetary Data System at https://atmos.nmsu.edu/pdsd/archive/data/co-s-rss-1-sagr18-v10/ and https://atmos.nmsu.edu/pdsd/archive/data/co-s-rss-1-sagr19-v10/. The MONTE space navigation code was obtained through a license agreement between NASA and the Italian Space Agency; the terms do not permit redistribution. MONTE licenses may be requested at https://montepy.jpl.nasa.gov/. Our output gravitational field model is listed in Table 1. Executable code for running the interior models is available at https://github.com/smwahl/CMS_Saturn (CMS method) and https://github.com/egalanti/TW_Saturn (thermal-wind method, requires MATLAB).