# Saturn’s Seismic Rotation Revisited

Christopher R. Mankovich  Division of Geological and Planetary Sciences,
Mailcode 150-21, California Institute of Technology, Pasadena, CA 91125, USA;
<EMAIL_ADDRESS>Janosz W. Dewberry  TAPIR, Walter Burke Institute for
Theoretical Physics, Mailcode 350-17, California Institute of Technology,
Pasadena, CA 91125, USA  Canadian Institute for Theoretical Astrophysics, 60
St. George Street, Toronto, ON M5S 3H8, Canada  Jim Fuller  TAPIR, Walter
Burke Institute for Theoretical Physics, Mailcode 350-17, California Institute
of Technology, Pasadena, CA 91125, USA

###### Abstract

Normal mode seismology is a promising means of measuring rotation in gas giant
interiors, and ring seismology presents a singular opportunity to do so at
Saturn. We calculate Saturn’s normal modes of oscillation and zonal gravity
field, using nonperturbative methods for normal modes in the rigidly rotating
approximation, and perturbative methods for the shifts that Saturn’s deep
winds induce in the mode frequencies and zonal gravity harmonics. The latter
are calculated by solving the thermo-gravitational wind equation in an oblate
geometry. Comparing many such models to gravity data and the frequencies of
ring patterns excited by Saturn normal modes, we use statistical methods to
estimate that Saturn’s cloud-level winds extend inward along cylinders before
decaying at a depth $0.125\textrm{--}0.138$ times Saturn’s equatorial radius,
or $7,530\textrm{--}8,320$ km, consistent with analyses of Cassini gravity and
magnetic field data. The seismology is especially useful for pinning down
Saturn’s poorly constrained deep rotation period, which we estimate at
$2\pi/\Omega_{\rm S}=634.7$ min (median) with a 5/95% quantile range
$633.8-635.5$ min. Outstanding residuals in mode frequencies at low angular
degree suggest a more complicated deep interior than has been considered to
date. Smaller but still significant residuals at high angular degree also show
that our picture for the thermal, composition, and/or rotation profile in
Saturn’s envelope is not yet complete.

##  1 Introduction

Precision gravity science from the Juno spacecraft at Jupiter and Cassini at
Saturn has provided evidence for deep differential flows in the atmospheres of
both planets (Kaspi et al., 2018; Guillot et al., 2018; Iess et al., 2019;
Galanti et al., 2019; Militzer et al., 2019) and has stimulated new attention
toward the heavy element gradients that may be a common feature of gas giant
interiors (Wahl et al., 2017; Debras & Chabrier, 2019; Stevenson, 2020;
Nettelmann et al., 2021; Miguel et al., 2022; Militzer et al., 2022). But much
is still unclear about these planets’ internal structures. In the case of
Saturn, the near-perfect alignment between the rotation axis and the magnetic
dipole axis (Cao et al., 2020) obscures the spin of Saturn’s deep interior, a
key input for structure models. Hence, alternative observational means of
disentangling Saturn’s rotation and composition structure are highly sought
after.

Parallel to the gravity science, Juno and Cassini have also renewed interest
in giant planet seismology. Anomalous spacecraft accelerations appear to
support the idea of acoustic (p-)mode seismicity in Jupiter (Durante et al.,
2022) and Saturn (Markham et al., 2020), although this offers no immediate
quantitative constraints on planetary structure. On the other hand, the unique
data set of ring seismology at Saturn has proven useful for constraining a
giant planet’s interior in new ways. Spurred on by observations by Rosen et
al. (1991) and Baillié et al. (2011) and detailed predictions laid out by
Marley (1991) and Marley & Porco (1993), a growing corpus of Cassini
detections of Saturn-associated ring waves (Hedman & Nicholson, 2013, 2014;
French et al., 2016, 2019; Hedman et al., 2019; French et al., 2021; Hedman et
al., 2022) has provided precise measurements of the frequencies of 34 ring
waves that can each be traced back to a distinct nonradial normal mode of
oscillation inside Saturn.

Most of these waves take the form of spiral density waves at outer Lindblad
resonances (OLRs) or spiral bending waves at outer vertical resonances
(OVRs)111In addition to these 34, French et al. (2021) report two Saturn-
associated inner Lindblad resonances (ILRs), but these stem from the same
Saturn modes responsible for known OLRs, and hence do not add new frequency
constraints.. The corresponding Saturn modes can typically222 Degree $\ell$
and radial order $n$ are not good quantum numbers for modes strongly affected
by rotation or mode-mode degeneracies. The Saturn modes at hand take the shape
of slightly distorted spherical harmonics and hence can still be assigned a
single $\ell$ describing the _dominant_ angular degree. Radial order $n$ can
lose its meaning entirely in cases involving degeneracy, e.g., in rosette
modes, which are not the focus of this work.  be categorized by their
azimuthal order $m$, dominant angular degree $\ell$, and number of radial
nodes $n$. The known Saturn-associated detections are prograde modes spanning
azimuthal orders $m=2-14$, and comparison to models (Mankovich et al., 2019,
“M19”) makes it clear that these data correspond to Saturn modes with
$\ell-m=0-8$ (French et al., 2021). (Because we deal here with exclusively
prograde patterns, we use the convention that $m>0$ denotes prograde planet
modes and ring patterns.) Most of the detected modes are the fundamental
(f-)modes with $n=0$, effectively surface gravity modes that are confined
increasingly close to the planet’s surface with increasing $\ell$. M19 found
the f-modes to be highly sensitive to Saturn’s rotation, using the 14 f-mode
frequencies that were known at the time to make the first seismic estimate of
Saturn’s bulk rotation period.

Beyond the f-modes, the observations have also uncovered internal gravity
(g-)modes with $n\geq 1$ at the lowest azimuthal orders $m=2,3$, where they
can induce finely separated mode multiplets via degenerate mixing with f-modes
(Fuller, 2014), or simply appear as isolated g-modes (Mankovich & Fuller,
2021, “MF21”). This striking result of ring seismology directly implies that
not all of Saturn’s interior is convective, as the classical picture would
have it. MF21 identified an instance of an isolated g-mode in the $m=2$
density wave W76.44, ascribing it to the
${}_{\,\ell}^{m}\cdot_{n}=\,_{2}^{2}g_{1}$ mode and using its frequency to
constrain the extent of the stably stratified g-mode cavity associated with
Saturn’s diffuse core. Other isolated low-order g-modes may be observable
through $m\lesssim 5$. Dewberry et al. (2021, “D21”) addressed the phenomenon
of mode mixing and fine-scale multiplicity in detail, considering the
additional influence of the differential rotation implied by the Cassini
gravity results. They found that Saturn’s deep winds affect the frequencies of
f-modes, which in turn can drive mixing with other modes by bringing their
frequencies toward resonance.

Critically, while Fuller (2014), M19, and MF21 used perturbation theory to
account for Saturn’s rapid rotation in the oscillation calculations, D21
employed a nonperturbative method for solving for modes in a two-dimensional
oblate geometry (see, e.g., Reese et al., 2006; Ouazzani et al., 2012; Xu &
Lai, 2017), the first time such a calculation has been applied to Saturn ring
seismology. This technical leap led to the identification of a host of
‘rosette’ modes, modes that are essentially coupled sequences of g-modes
strongly influenced by rotation, named for their rosette-shaped eigenfunctions
(see D21, Figures 3-5; see also Ballot et al. 2012 and Takata & Saio 2013).

More germane to the question of Saturn’s rotation state is the major accuracy
advantage of nonperturbative methods compared to the approximate nature of
perturbative methods: The latter (e.g., M19) have generally included effects
to second order in the ratio $\Omega_{\rm S}/\Omega_{\rm dyn}$ where
$\Omega_{\rm S}$ is spin frequency and $\Omega_{\rm dyn}\equiv(GM_{\rm
S}/R_{\rm eq}^{3})^{1/2}$ is Saturn’s equatorial dynamical frequency. For
Saturn, $\Omega_{\rm S}/\Omega_{\rm dyn}\sim 0.4$ so that mode frequencies
from second-order theories are uncertain at the level of $(\Omega_{\rm
S}/\Omega_{\rm dyn})^{3}\sim 5\%$, a major limitation considering that ring
observations yield frequencies to precisions better than $0.01\%$ (French et
al., 2021). A realistic model certainly must also abandon the rigid rotation
approximation and address Saturn’s deep winds, which at
$|\delta\Omega/\Omega_{\rm S}|\sim 1\%$ may modify mode frequencies at the
level of $\delta\Omega/\Omega_{\rm dyn}\sim 0.4\%$. In contrast to the rapid
bulk rotation, the rather weak differential rotation lends itself well to a
perturbative treatment in mode calculations, an approximation whose validity
Dewberry et al. (2022) quantified by comparing with fully nonperturbative
calculations for differentially rotating models.

These theoretical tools and the clearer understanding of Saturn’s deep wind
structure warrant a reappraisal of the seismic rotation rate inferred by
Mankovich et al. (2019). Here we use the zonal gravity field and frequencies
of f-modes observed in the rings to jointly constrain the flow profile in
Saturn’s envelope, making the common assumption that the flow at depth follows
the alternating eastward and westward zonal jet streams observed at the cloud
level, with a simple decay prescription as a function of depth. From Markov
chain Monte Carlo (MCMC) analysis we find Saturn’s wind depth and deep spin to
both be consistent with previous estimates, and use precision ring seismology
through a nonperturbative numerical method to place a more precise constraint
on the deep spin than has been possible with approximate perturbative methods.

Section 2 describes our interior model and shape calculation, our method for
computing modes in the presence of rapid rotation and their frequency shifts
arising from differential rotation, the wind-induced contributions to zonal
gravity, and finally our MCMC sampling procedure. We present our findings in
Section 3, discuss in Section 4, and summarize in Section 5.

##  2 Methods

In brief our method is as follows: we

  1. 1.

Build rigidly rotating models for Saturn’s oblate background structure and
even zonal gravity field $J_{2n}$.

  2. 2.

Calculate each model’s adiabatic oscillation modes (and frequencies
$\sigma_{i}$) including non-perturbative treatment of rotational forces,
considering many azimuthal orders $m$.

  3. 3.

Prescribe a profile for the decay of Saturn’s surface winds toward the
planet’s center.

  4. 4.

Use a model based on geostrophic (thermal wind) balance to estimate the
density perturbations and resulting shifts $\delta J_{n}$ to the zonal gravity
stemming from the differential rotation.

  5. 5.

Perform a perturbative calculation of the frequency shift each mode attains as
a result of the differential rotation.

  6. 6.

Finally evaluate each model against the observed $J_{2n}^{\rm
obs}=J_{2n}+\delta J_{2n}$ and ring wave pattern speeds $\Omega_{{\rm
p},i}=\sigma_{i}/m$ through a multivariate normal likelihood function $L$. A
Markov Chain Monte Carlo algorithm samples interior models from this
likelihood function, providing posterior distributions of the model parameters
of interest.

Each of these steps is detailed in the following.

###  2.1 Rigidly rotating background

We begin with rigidly rotating models for Saturn’s interior structure. The
initial step of solving for the rotating hydrostatic equilibrium with
prescribed compositional and thermal structure and tabulated equations of
state (EOS) is an iterative process, making the fully two-dimensional methods
of D21 an impractical starting point. Instead we start by solving for the
oblate equilibrium structure using a seventh-order perturbative theory of
figures (ToF7; Zharkov & Trubitsyn 1978; Nettelmann et al. 2021). The material
is modeled as a mixture of H (mass fraction $X$), He (mass fraction $Y$),
silicates, and ices by combining the “MH13-SCvH” EOS for H/He (Militzer &
Hubbard, 2013; Miguel et al., 2016) with the ANEOS (Thompson, 1990) EOS for
serpentine and water ice. Constituent species are combined assuming an
idealized mixture of noninteracting species, such that intensive quantities
like the specific volume $1/\rho$ are sums over that of each species, weighted
by the species’ mass fractions. Note however that nonideal contributions to
the entropy of H/He mixtures are inherently included in the _ab initio_ based
Militzer & Hubbard (2013) EOS.

We take Saturn’s composition structure to include a stably stratified dilute
core, the simplest interior structure permitted by the low degree ring
seismology (MF21, D21). We take the further step of assuming that this
gradient extends to Saturn’s center, leaving its outer radial extent $r_{\rm
out}$ as a free parameter such that the heavy element mass fraction follows

$Z(r)=\left\\{\begin{array}[]{l}Z_{\rm out}+(Z_{\rm in}-Z_{\rm
out})\sin^{2}\left[\frac{\pi}{2}\left(1-\frac{r}{r_{\rm
out}}\right)\right],\\\ \hfill r<r_{\rm out};\\\ Z_{\rm out},\hfill r\geq
r_{\rm out}.\end{array}\right.$ (1)

While MF21 could not rule out the existence of a small ($\lesssim\\!5\ M_{\rm
E}$) chemically homogeneous, adiabatically or sub-adiabatically stratified
inner core, much of the constraining power in the present analysis is in the
high-degree ($\ell\gtrsim 5$) f-modes whose eigenfunctions are localized close
to Saturn’s surface ($r/R_{\rm eq}\gtrsim 0.8$) and so the details of the
inner core are not crucial. Finally, in lieu of direct knowledge of Saturn’s
helium distribution, we presume as in MF21 that the helium abundance is
stratified over the same dilute core region such that an equation analogous to
Eq. 1 holds with $Y^{\prime}\equiv Y/(X+Y)=Y/(1-Z)$ appearing in place of $Z$
and $Y_{\rm in}^{\prime}$ and $Y_{\rm out}^{\prime}$ defined analogous to
$Z_{\rm in}$ and $Z_{\rm out}$. This parameterization allows for a range of
possible helium profiles, from homogeneous to bottom-heavy distributions
consistent with hydrogen-helium phase separation, the quantitative
understanding of which is still greatly uncertain (Helled et al., 2020; Brygoo
et al., 2021).

The condition that the mean of $Y^{\prime}$ over Saturn’s volume satisfies the
protosolar helium abundance $Y^{\prime}_{\rm proto}=0.275\pm 0.001$ (e.g.,
Lodders 2003; Asplund et al. 2009) eliminates $Y_{\rm in}^{\prime}$. Tests
treating $Y^{\prime}_{\rm proto}$ as a free parameter confirm that this
uncertainty, here estimated from the range of available literature values,
bears negligibly on our findings. Fixing the equatorial radius at $P=1\ {\rm
bar}$ to 60,268 km (Lindal et al., 1985) and demanding that the model’s total
mass satisfies Saturn’s total mass in turn eliminates $Z_{\rm in}$. We assume
that the _thermal_ stratification is adiabatic throughout Saturn’s interior,
including the diffuse core, and fix the temperature at $P=1\ {\rm bar}$ to be
$T_{1}=135\ {\rm K}$ (Lindal et al., 1985). Although the diffuse core likely
has a superadiabatic temperature gradient associated with double-diffusive
convection, its effect on fluid stability is overwhelmed by that of the
composition gradient (MF21) and is hence negligible in the present analysis.
We are left with a handful of free parameters $Z_{\rm out},\ Y_{\rm
out}^{\prime},\ r_{\rm out},$ and the angular spin frequency $\Omega_{\rm S}$
of the rigidly rotating deep interior.

Figure 1:  Errors in zonal gravity harmonics for the $n=1$ polytrope at
$q\equiv(\Omega^{2}R_{\rm eq}^{3}/GM)_{\rm Jup}=0.089195487$, compared to the
exact Bessel solution given by Wisdom & Hubbard (2016). Open symbols indicate
ToF4 results and filled symbols indicate ToF7 results. ToF4 only yields even
coefficients up to $J_{8}$. Shaded regions indicate our working resolution of
4,096 level surfaces and the corresponding $J_{2n}$ errors in ToF7 dictate the
variances adopted in the gravity likelihood function $L_{\rm grav}$ (see
Eq.8).

The fidelity of the calculated rigid-body $J_{2n}$ is limited by the accuracy
of our ToF7 implementation, which we benchmark using an $n=1$ polytrope for
which a precise solution for the $J_{2n}$ is known (Wisdom & Hubbard, 2016;
Nettelmann et al., 2021). The resulting errors are shown as a function of
spatial resolution in Figure 1. From the results of the benchmark we conclude
that 4,096 level surfaces offers an acceptable accuracy, and adopt the
corresponding fractional errors $(\Delta J_{2}/J_{2},\Delta J_{4}/J_{4},\Delta
J_{6}/J_{6})=(8\times 10^{-5},7\times 10^{-4},1\times 10^{-3})$ as our
effective uncertainties. These errors, significantly larger than the
measurement uncertainties (Iess et al., 2019; Jacobson, 2022), are similar to
the residuals found at comparable resolution in the independent ToF7
implementations of Nettelmann and Movshovitz (Nettelmann et al., 2021).

Although Figure 1 shows typically one-sided errors for each $J_{2n}$, we find
that in general, these errors can be positive or negative depending on the
particular model in question. Hence, we conservatively treat the $\Delta
J_{2n}/J_{2n}$ observed from this figure as symmetric (two-sided)
uncertainties.

###  2.2 Solving for normal modes

Normal modes are first computed for the rigidly rotating background using the
method described by D21 that includes a fully nonperturbative treatment for
the effects of the Coriolis force and centrifugal distortion arising from
Saturn’s rapid rotation. We refer to that work for the full details of the
calculation. These pseudospectral calculations employ nonspherical coordinates
$(\zeta,\mu=\cos\theta)$ with $2N_{\mu}+1$ grid points in latitude and
$N_{\zeta}$ grid points in the quasi-radial coordinate $\zeta$; see Appendix A
in D21 for an explicit definition of these coordinates. The density $\rho$,
pressure $P$, Brunt-Väisälä (buoyancy) frequency $N$, and adiabatic sound
speed $c_{s}$ are interpolated from the ToF model onto the much lower radial
resolution oblate model used for the mode calculations. Aiming to balance
frequency fidelity against computational expense, we find that a rather modest
radial resolution $N_{\zeta}=40$ is sufficient for pattern speeds of the
relevant modes to be converged to within 0.1 ${\rm deg\ d}^{-1}$ (a fractional
frequency error $\delta\sigma/\sigma\lesssim 0.01\%$), matching a conservative
estimate of the precision on measured ring pattern speeds (French et al.,
2021). Computation time scales less steeply with $N_{\mu}$, which we fix at
61. Finally, spectral convergence is critical: rotation endows any
equatorially symmetric mode333The same truncation parameter $N_{\tilde{\ell}}$
applies to antisymmetric modes, which get contributions from
$\tilde{\ell}=m+1,m+3,m+5,\ldots$. of a given $m$ value with contributions
from angular degrees $\tilde{\ell}=m,m+2,m+4,\ldots$ and the number
$N_{\tilde{\ell}}$ of terms retained in the spectral expansion must be chosen
judiciously. We find that for high-degree f-modes, up to 12 values of
$\tilde{\ell}$ contribute to an extent that affects the frequencies at
$\delta\sigma/\sigma\lesssim 0.01\%$. Hence when searching a given value of
$m$ and modes up to a given $\ell-m$ are desired, retaining
$N_{\tilde{\ell}}=12+\ell-m$ terms is sufficient. The lower degree ($\ell\leq
5$) modes that we fit in this work are more spectrally rich and require
$N_{\tilde{\ell}}=20-24$. In what follows, the modes under discussion always
have an unambiguous _dominant_ angular degree that we label simply as $\ell$.

Figure 2:  _Left:_ Residuals (model minus data) of C ring density and bending
wave pattern speeds for the best fitting rigidly rotating Saturn model.
Different Saturn modes are distinguished by their azimuthal wavenumber $m$
(vertical axis) and latitudinal wavenumber $\ell-m$ (see colors in legend).
Solid circles are fully nonperturbative calculations for mode frequency; open
circles highlight modes calculated on the fly and fit to data during MCMC
sampling (see Section 2.6). Error bars on these modes indicate 5-95% quantiles
of their residuals across the sample of rigid models. The systematic offsets
seen for this optimized rigid model, most obvious for the sectoral
($\ell-m=0$) modes, betray the influence of Saturn’s deep winds. _Right:_
Gravitational moments $J_{2n}$ of the rigid model compared to data. Top panel
plots the negative $J_{2n}$ as open circles, positive coefficients as filled
circles. Bottom panel plots their absolute residuals with respect to the
observed $J_{2n}$ (Iess et al., 2019).

Because our present emphasis is on Saturn’s rotation profile rather than deep
interior structure, and we require a lightweight enough calculation to
efficiently sample a large number of models, we sidestep some of the lowest
$m$ modes where multiplicities of observed ring waves arise from mixtures of
fundamental and gravito-inertial modes (Fuller 2014, MF21, D21). A complete
description of this phenomenon arising in the $m=2$ and $m=3$ wave detections
may require more than one g-mode cavity in the interior. Instead we limit real
time mode calculations to f-modes with $m\geq 5$, where all Saturn-associated
ring waves appear to be excited by well isolated f-modes residing mostly in
the envelope (M19, French et al. 2021). To these we add a single deep-seated
mode at each of $m=2,3,4$. At $m=3$ and $m=4$ these are $\ell=m$ f-modes that
interact somewhat with nearby g-modes. The more nuanced case of $m=3$ involves
three closely spaced detections that likely arise from interactions between
${}_{3}^{3}f$ and nearby g-modes or rosette modes (Fuller 2014, D21); in this
work we simply compute ${}_{3}^{3}f$ and compare to the strongest $m=3$ wave
W82.06. At $m=2$ we focus on the deep-seated mode ${}_{2}^{2}g_{1}$ that MF21
ascribed to the fastest $m=2$ density wave W76.44. As in MF21 this mode
provides a crucial anchor for $N(r)$, mitigating otherwise major degeneracies
between $Z_{\rm out}$, $r_{\rm out}$, and $Y_{\rm out}^{\prime}$ seen when
fitting only the zonal gravity field. ${}_{2}^{2}g_{1}$ also has the major
practical advantage of not mixing strongly with higher degree ($\ell\gg 2$)
g-modes, unlike other $m=2$ and $m=3$ waves where strong interactions between
the $\ell=m$ f-modes and $\ell\geq m$ modes may necessitate $N_{\ell}$ and
$N_{\zeta}$ of order $100$, the expense of which puts these outside our
present scope.

Wave | Mode (${}_{\,\ell}^{m}\cdot_{n}$) | Fit in MCMC | Reference
---|---|---|---
W76.44 | ${}_{2}^{2}g_{1}$ | y | F19
W82.06 | ${}_{3}^{3}f$ | y | H19
W74.93 | ${}_{5}^{4}f$ | - | F19
W80.98 | ${}_{4}^{4}f$ | y | H19
W81.02a | ${}_{5}^{5}f$ | y | H19
W81.43 | ${}_{6}^{6}f$ | y | H19
W81.96 | ${}_{7}^{7}f$ | y | H19
W76.46 | ${}_{9}^{7}f$ | y | F19
W74.67 | ${}_{10}^{\ 7}f$ | - | F19
W82.53 | ${}_{8}^{8}f$ | y | H19
W76.24 | ${}_{11}^{\ 8}f$ | - | F19
W74.51 | ${}_{12}^{\ 8}f$ | y | F21
W74.94 | ${}_{14}^{\ 9}f$ | - | F19
W83.09 | ${}_{9}^{9}f$ | - | H19
W79.04 | ${}_{11}^{\ 9}f$ | - | F21
W76.02 | ${}_{13}^{\ 9}f$ | - | F19
W83.63 | ${}_{10}^{10}f$ | - | H19
W77.34 | ${}_{14}^{10}f$ | - | F21
W81.87 | ${}_{11}^{10}f$ | - | F21
W75.14 | ${}_{16}^{10}f$ | - | F21
W78.51 | ${}_{15}^{11}f$ | y | F21
W81.02b | ${}_{13}^{11}f$ | y | H19
W84.15 | ${}_{11}^{11}f$ | y | F21
W79.84 | ${}_{14}^{11}f$ | - | F21
W74.76 | ${}_{19}^{11}f$ | - | F21
W79.55 | ${}_{16}^{12}f$ | - | F21
W82.61 | ${}_{15}^{13}f$ | y | F21
W80.49 | ${}_{17}^{13}f$ | y | F21
W81.33 | ${}_{18}^{14}f$ | - | F21
Table 1:  Associations between Saturn modes and resonant features observed in
Saturn’s C ring. F19: French et al. (2019); H19: Hedman et al. (2019); F21:
French et al. (2021).

For the model under consideration, many of the mode pattern speeds are
strongly covariant, meaning that little information is gained from the arduous
process of calculating and fitting the full set of observed modes for every
model. Instead useful subsets are chosen to represent a diversity of angular
degrees $\ell$ and latitudinal wavenumbers $\ell-m$, and models of interest
are followed up with complete calculations.

Figure 2 summarizes the best fit we find for a rigidly rotating Saturn fit to
$J_{2}$, $J_{4}$, $J_{6}$, and the pattern speeds of 14 C ring waves. The full
list of wave-mode associations used in this work, as well as the subset fit on
the fly, is given in Table 1. The full model sampling process is detailed in
Section 2.6 below; here it suffices to say that this rigid model is for an
optimized combination of $Z_{\rm out}$, $r_{\rm out}$, $Y_{\rm out}^{\prime}$,
and $2\pi/\Omega_{\rm S}$. Following up on the full set of modes with
corresponding density wave detections, we find pattern speed residuals ranging
from $-2$ to $+3\ {\rm deg\ d}^{-1}$, amounting to an overall spread of
approximately 50 times the fiducial observational uncertainty $0.1\,{\rm deg\
d}^{-1}$. Clearly rigid rotation is not an adequate description, even with a
considerable amount of freedom in the composition structure and planetary
spin.

The optimized rigid model’s 632.9 min spin period agrees with the initial ring
seismology estimate of $633.6^{+1.9}_{-1.3}\ {\rm min}$ reported by M19. But
in contrast to that work, the improved accuracy of the nonperturbative method
reveals significant systematic frequency residuals, most noticeable for the
sectoral ($\ell-m=0$) modes, that betray the influence of Saturn’s deep winds.
The sectoral mode frequencies are generally underestimated relative to the
others, suggesting a prograde differential rotation concentrated near the
equator where the sectoral modes are predominantly confined. Furthermore this
systematic underestimate gets worse with increasing angular degree $\ell$,
notable because f-modes with higher $\ell$ are confined closer to Saturn’s
surface. Together these outcomes point toward a prograde differential rotation
confined mostly to the equator and strongest near Saturn’s surface, consistent
with the dominant feature of Saturn’s zonal winds. From the perspective of the
f-mode frequencies there is an inherent tradeoff between the depth or
amplitude of this prograde feature and the unknown spin frequency of the deep
interior, and hence in a complete model addressing Saturn’s deep winds, we
expect slower deep spin than in the rigid model shown in Figure 2 or in M19.

The frequency residuals in Figure 2 differ from those previously obtained in
M19 (their Figure 5) and French et al. (2021) (their Figure 10) for several
reasons. First, the new model is stably stratified at depth, fundamentally
altering the spectrum by introducing g-modes and strongly affecting the
frequencies of all $m\lesssim 5$ modes. Second, the new model’s improved
treatment for rotation in the mode calculations affects all modeled
frequencies, most strongly those with larger values of $m$, leading to clear
differences compared to the perturbative calculations shown in M19 and French
et al. (2021). The present model’s improved $J_{2n}$ accuracy and inclusion of
$J_{6}$ in the fit also contribute to the slightly faster bulk rotation
favored here, driving the more obvious systematic trend in frequency residuals
with $m$ than was seen in M19 or French et al. (2021). In any case, we
emphasize that this optimized rigid model is a poor fit to seismology and
gravity both; understanding Saturn’s differential rotation is essential.

Figure 3:  Example profiles of baroclinic differential rotation satisfying
Saturn’s observed cloud-level winds (García-Melendo et al., 2011), shown for
fast, medium, and slow deep spin $\Omega_{\rm S}$ (left, center, and right
columns) and with wind decay depth set to $0pt=0.12$. Depending on the deep
spin, the winds can range from predominantly retrograde to predominantly
prograde.

###  2.3 Baroclinic wind model

The winds observed at Saturn’s surface are taken to extend inward on surfaces
parallel to the rotation axis, before multiplication by a radial decay
function centered at $r=r_{\rm surf}-d$ smoothly transitions to a rigidly
rotating interior. In particular the spin frequency is given by
$\boldsymbol{\Omega}(r,\mu)=\hat{\boldsymbol{z}}\left[\Omega_{\rm
S}+\delta\Omega(r,\mu)\right]$ where

$\begin{split}\delta\Omega(r,\mu)&=\delta\Omega_{\rm surf}(\tilde{\mu})\\\
&\times\frac{1}{2}\left(1+\tanh\left[\frac{r-\left[r_{\rm
surf}(\mu)-d\right]}{w}\right]\right).\end{split}$ (2)

Here $\mu=\cos\theta$, $r_{\rm surf}=r_{\rm surf}(\mu)$ is the spherical
radius coordinate of the oblate $P=1$ bar surface, and $\tilde{\mu}$
corresponds to the latitude on this surface sharing the same cylindrical
radius as an interior point $(r,\mu)$. The depth parameter $d$ controls the
centroid of the decay function, and $w$ is the decay width. Tests reveal that
wind-induced perturbations to mode frequencies and gravity harmonics are only
weakly sensitive to the decay width $w$; in what follows we fix $w=0.02$.

The assumed surface winds $\delta\Omega_{\rm surf}(\mu)$ follow from a
Legendre polynomial expansion of the average wind speeds observed in the ISS
narrowband continuum filters CB2/CB3 (García-Melendo et al., 2011), normalized
to a given assumed deep rotation rate $\Omega_{\rm S}$. We also consider a
slightly modified surface wind profile (“Galanti reconstructed”; E. Galanti,
private communication) whose main departure from the García-Melendo winds is
its enhanced retrograde rotation in the region $0.4<|\mu|<0.6$. This is
anticipating that a rotation law of the form 2 and adhering strictly to the
García-Melendo cloud-level velocity profile may have difficulty reproducing
Saturn’s observed odd and higher degree even ($J_{6+}$) gravity harmonics
(Galanti et al., 2019; Galanti & Kaspi, 2021).

In both cases we find that a 60th order polynomial and our 123 zone
latitudinal grid resolution (Section 2.2) are sufficient to resolve the
$J_{n}$ and $\Omega_{{\rm p},i}$ to within our fit tolerances. We retain both
equatorially symmetric and antisymmetric terms in the expansion. Although we
will see that the antisymmetric components have a negligible effect on the
frequencies of isolated f-modes, we do wish to consider the asymmetric gravity
moments ($J_{3},\ J_{5},\ \ldots$) induced by the winds. In all cases wind
speeds are forced to vanish near the poles using an exponential decay at
$|\mu|>0.975$, and antisymmetric terms are similarly attenutated close to the
midplane to avoid discontinuities; for further details see D21. Sample
profiles of the differential rotation model are depicted in Figure 3.

###  2.4 Wind-induced perturbations to zonal gravity

Saturn’s winds modify the planet’s equilibrium structure and gravity. We
calculate shifts in the gravitational coefficients $\delta J_{n}$ by solving a
“thermal gravitational wind equation” (TGWE; Zhang et al. 2015; Cao &
Stevenson 2017a) that enforces geostrophic balance. Appendix A describes the
numerical approach we take to this problem. In short, we solve the coupled
system of Euler and Poisson equations governing small-amplitude, time-
independent density and gravity perturbations $\rho^{\prime}({\bf r})$,
$\Phi^{\prime}({\bf r})$ that result from imposing a wind with a given decay
profile.

These calculations build on the common practice of solving the “thermal wind
equation” (TWE) in a spherical geometry (Galanti et al., 2019; Galanti &
Kaspi, 2021) with the inclusion of self-gravity and the already oblate
background of the rigidly rotating model. Like Cao & Stevenson (2017a), we do
not find that self-gravity hugely affects the $\delta J_{n}$ for $n\geq 6$,
but including it does allow for the simple imposition of the constraint that
the total mass remain unchanged (i.e., that $\delta J_{0}=0$). On the other
hand, surface oblateness induced by Saturn’s rapid bulk rotation can have a
larger effect, depending on the differential rotation profile (see Figure 14).

Note that by construction both the TWE and TGWE neglect irrotational
contributions to the Euler equation (Cao & Stevenson, 2017a). While our
comparisons with fully nonlinear calculations for differentially rotating
polytropes (see Figure 14) suggest that the TGWE provides accurate $\delta
J_{n}$ for $n\geq 6$, fully nonlinear calculations of differentially rotating
equilibria are desirable.

Figure 4:  Wind-induced zonal gravity perturbations $\delta J_{2n}$ calculated
assuming thermal-gravitational wind balance for baroclinic rotation laws of
the form given in Equation 2, as a function of assumed deep spin and wind
decay depth. These points are for our baseline sample that adopts a modified
cloud-level wind (see Section 2.4). Color bars are labeled in parts per
million.

Figure 4 summarizes the $\delta J_{2n}$ predicted by the oblate TGWE over the
baseline sample of differentially rotating models that we obtain below. These
models assume the Galanti reconstructed cloud-level winds. Winds enhance
$J_{6}$ by typically $\delta J_{6}\approx 3$ ppm.

###  2.5 Wind-induced perturbations to mode frequencies

####  2.5.1 Method

Figure 5:  Wind-induced perturbations to f-mode frequencies, calculated from
Equation 3, expressed here in terms of inertial frame frequency
$\sigma_{i}=\omega_{i}+\delta\omega_{i}$. After Figure 4 in Dewberry et al.
(2022).

Saturn’s $\lesssim 300$ m s-1 wind speeds at the cloud level translate to a
modest differential rotation amplitude $|\delta\Omega/\Omega_{\rm S}|\lesssim
3\%$. The order of magnitude influence on frequencies of oscillation modes is
even smaller at $|\delta\Omega/\Omega_{\rm dyn}|\lesssim 1\%$, motivating the
use of a perturbation theory to calculate the leading order frequency shifts
from differential rotation. We accordingly calculate these shifts accounting
for the self-coupling induced in each mode by the winds, a combination of
Coriolis and centrifugal effects, as (D21 Equation 24)

$\delta\omega_{i}\simeq\frac{-(2\omega_{i}\mathcal{W}_{ii}^{D}+\mathcal{U}_{ii}^{D})}{2(\omega_{i}\mathcal{T}_{ii}+\mathcal{W}_{ii}^{R}+\mathcal{W}_{ii}^{D})}$
(3)

where $\omega_{i}$ is the frequency of the $i$th mode under rigid rotation,
$\mathcal{W}_{ii}^{D}$ and $\mathcal{W}_{ii}^{R}$ are inner products involving
the Coriolis operators associated with differential and rigid rotation
respectively, $\mathcal{T}_{ii}$ corresponds to kinetic energy and is equal to
the identity matrix in our chosen normalization, and $\mathcal{U}_{ii}^{D}$
corresponds to a potential energy operator involving cross terms between the
rigid and differential parts of the rotation. D21 gives full integral
expressions for these operators. Figure 5 gives the frequency shifts
calculated from Equation 3 for a subset of the $\ell\geq 6$ modes in a
representative model.

This approach is accurate insomuch as (i) the frequency shifts are not large
enough to cause accidental degeneracies between different modes with the same
$m$, which D21 showed would lead to wind-induced mode mixing, and (ii) higher-
order consequences of the differential rotation, such as feedback from wind-
induced centrifugal flattening on the mode frequency, are negligible.
Condition (i) is an excellent approximation for the high-degree f-modes, which
as we have mentioned are extremely well isolated in frequency within a single
$m$ family444The co-rotating f-mode frequencies scale roughly as $\omega_{\rm
f}/\Omega_{\rm dyn}\sim\ell^{1/2}$ so that within a given $m$, each f-mode is
separated from the next by at least $\Delta\omega_{\rm f}/\Omega_{\rm
dyn}\sim(\ell+1)^{1/2}-\ell^{1/2}\approx 20\%$ at $\ell=6$ and $10\%$ at
$\ell=20$, safely an order of magnitude larger than the wind-induced shifts..
Condition (ii) is more dubious, and with this question in mind Dewberry et al.
(2022) apply a complete nonperturbative treatment for Saturn-like differential
rotation in an $n=1$ polytrope to conclude that the simpler perturbative
approach adopted here and in D21 overestimates the wind-induced frequency
shifts by up to 0.1-0.2%. This small but significant offset is most noticeable
for the highest degree sectoral ($\ell-m=0$) modes and can be attributed
mostly to the perturbative method’s neglect of the wind-induced modification
to the equilibrium shape.

Because fully nonperturbative calculations along the lines of Dewberry et al.
(2022) are not readily possible for the more realistic (non-polytropic)
equation of state based models considered here, we caution that our
perturbative calculation (3) may slightly overestimate the true mode
frequencies. We expect however that the frequency errors introduced by the
perturbation theory will be smaller than in Dewberry et al. 2022, mainly
because of differences in the decay profile. Below we show that a suitable fit
to the gravity field (particularly $J_{6}$) requires a modified cloud-level
wind similar to that of Galanti et al. (2019) (see Section 2.4). This leads to
a more strongly sub-corotating region at latitudes $0.4<|\mu|<0.6$ (Figure 7),
mitigating somewhat the effect of the predominant prograde equatorial jet.
Hence we expect smaller frequency perturbations from the winds in general,
especially for the tesseral ($\ell-m>0$) modes that sample higher latitudes.
We also expect weaker feedback on frequencies from the wind-induced
perturbation to the planetary shape, the predominant effect responsible for
the perturbative/nonperturbative frequency offsets observed by Dewberry et al.
(2022). Futhermore half the frequencies we fit on the fly correspond to
$\ell-m>0$ modes, which Dewberry et al. (2022) found are dramatically less
affected by the perturbative approximation, their errors insignificant at
$\delta\sigma/\sigma\lesssim 0.05\%$. We finally note that our formulation
takes Saturn’s equatorial radius as the reference length scale for non-
dimensionalization, and hence excluding the shape change amounts to ignoring
wind-induced displacement only at higher latitudes.

Figure 6:  Meridional slices of the first-order Coriolis kernels (left
hemispheres; Equation 7) and the resulting integrands of the first-order
frequency shift (right hemispheres; Equation 6) for the four f-modes with
$m=11$ and $\ell-m=0,\ 2,\ 4,\ 6$ in our favored Saturn model. The dashed grey
boundary marks the centroid of the wind decay function. All four modes are
shifted to higher frequency by the prograde equatorial jet, although modes
with larger $\ell-m$ acquire smaller shifts because they sample higher
latitudes where they suffer self-cancellation. Color bars are normalized
arbitrarily; values of the frequency shift are given in the text. Note that
the perturbative treatment of winds in our main calculations captures higher-
order rotation effects than the heuristic first-order perturbations shown
here.  Figure 7:  Examples of a few modes’ latitudinal sensitivity to
rotation: as in Figure 6, but comparing the Legendre expansion of the cloud-
level wind angular frequency with the radial peak amplitude reached by
$K(r,\mu)$ (Equations 6-7) as a function of planetocentric latitude for each
of the three modes (black curves with labels, normalized to unity). Small
orange points are the mean measured wind profile from Casini imaging (García-
Melendo et al., 2011) and the solid blue curve is a gravity-optimized profile
from Galanti & Kaspi (2021); both are scaled to 635 min baseline rotation
here. Modes with $\ell-m=0$ essentially only sample the prograde equatorial
jet; $\ell-m=2$ and $4$ are more sensitive to the strongly retrograde regions
from $30^{\circ}$ to $40^{\circ}$ latitude.

####  2.5.2 Which modes probe which regions?

A useful demonstration of the ring seismology’s spatial sensitivity to
Saturn’s differential rotation can be gleaned from the perturbative shifts
(Equation 3). Inspection of the inner products appearing in that equation for
high-degree ($\ell\gtrsim 5$) f-modes in a Saturn model shows that to a fair
approximation, we can neglect $\mathcal{U}_{ii}^{D}$ relative to
$2\omega_{i}\mathcal{W}_{ii}^{D}$ in the numerator, and neglect
$\mathcal{W}_{ii}^{R}+\mathcal{W}_{ii}^{D}$ relative to
$\omega_{i}\mathcal{T}_{ii}$ in the denominator, leaving

$\displaystyle\delta\omega_{i}$
$\displaystyle\simeq-\mathcal{W}_{ii}^{D}/\mathcal{T}_{ii}$ (4)
$\displaystyle=\frac{\int\rho\,\delta\Omega\,\boldsymbol{\xi}_{i}^{*}\cdot\left(m\boldsymbol{\xi}_{i}-{\rm
i}\hat{\boldsymbol{z}}\times\boldsymbol{\xi}_{i}\right)\,dV}{\int\rho\,\boldsymbol{\xi}_{i}^{*}\cdot\boldsymbol{\xi}_{i}\,dV}$
(5)

with $\boldsymbol{\xi}_{i}$ denoting the Lagrangian displacement eigenfunction
for the $i$th mode. Then the frequency shift can be succinctly written as
(dropping the subscripts $i$)

$\delta\omega\simeq\int\delta\Omega(r,\mu)K(r,\mu)\,dV$ (6)

by defining the first-order Coriolis kernel

$K=\frac{{\rho\left[m|\boldsymbol{\xi}|^{2}-{\rm
i}\boldsymbol{\xi}^{*}\cdot(\boldsymbol{\hat{z}}\times\boldsymbol{\xi})\right]}}{\int\rho|\boldsymbol{\xi}|^{2}\,dV},$
(7)

a field independent of the assumed rotation model (and equivalent to, e.g.,
Equation 3.332 of Aerts et al. 2010). Here we use this simplified rotation
kernel for heuristic purposes; our main results will use the full shifts given
by Equation 3.

Figure 6 displays meridional slices of $K(r,\mu)$ and the $\delta\omega$
integrand (Equation 6) for the four $m=11$ modes ${}_{11}^{11}f$,
${}_{13}^{11}f$, ${}_{15}^{11}f$, and ${}_{19}^{11}f$ that have been detected
as OLRs in the rings. The high-$\ell$ f-modes are confined close to the
surface, and lower latitudinal wavenumbers $\ell-m$ are confined closer to the
equator. Consequently the sectoral modes predominantly sample Saturn’s
prograde equatorial jet, with tesseral ($\ell-m>0$) modes sampling rotation
both at the prograde equator and alternating prograde/retrograde rotation at
higher latitudes (D21). Furthermore the fact that $K(r,\mu)$ is strictly
positive means that prograde regions always increase mode frequency, while
retrograde regions decrease mode frequency. The frequency shifts induced by
winds are thus typically positive for high-degree f-modes, albeit with smaller
magnitudes toward larger $\ell-m$ as a result of offsetting contributions from
prograde and retrograde winds. Note also that the equatorial symmetry
$K(r,\mu)=K(r,-\mu)$ implies that the small antisymmetric components of the
winds have _zero_ first order effect on the frequency of any mode thanks to
the vanishing integral 6.

In our full calculation (Equation 3), the term $\mathcal{U}_{ii}^{D}$ can be
nonzero and is largest for high $\ell$, high $m$ modes, where it nonetheless
does not exceed magnitudes $\sim 10^{-6}-10^{-5}~{}\Omega_{\rm dyn}$ and can
safely be ignored. A complication would arise should the asymmetric wind cause
an avoided crossing between a mode of interest and another mode with opposite
equatorial parity. As shown by D21, antisymmetric components of the wind can
in fact mix low-degree f-modes with nearby low-degree g-modes with opposite
parity, a phenomenon that does not appear to play a role for any of the modes
that we focus on here.

Figure 7 presents a one-dimensional view of these ideas for the same three
modes, emphasizing the latitudinal overlap between these modes’ rotation
kernels and the cloud-level wind speed.

###  2.6 Parameter estimation

We proceed much as in MF21, calculating the total likelihood as the product
$L=L_{\rm seis}\times L_{\rm grav}$ of two multivariate Gaussians such that

$\ln L=-\sum_{i}\frac{(\Omega_{\rm p}^{\rm calc}-\Omega_{\rm p}^{\rm
obs})^{2}}{2\,\Delta\Omega_{\rm p}^{2}}-\sum_{n=1}^{3}\frac{(J_{2n}^{\rm
calc}-J_{2n}^{\rm obs})^{2}}{2\,\Delta J_{2n}^{2}}.$ (8)

$\Delta\Omega_{\rm p}$ is taken to be uniform at $0.3\ {\rm deg\ d}^{-1}$,
approximately $0.02\%$ of the pattern speeds themselves, slightly larger than
the fiducial observational uncertainty $\approx 0.1\ {\rm deg\ d}^{-1}$
(Section 2.2). The slightly larger uncertainty aids convergence of the
sampling algorithm by retaining a larger fraction of proposed models. We also
compare against an alternative sample with an inflated uncertainty
$\Delta\Omega_{\rm p}=3\ {\rm deg\ d}^{-1}$ ($0.2\%$), and a sample that
neglects the ring seismology altogether, to emphasize the value of precision
seismology for constraining the rotation. The adopted gravity uncertainties
$\Delta J_{2n}$ are given in Section 2.1.

We consider only uniform distributions for prior probability, so the posterior
probability is simply proportional to $L$ inside the allowed prior volume and
zero outside of it. Aiming to consider the largest plausible prior volume, we
adopt the following:

  1. 1.

$Z_{\rm out}\in[0.015,0.3]$, i.e., the outer envelope metallicity is at least
solar. Since $Z$ is assumed to increase with depth, $Z_{\rm out}$ is
inevitably bounded from above by the mean metallicity consistent with Saturn’s
mean density, approximately $\langle Z\rangle\lesssim 0.2$ (e.g., Guillot,
1999) or very conservatively 0.3.

  2. 2.

$r_{\rm out}\in[0,1]$; in practice $J_{2n}$ and the $m=2$ ring seismology
constraint are expected to strongly favor $r_{\rm out}\sim 0.6$ (MF21).

  3. 3.

$Y_{\rm in}^{\prime}\in[0.275,1]$; the helium distribution lies between the
limiting cases of homogeneous at the protosolar He/H fraction ($Y_{\rm
in}^{\prime}=Y_{\rm out}^{\prime}=0.275$) and maximally stratified ($Y_{\rm
in}^{\prime}=1$; $Y_{\rm out}^{\prime}<0.275$).

  4. 4.

$2\pi/\Omega_{\rm S}\,({\rm min})\in[630,640]$, spanning the range of
estimates for Saturn’s deep spin from various sources (Desch & Kaiser, 1981;
Read et al., 2009; Helled et al., 2015; Militzer et al., 2019). This rules out
the longer periods measured in Saturn’s kilometric radiation (e.g., Giampieri
et al., 2006) as being directly tied to Saturn’s deep rotation; our results
will validate this choice.

  5. 5.

$d/R_{\rm eq}\in[0.05,0.2]$. The lower limit is nonzero to ensure that our
radial decay function does not artificially modify the cloud-level winds.
Models near this lower limit have f-mode spectra very similar to the the
rigidly rotating case (Figure 2), suffering a clear systematic mismatch in the
sectoral $(\ell=m)$ modes and precluding any viable solution with $d/R_{\rm
eq}<0.05$. The maximum decay centroid depth is chosen to well encompass the
values $\approx 0.12\,R_{\rm eq}\approx 7,000\,{\rm km}$ estimated from
analysis of zonal gravity (Iess et al., 2019; Galanti et al., 2019) and
magnetic field strength (Galanti & Kaspi, 2021). This depth is approximately
where the electrical conductivity, increasing steeply with depth, becomes
sufficiently large (Liu et al., 2008) that Lorentz forces eliminate any
differential flows with velocities comparable the cloud-level winds (Cao &
Stevenson, 2017b; Guillot et al., 2018).

We use `emcee` (Foreman-Mackey et al., 2013) version 2.2.1 to sample these
five parameters against the likelihood of Equation 8 given the first three
even zonal gravity harmonics $J_{2}-J_{6}$ from Iess et al. (2019), normalized
to $R_{\rm eq}=60,268$ km, and the pattern speeds of the same 14 patterns that
are fit by the rigid model in Figure 2. These patterns and the additional
mode-wave associations used to compare models to observations post-hoc are
listed in Table 1. Although the odd gravity harmonics ($J_{3}$, $J_{5}$,
$\ldots$) and higher order even harmonics ($J_{8}$, $J_{10}$, $\ldots$)
dominated by the winds are readily available for all models, these are ignored
in the fitting procedure: our wind model only includes a single parameter
(decay depth) and hence lacks the flexibility to modulate these harmonics
independently from one another (e.g., to modify $J_{8}$ without sacrificing
the fit to $J_{2}-J_{6}$). Instead, models are directly constrained by just
the first three even harmonics, the remainder compared to after the fact.

Additional samples introduce a parameter $f_{\rm ice}\in[0,1]$ controlling the
mass ratio of ices to silicates (see Section 2.1), but we find that silicate-
rich models $f_{\rm ice}\ll 1$ are heavily disfavored, significantly
overestimating the pattern speed of the ${}_{3}^{3}f$ mode by several ${\rm
deg\ d}^{-1}$ and underestimating $|J_{4}|$ by typically 2 ppm. Hence in what
follows we focus on the successful water ice models with $f_{\rm ice}=1$. When
discussing distributions (e.g., of model parameters), all numerical ranges are
given in terms of their 5% and 95% quantiles.

##  3 Results

###  3.1 The cloud level wind profile

Joint analysis of gravity and seismology strongly favor the Galanti modified
cloud-level winds (Section 2.3) over the mean profile from García-Melendo et
al. (2011). First, ignoring ring seismology entirely and fitting just
$J_{2}-J_{6}$, a sample strictly adopting the García-Melendo profile tends to
systematically overestimate $J_{6}$. This overestimate drives the most
successful models to minimize the wind contribution $\delta J_{6}$ by adopting
the shallowest wind depths allowed by our prior $0pt\geq 0.05$. This can lead
to plausible fits to $J_{6}$, but poor fits to the higher degree harmonics
$J_{8}$ and $J_{10}$, and post hoc comparisons to the seismology find poor
fits that are akin to the rigid model in Figure 2 at best.

Figure 8:  Gravity harmonics $J_{n}$ (top) and their residuals (bottom) for
the best model that strictly satisfies the García-Melendo et al. (2011) cloud
level wind profile, and is fit to $J_{2}-J_{6}$ and the ring seismology.
Dashed and solid vertical lines represent uncertainty in the ToF calculation
and the Cassini measurement uncertainty respectively (see legend). The
$J_{2n}$ are overplotted with plus signs. This model requires shallow
($d\approx 3,400\ {\rm km}$) winds to even marginally fit $J_{6}$. The poor
fit to the odd harmonics, even harmonics $J_{6+}$, and seismic constraints
(not pictured) all strongly favor a model with somewhat different surface
winds (see text).

Repeating this exercise but folding in the seismic constraints yield models
that compromise the gravity fit further. The gravity spectrum of the best such
model is shown in Figure 8, where the fit to the higher-degree even harmonics
$J_{6}$, $J_{8}$ and $J_{10}$ is especially poor. Hence our results support
what Galanti et al. (2019) and Galanti & Kaspi (2021) found by applying a
similar functional form for the wind decay: the gravity favor a somewhat
different latitudinal profile than the García-Melendo observations. The
seismology strengthens this finding.

The constancy of Saturn’s cloud level winds over multiple observations and
long time baselines (García-Melendo et al., 2011) suggests that this
discrepancy arises not from some observational bias but from the simplified
decay profiles, including our Equation 2, that have so far been used to
interpret the gravity. The slowly rotating features in Saturnian gravity
uncovered by Hedman et al. (2022) also point toward stronger retrograde
rotation at depth than the surface winds alone would suggest.

In what follows we limit our attention to the much more successful models that
assume the slightly modified “Galanti reconstructed” cloud level wind profile
depicted in Figure 7. As we have mentioned, the principal difference from
observations is in the more pronounced retrograde rotation near
$|\mu|=0.4-0.6$ ($30\textrm{--}45^{\circ}$ planetographic latitude), which
predominantly affects the frequencies of modes with $\ell-m\geq 2$ (Figures
6-7). The modified profile also has a more subtly exaggerated prograde
rotation in the next zone at $|\mu|=0.65-0.75$ ($45\textrm{--}55^{\circ}$),
which the same figures show is probed most strongly by modes with $\ell-m\geq
4$. The wind structure at higher latitudes does not greatly affect the
spectrum thanks to the limited number of detections with large latitudinal
wavenumbers, a point that we return to in this paper’s Discussion.

Figure 9:  Wind model (left), pattern speed residuals for waves observed at
Lindblad or vertical resonances in the C ring (center), and gravity harmonics
(top right) and residuals (bottom right) for our best fitting Saturn model.
Different modes are distinguished by their azimuthal wavenumber $m$ (vertical
axis) and latitudinal wavenumber $\ell-m$ (see colors and markers in legend).
The circled modes are the high-degree modes calculated on the fly and fit to
data during MCMC sampling, their error bars showing the 5-95% quantiles of
residuals obtained by all models across the sample of tens of thousands.

###  3.2 The seismic imprint of Saturn’s winds

Deep zonal winds bring the f-mode frequencies into substantially better accord
with observations. Figure 9 summarizes the best model’s zonal wind structure
and displays pattern speed residuals (calculated minus observed) for Saturn-
associated spiral density and bending waves in the C ring. The winds
ameliorate the worst systematic issues identified in the rigid model, mainly
the strong tendency to underestimate the frequencies of sectoral ($\ell=m$)
modes relative to higher latitudinal wavenumbers $\ell-m>0$, and the tendency
for this bias to get worse toward higher angular degrees (the steep slope of
blue points in Figure 2). For the high-degree ($\ell=m\geq 6$) modes, the
differentially rotating model’s pattern speed residuals span $-0.2$ to $+1.2\
{\rm deg\ d}^{-1}$, a marked improvement from the rigid model’s span of $-2$
to $+3\ {\rm deg\ d}^{-1}$ in Figure 2. This model achieves this with a wind
decay depth $0pt=0.132$ and deep spin period $2\pi/\Omega_{\rm S}=634.8$ min,
the deep spin some 2 min slower than the rigid model as necessary to
compensate for the tendency of the predominantly prograde deep winds to shift
all modes toward higher frequency.

The 6 antisymmetric modes with $\ell-m=1,3,5$ detected via bending waves show
a similar quality of fit, their residuals ranging from 0 to $+1.1\ {\rm deg\
d}^{-1}$. This improvement over the best rigid model’s $+0.2$ to $+1.9$ span
in the $\ell-m=1,3,5$ residuals (Figure 2) is not a result the N-S asymmetric
component of Saturn’s winds, which have a negligible effect on the frequency
of isolated f-modes (see Section 2.5). Rather, the improvement is a direct
result of the symmetric component of the winds, as well as an indirect result
of the same symmetric winds’ impact on the $J_{2n}$ and $\ell=0,2,4$ modes
that these models optimize.

The situation is different for the low-degree ($\ell=m<6$) modes, whose
residuals are very similar to the rigidly rotating case (Figure 2), reflecting
these deep-seated modes’ relatively low sensitivity to the zonal winds and
greater sensitivity to the structure of the deep interior. Their outstanding
$1-6\ {\rm deg\ d}^{-1}$ residuals are much larger than the fiducial pattern
speed uncertainty $0.1\ {\rm deg\ d}^{-1}$, likely a result of the composition
structure posited here not being perfectly accurate, as discussed in MF21; we
revisit this in Section 4. On this topic we also note that the residuals are
overwhelmingly positive, a fact that follows in part from the simultaneous
fitting of the $J_{2n}$ and in part from the outlying ${}_{4}^{4}f$ mode.
Repeating this sample but neglecting the less well understood ${}_{3}^{3}f$,
${}_{4}^{4}f$, and ${}_{5}^{5}f$ modes yields consistent deep spin and wind
decay depths, and does little to alleviate the systematic trends seen in the
residuals at higher $m$ (see Appendix B).

Figure 10:  Decay depth of the zonal winds as a function of Saturn’s deep spin
period. The light blue distribution fits the ring seismology with an
exaggerated pattern speed uncertainty $\Delta\Omega_{p,i}=3$ deg/d while the
dark blue distribution uses $\Delta\Omega_{p,i}=0.3$ deg/d, slightly larger
than the observational uncertainty (see Section 2.2). Plus signs designate the
most likely single model from each sample. Histograms are marginalized
posterior probability distributions. The faint green distribution in the
background ignores the seismology; histograms are omitted for this case. The
error bar laid over the histograms (top) gives the earlier M19 estimate for
Saturn’s spin period based on rigid models and fewer data.

Figure 10 summarizes the posterior distribution of models in deep spin – wind
decay depth space. Here we compare three samples that all fit
$J_{2}\textrm{--}J_{6}$ but incorporate progressively more information: the
first neglects ring seismology, the second fits ring seismology but with
inflated pattern speed uncertainties $\Delta\Omega_{\rm p}=3\ {\rm deg\
d}^{-1}$, and the final favored sample fits the ring seismology with
$\Delta\Omega_{\rm p}=0.3\ {\rm deg\ d}^{-1}$, a conservative value for the
observational uncertainty (see Section 2.2). All three distributions occupy a
similar range in $0pt$, reflecting the fact that the decay depth is most
strongly constrained by $J_{6}$. However, the $0pt$ distribution is both
systematically lower and modestly better constrained in the favored (“Precise
seismology”) sample, demonstrating that the ring pattern speeds also constrain
the depth of the differential rotation. For the favored sample we find a
median depth $0pt=0.132$ (consistent with the maximum likelihood estimate
appearing in Figure 9) and width $0pt=0.125\textrm{--}0.138$.

The other main result in Figure 10 is that compared to models indifferent to
or weakly informed by the ring seismology, the favored sample has a much
better constrained deep spin period. From this sample we find a median deep
spin $2\pi/\Omega_{\rm S}=634.7$ min (close to the maximum-likelihood $634.8$
min period of Figure 9) and width $633.8\textrm{--}635.5$ min. In the case of
the more loosely constrained sample, the distribution peaks at similar deep
spins but is unhelpfully broadly distributed, emphasizing the critical
importance of the frequency precision afforded by Cassini ring seismology. The
final sample that ignores ring seismology is nearly uniformly distributed in
$2\pi/\Omega_{\rm S}$, although the paucity of constraints in this case (3
$J_{2n}$ compared to 5 free parameters) lead to a sample that is not
statistically well converged even after $\sim 10,000$ MCMC iterations. For
this reason we refrain from plotting histograms for this case, but it is
noteworthy that high likelihood models are found throughout the allowed range
of $2\pi/\Omega_{\rm S}$. Neglecting seismology, the single most likely model
of the pure-gravity sample attains a slow deep spin $2\pi/\Omega_{\rm
S}=639.3$ min. This underscores the value of ring seismology for eliminating
much of this parameter space for Saturn’s deep rotation, including periods
$\approx 10^{\rm h}39^{\rm m}$ and longer that appear in Saturn’s kilometric
radiation (e.g., Desch & Kaiser, 1981; Giampieri et al., 2006; Ye et al.,
2018).

Figure 11:  Even gravity harmonics through $J_{12}$ for the three samples
indicated in the legend. Open circles are rigid body components and filled
circles are total harmonics (rigid plus wind contributions). “Strict winds”
refers to the García-Melendo et al. (2011) cloud level velocity profile; the
other two samples assume the Galanti reconstructed profile (see Figure 7 and
Section 2.3). Diamonds with plus signs show the best model from each sample.
Black crosses show the Iess et al. (2019) $3\sigma$ measurement uncertainty
and the red crosses show the ToF7 uncertainty (Section 2.1), slightly offset
here for readability. For ease of comparison, this figure normalizes harmonics
to the conventional reference radius 60,330 km.

Figure 11 shows the full distribution of the $J_{2n}$ in three of our samples,
separating the contributions from the rigid background and the differential
rotation. Here it is clear that even across the entire parameter space
available to the models, models with strict cloud-level winds struggle to fit
$J_{6}$ (by going to very shallow winds; see Section 3.1) and failing more
acutely to fit $J_{8+}$.

Figure 12:  Resonance locations predicted for the best model with winds (small
filled circles) and for the same model but ignoring winds (small open
circles), compared to all observed C ring outer Lindblad resonances (open
black circles) and outer vertical resonances (open black diamonds) excited by
Saturn modes (French et al. 2021 and references therein). Solid lines connect
modes with even $\ell-m$ capable of producing density waves at Lindblad
resonances; dashed lines connect modes with odd $\ell-m$ capable of producing
bending waves at vertical resonances. The shaded region is the D ring where
resonances are not readily observable. Rigid and differential rotation are
indistinguishable at this scale for the $\ell=m$ g1 modes and $m\lesssim 4$
modes.

A summary of the best model’s fit to observations, in terms of absolute
resonance locations in Saturn’s rings, is displayed in Figure 12. Here all
model f-modes with $m=5-14$ and $\ell-m\leq 8$ are shown; for $m=2-4$ the
distinction between f- and g-modes becomes less clear (see MF21, D21) and we
plot only modes with strong surface gravitational potential perturbations,
with lines connecting the modes with the strongest perturbation in each $m$.
Among these $m=2-4$ modes only a subset of detections are used to constrain
the model, and the remainder of these low-$m$ modes do not represent an
exhaustive search, which requires a great deal of radial and spectral
resolution (see Section 2.2). Ongoing work will revisit these in detail.
Toward higher $m$, this broad-scale view shows the overall success of the
model; in particular, the wind model brings the sectoral ($\ell=m$) modes into
excellent agreement with observations.

Figure 13:  Posterior probability distributions of model parameters,
normalized to equal area, compared among the samples appearing in Figure 11.
The grey distribution shows the rigidly rotating sample corresponding to
Figure 2; its narrow $2\pi/\Omega_{\rm S}$ distribution is cut off somewhat
for readability of the other distributions. The favored (dark blue) sample is
the one represented in Figures 9$-$12.

###  3.3 Composition

Figure 13 reports marginalized posterior distributions for the same three
samples as in Figure 11, plus the rigidly rotating sample. In all cases, zonal
gravity constrains these models to 1-3 times protosolar metallicity (here
$Z=0.014$; e.g., Asplund et al. 2009) in the envelope, at odds with the
$9\times$ solar carbon enrichment from spectroscopy (Fletcher et al., 2009).
Higher envelope metallicities lead to unrealistically large magnitudes of the
$J_{2n}$. Part of this discrepancy may be attributable to uncertainties in the
hydrogen EOS in the vicinity of metallization, where fractional density errors
of order a few percent are plausible and can lead to substantial increases in
the envelope metallicity (Militzer et al., 2022). Higher envelope entropy
(parameterized in adiabatic models as an artificially high 1-bar temperature)
also tends to allow more realistic metallicities (Nettelmann et al., 2021;
Miguel et al., 2022) by increasing the mean molecular weight necessary to
achieve the density required by the $J_{2n}$. The source of this tension may
ultimately be a combination of errors in the hydrogen EOS and some degree of
superadiabaticity in the outer envelope.

As in MF21 we find no solutions with uniform $Y^{\prime}$: in all models the
zonal gravity, particularly $J_{4}$, favors a moderately helium-depleted
envelope $Y_{\rm out}^{\prime}=0.15-0.20$, consistent with the measurement by
Koskinen & Guerlet (2018) combining Cassini UVIS and CIRS data but not the
stronger depletion implied by Achterberg & Flasar (2020). Models with $Y_{\rm
out}^{\prime}\gtrsim 0.20$ are generally compatible with gravity but would
require sub-solar $Z_{\rm out}$ and lie outside our prior volume (see Section
2.6). Models with a stronger helium stratification $Y_{\rm out}^{\prime}<0.14$
accommodate more strongly super-solar envelope heavy element enrichments, but
acquire unrealistically high $|J_{4}|$. In light of Saturn’s $9\times$ solar
carbon enrichment, it may be that refinements in the EOS will lead to a
reduced density for hydrogen at depth, allowing for a more metal-rich and
helium-poor outer envelope than the models have achieved to date.

###  3.4 Moment of inertia

Saturn’s moment of inertia is of particular interest from the perspective of a
proposed resonance between Saturn’s spin axis precession and Neptune’s orbital
precession, a commensurability that could help to explain Saturn’s high
obliquity (Ward & Hamilton, 2004; Saillenfest et al., 2021) and the recent
formation of the rings (Wisdom et al., 2022). Our differentially rotating
models have a narrow range of polar moments of inertia $C/MR_{\rm
eq}^{2}=0.2188-0.2191$, consistent with the upper end of estimates from the
gravity modeling from Wisdom et al. (2022), and approximately $0.5\%$ smaller
than the value required for the resonance to be active today. This supports
the notion that the resonance is inactive today, but may have been active in
the past. The significantly higher values estimated by Movshovitz et al.
(2020) may be attributable in part to that work’s assumption that the
magnitude of the $\delta J_{2n}^{\rm winds}$ do not exceed the Galanti & Kaspi
(2017) predictions in advance of the Cassini Grand Finale. The now compelling
evidence for deep zonal winds and their considerable perturbations to zonal
gravity (see Figure 4) warrants a revisiting of the probabilistic Saturn
interior modeling advanced by Movshovitz et al. (2020).

##  4 Discussion

Under our assumed functional form (Equation 2) for the wind decay, we find
that a somewhat modified profile for the surface winds was required for
satisfactory solutions to the gravity field, a requirement made even more
acute when fitting gravity and seismology jointly (Section 3.1). Meanwhile, it
appears that the observed winds (García-Melendo et al., 2011) are robust
across imaging data sets and observational epochs. Other studies have
addressed similar tension by appealing to latitudinally dependent wind depths,
for instance in the Jupiter models of Militzer et al. (2022). We suggest that
Equation 2 is failing to capture some nontrivial spatial dependence below the
cloud level that tends to exaggerate the retrograde rotation of the latitudes
just off the main prograde jet, such as a deeper penetrating wind localized to
those latitudes. Future work exploring these possibilities will benefit from
directly fitting the full set of even zonal harmonics including $J_{8+}$.

Of Saturn’s odd zonal gravity harmonics, at least $J_{3}$, $J_{5}$, and
$J_{9}$ are statistically different from zero (Iess et al., 2019; Jacobson,
2022), a signal arising entirely from the antisymmetric component of the
differential rotation. We have not fit these harmonics directly in our
calculations because our simplified wind model (Equation 2) lacks the
flexibility required to alter these greatly, and more importantly the odd
component of the winds has a negligible effect on the frequency of isolated
Saturn modes (Section 2.5), meaning that gravity remains the best probe of N-S
asymmetry in the differential rotation. Nonetheless we find a plausible fit to
the odd $J_{n}$ (Figure 9) when using the Galanti reconstructed winds (Section
3.1). This is expected since that profile was allowed to deviate from the
García-Melendo profile slightly in order to match the $J_{n}$. We note that
the strict winds yield a much worse fit to the higher degree harmonics (Figure
8), confirming the finding from Galanti et al. (2019) and Galanti & Kaspi
(2021) that the mean cloud level wind profile from García-Melendo et al.
(2011) and a simple decay law like Equation 2 are not a complete description
for the flow within Saturn. Nonetheless the joint analysis of seismology and
gravity clearly strengthens the case for deep winds in Saturn, and the
seismology in particular is an effective anchor for the deep spin period
(Figure 10).

New estimates for Saturn’s gravity moments obtained by Jacobson (2022) depart
from the values reported by Iess et al. (2019), in many cases by an amount
greater than the formal statistical $3\sigma$ uncertainty associated with each
measurement. While the deviations for $J_{2}$, $J_{4}$ and $J_{6}$ are too
small to fundamentally alter our results here, they do suggest that the
$J_{n}$ measurements can be affected by systematic uncertainties that we hope
may be better quantified in the future. In light of this, we encourage
modeling efforts in the meantime to allow for the broadest plausible range of
gravity moments.

Our estimate for the deep spin period $2\pi/\Omega_{\rm S}=633.8-635.5\ {\rm
min}$ is consistent with estimates from atmospheric Rossby waves (Read et al.,
2009) and Saturn’s gravity field and oblateness (Militzer et al., 2019). Our
range of 1.7 min in retrieved spin periods is approximately twice as precise
as the 3.2 min estimate from M19. The improvement follows from the application
of more precise and accurate methods for computing the $J_{2n}$ and the
additional constraint from $J_{6}$ (Section 2.1$-$2.4), as well as more
precise methods for computing mode frequencies and the many additional modes
frequencies used to constrain the model (Section 2.2). That the improvement
was so modest at a factor of just $\approx 2$ follows from the additional
degeneracies that our model encounters, allowing 5 degrees of freedom compared
to M19’s 2 to 3. Our precision is also limited in part by the broad tolerances
we must allow for $J_{2n}$ owing to the imperfect accuracy of the theory of
figures (see Section 2.1 and Figure 9, bottom right panel). It may be
worthwhile to bring more accurate treatments of the rigid body component of
the gravity field, such as the concentric Maclaurin spheroid method (Hubbard,
2013), to bear on this line of seismological inference.

Even more critical are the statistically significant frequency residuals we
find even in our most successful models (Figure 9). These residuals are worst
at low degree $\ell<6$, where Saturn’s less certain deep interior structure
plays a large role (MF21, D21). Here aspects of the spectrum (e.g., the $m=2$
spectrum; see MF21) may be better explained by distinct g-mode cavities
arising from disconnected gradients in heavy element and helium abundance
respectively. More work is required to explore the consequences of such a
structure on the other $\ell<6$ modes; we return to this below. In our model
the worst fitting pattern speed, that of the mode ${}_{4}^{4}f$, is influenced
by an avoided crossing with the mode ${}_{4}^{4}g_{1}$ (MF21, D21, Figure 12),
effectively causing a frequency repulsion between the two (Fuller, 2014).
Resolving this mismatch would require a weaker interaction between the two,
and this may be a useful clue for future work on the deep interior structure.

Our best pattern speed residuals are also significant at $\ell\geq 6$, where
modes are confined to the envelope, indicating outstanding issues with the
thermal, composition, or rotation profile in the outer regions ($r>r_{\rm
out}\sim 0.6$) that are challenging to resolve given the inherent degeneracy
of these different factors. This situation is reminiscent of the persistent
difficulty encountered in fitting Jupiter and Saturn interior models to
gravity while retaining the super-solar envelope heavy element enrichments
indicated by spectroscopic measurements (Guillot et al., 2022). Ad hoc
resolutions have been framed in terms of either (1) unrealistically high 1-bar
temperatures, a simple means of parameterizing some larger interior entropy
than adiabats satisfying these planet’s 1-bar temperatures would predict; or
(2) a systematic error in the density of hydrogen expected from high-pressure
physics, especially at pressures 0.1 to 1 Mbar (Nettelmann et al., 2021;
Miguel et al., 2022; Militzer et al., 2022). Jupiter models proposed by Debras
& Chabrier (2019) jettison these possibilities in favor of a more physically
consistent, albeit intricate, interpretation that appeals to composition and
entropy gradients. Even with the availability of precision Juno and Cassini
gravity the overall picture remains unclear, but composition gradients and
nonadiabaticity will certainly need to be considered as the interpretation of
the seismology matures.

The ring seismology carries more information about low latitudes than high
latitudes, a consequence of the fact that modes with latitudinal wavenumbers
$\ell-m=0,2,4$ are apparently more readily detected than higher wavenumbers.
This is partly a geometric effect: the torques exerted by a resonance decrease
with $\ell-m$ (M19). Preferential excitation of low wavenumbers may also play
a role. French et al. (2021) uncovered the highest latitudinal wavenumber mode
to date in ${}_{19}^{11}f$ (W74.76). We note that this mode fits well for the
rigid model in Figure 2 but is somewhat overestimated in our favored model
including winds (Figure 9). This indicates that in reality this mode must
suffer more self-cancellation in the $\delta\omega$ integrand (Equation 3)
than in our wind model, perhaps suggesting that the prograde rotation at
$|\mu|=0.65-0.85$ is exaggerated in the “Galanti reconstructed” profile.
Future work that incorporates more freedom in the wind model will be able to
use this and all the available constraints spanning $m$ and $\ell-m$, and take
full advantage of the high-degree even ($J_{8+}$; Figure 11) and odd ($J_{3}$,
$J_{5}$, $J_{9}$; Figure 9) gravity harmonics, to refine our understanding of
the flow structure.

The present models, like those in MF21, assume for simplicity that the helium
and heavy element gradients occupy the same region of the interior. In reality
an adiabatic convection zone may separate the stably stratified dilute core
from a shallower, stably stratified, helium gradient. Such a convection zone
would not necessarily fundamentally alter the ${}_{2}^{2}g_{1}$ mode that we
identify with the wave W76.44, but may aid in interpreting other aspects of
the low-degree spectrum, such as the proximity of the $m=2$ pattern W84.64 to
the W87.19/Maxwell pair that MF21 (or our present models; see Figure 12) could
not fully explain. This proximity could be achieved more easily with two
distinct g-mode cavities, where for instance an $n=2$ g-mode confined mostly
to the diffuse core would interact more weakly (and suffer less frequency
repulsion) with the f-mode than in our present model. This multiple-cavity
structure is also attractive from the perspective of Saturn’s magnetic field:
the intervening convective region could be the seat of Saturn’s magnetic
dynamo, with an overlying stably stratified helium gradient helping to explain
the axisymmetry (Stevenson, 1982; Stanley, 2010; Yan & Stanley, 2021) and lack
of secular variability (Moore et al., 2021) of the external field. Future work
focusing on the low-degree ring seismology will be able to gauge the viability
of such a structure.

On a different note, a future mission to the ice giants could lead to the
discovery of Uranian seismicity via resonances its ring system (A’Hearn et
al., 2022), where even a single mode detection would offer a precious anchor
for Uranus’ poorly constrained interior rotation and stratification (Helled et
al., 2010; Movshovitz & Fortney, 2022). The same can be said for Neptune
oscillations, which may be detectable in Neptune’s tenuous ring system, albeit
a more remote possibility. We note that in comparison to Saturn, the ice
giants’ slower fractional rotation $\Omega/\Omega_{\rm dyn}\lesssim 20\%$
implies that frequencies measured in inertial space receive a smaller
contribution from the Doppler shift, resulting in a larger degeneracy between
rotation and the planet’s uncertain interior structure.

##  5 Summary

We have applied oblate, partially stably stratified Saturn interior models to
the wealth of gravity and ring seismology data made available by Cassini, to
conclude that

  1. 1.

Saturn’s cloud-level winds penetrate deeply before decaying rapidly at a depth
of $0.125-0.138\ R_{\rm eq}$, or $7,530-8,320$ km, with a median depth $0.132\
R_{\rm eq}=7,960$ km. These are consistent with estimates from gravity and
magnetic field science (Galanti et al., 2019; Galanti & Kaspi, 2021).

  2. 2.

Minor changes to the García-Melendo et al. (2011) cloud-level wind profile
appear to be necessary, likely indicating that Saturn’s wind decay is more
complicated than the functional form we have assumed in Equation 2.

  3. 3.

Saturn’s rigidly rotating deep interior has a spin period of $633.8-635.5$ min
($10^{\rm h}33^{\rm m}48^{\rm s}-10^{\rm h}35^{\rm m}30^{\rm s}$), or a
pattern speed $815.7-818.0\ {\rm deg\ d}^{-1}$. The median spin period is
634.7 min ($10^{\rm h}34^{\rm m}48^{\rm s}$) for a pattern speed $816.8\ {\rm
deg\ d}^{-1}$. The new estimate is consistent with Mankovich et al. (2019) but
more precise thanks to additional ring seismology constraints and improved
methods for calculating the gravity and normal modes, both in the rigidly
rotating limit and in the presence of differential rotation.

Despite these gains, the persistent pattern speed residuals seen for our most
successful model, at both low and high angular degree (Figure 9), show that
the Cassini ring seismology data are still not perfectly understood. Future
analysis of the same data but allowing for more freedom in the composition,
temperature, and rotation profiles should bring Saturn’s interior into clearer
focus.

## Acknowledgments

We thank the anonymous referees for their helpful critiques and Eli Galanti
for sharing his Saturn wind solution. We are grateful for support from the
Caltech Center for Comparative Planetary Evolution. JWF gratefully
acknowledges support from the Sloan Foundation through grant FG-2018-10515,
and JWD acknowledges support from the Natural Sciences and Engineering
Research Council of Canada (NSERC) [funding reference #CITA 490888-16].

## References

  * Achterberg & Flasar (2020) Achterberg, R. K., & Flasar, F. M. 2020, , 1, 30, doi: 10.3847/PSJ/ab9cb6
  * Aerts et al. (2010) Aerts, C., Christensen-Dalsgaard, J., & Kurtz, D. W. 2010, Asteroseismology (Dordrecht: Springer), doi: 10.1007/978-1-4020-5803-5
  * A’Hearn et al. (2022) A’Hearn, J. A., Hedman, M. M., Mankovich, C. R., Aramona, H., & Marley, M. S. 2022, , 3, 194, doi: 10.3847/PSJ/ac82bb
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481, doi: 10.1146/annurev.astro.46.060407.145222
  * Baillié et al. (2011) Baillié, K., Colwell, J. E., Lissauer, J. J., Esposito, L. W., & Sremčević, M. 2011, Icarus, 216, 292, doi: 10.1016/j.icarus.2011.05.019
  * Ballot et al. (2012) Ballot, J., Lignières, F., Prat, V., Reese, D. R., & Rieutord, M. 2012, in Astronomical Society of the Pacific Conference Series, Vol. 462, Progress in Solar/Stellar Physics with Helio- and Asteroseismology, ed. H. Shibahashi, M. Takata, & A. E. Lynas-Gray, 389.  https://arxiv.org/abs/1109.6856
  * Bonazzola et al. (1998) Bonazzola, S., Gourgoulhon, E., & Marck, J.-A. 1998, Phys. Rev. D, 58, 104020, doi: 10.1103/PhysRevD.58.104020
  * Brygoo et al. (2021) Brygoo, S., Loubeyre, P., Millot, M., et al. 2021, Nature, 593, 517, doi: 10.1038/s41586-021-03516-0
  * Cao et al. (2020) Cao, H., Dougherty, M. K., Hunt, G. J., et al. 2020, Icarus, 344, 113541, doi: 10.1016/j.icarus.2019.113541
  * Cao & Stevenson (2017a) Cao, H., & Stevenson, D. J. 2017a, Journal of Geophysical Research (Planets), 122, 686, doi: 10.1002/2017JE005272
  * Cao & Stevenson (2017b) —. 2017b, Icarus, 296, 59, doi: 10.1016/j.icarus.2017.05.015
  * Debras & Chabrier (2019) Debras, F., & Chabrier, G. 2019, ApJ, 872, 100, doi: 10.3847/1538-4357/aaff65
  * Desch & Kaiser (1981) Desch, M. D., & Kaiser, M. L. 1981, Geophys. Res. Lett., 8, 253, doi: 10.1029/GL008i003p00253
  * Dewberry et al. (2022) Dewberry, J., Mankovich, C. R., & Fuller, J. 2022, MNRAS, 516, 1, doi: 10.1093/mnras/stac1957
  * Dewberry et al. (2021) Dewberry, J. W., Mankovich, C. R., Fuller, J., Lai, D., & Xu, W. 2021, , 2, 198, doi: 10.3847/PSJ/ac0e2a
  * Durante et al. (2022) Durante, D., Guillot, T., Iess, L., et al. 2022, Nature Communications, 13, 4632, doi: 10.1038/s41467-022-32299-9
  * Fletcher et al. (2009) Fletcher, L. N., Orton, G. S., Teanby, N. A., Irwin, P. G. J., & Bjoraker, G. L. 2009, Icarus, 199, 351, doi: 10.1016/j.icarus.2008.09.019
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * French et al. (2021) French, R. G., Bridges, B., Hedman, M. M., et al. 2021, Icarus, 370, 114660, doi: 10.1016/j.icarus.2021.114660
  * French et al. (2019) French, R. G., McGhee-French, C. A., Nicholson, P. D., & Hedman, M. M. 2019, Icarus, 319, 599, doi: 10.1016/j.icarus.2018.10.013
  * French et al. (2016) French, R. G., Nicholson, P. D., Hedman, M. M., et al. 2016, Icarus, 279, 62, doi: 10.1016/j.icarus.2015.08.020
  * Fuller (2014) Fuller, J. 2014, Icarus, 242, 283, doi: 10.1016/j.icarus.2014.08.006
  * Galanti & Kaspi (2017) Galanti, E., & Kaspi, Y. 2017, ApJ, 843, L25, doi: 10.3847/2041-8213/aa7aec
  * Galanti & Kaspi (2021) —. 2021, MNRAS, 501, 2352, doi: 10.1093/mnras/staa3722
  * Galanti et al. (2019) Galanti, E., Kaspi, Y., Miguel, Y., et al. 2019, Geophys. Res. Lett., 46, 616, doi: 10.1029/2018GL078087
  * García-Melendo et al. (2011) García-Melendo, E., Pérez-Hoyos, S., Sánchez-Lavega, A., & Hueso, R. 2011, Icarus, 215, 62, doi: 10.1016/j.icarus.2011.07.005
  * Giampieri et al. (2006) Giampieri, G., Dougherty, M. K., Smith, E. J., & Russell, C. T. 2006, Nature, 441, 62, doi: 10.1038/nature04750
  * Guillot (1999) Guillot, T. 1999, Planet. Space Sci., 47, 1183, doi: 10.1016/S0032-0633(99)00043-4
  * Guillot et al. (2022) Guillot, T., Fletcher, L. N., Helled, R., et al. 2022, arXiv e-prints, arXiv:2205.04100.  https://arxiv.org/abs/2205.04100
  * Guillot et al. (2018) Guillot, T., Miguel, Y., Militzer, B., et al. 2018, Nature, 555, 227, doi: 10.1038/nature25775
  * Hedman & Nicholson (2013) Hedman, M. M., & Nicholson, P. D. 2013, AJ, 146, 12, doi: 10.1088/0004-6256/146/1/12
  * Hedman & Nicholson (2014) —. 2014, MNRAS, 444, 1369, doi: 10.1093/mnras/stu1503
  * Hedman et al. (2022) Hedman, M. M., Nicholson, P. D., El Moutamid, M., & Smotherman, S. 2022, , 3, 61, doi: 10.3847/PSJ/ac4df8
  * Hedman et al. (2019) Hedman, M. M., Nicholson, P. D., & French, R. G. 2019, AJ, 157, 18, doi: 10.3847/1538-3881/aaf0a6
  * Helled et al. (2010) Helled, R., Anderson, J. D., & Schubert, G. 2010, Icarus, 210, 446, doi: 10.1016/j.icarus.2010.06.037
  * Helled et al. (2015) Helled, R., Galanti, E., & Kaspi, Y. 2015, Nature, 520, 202, doi: 10.1038/nature14278
  * Helled et al. (2020) Helled, R., Mazzola, G., & Redmer, R. 2020, Nature Reviews Physics, 2, 562, doi: 10.1038/s42254-020-0223-3
  * Hubbard (2013) Hubbard, W. B. 2013, ApJ, 768, 43, doi: 10.1088/0004-637X/768/1/43
  * Iess et al. (2019) Iess, L., Militzer, B., Kaspi, Y., et al. 2019, Science, 364, aat2965, doi: 10.1126/science.aat2965
  * Jacobson (2022) Jacobson, R. A. 2022, AJ, 164, 199, doi: 10.3847/1538-3881/ac90c9
  * Kaspi et al. (2018) Kaspi, Y., Galanti, E., Hubbard, W. B., et al. 2018, Nature, 555, 223, doi: 10.1038/nature25793
  * Koskinen & Guerlet (2018) Koskinen, T. T., & Guerlet, S. 2018, Icarus, 307, 161, doi: 10.1016/j.icarus.2018.02.020
  * Lindal et al. (1985) Lindal, G. F., Sweetnam, D. N., & Eshleman, V. R. 1985, AJ, 90, 1136, doi: 10.1086/113820
  * Liu et al. (2008) Liu, J., Goldreich, P. M., & Stevenson, D. J. 2008, Icarus, 196, 653, doi: 10.1016/j.icarus.2007.11.036
  * Lodders (2003) Lodders, K. 2003, ApJ, 591, 1220, doi: 10.1086/375492
  * Mankovich et al. (2019) Mankovich, C., Marley, M. S., Fortney, J. J., & Movshovitz, N. 2019, ApJ, 871, 1, doi: 10.3847/1538-4357/aaf798
  * Mankovich & Fuller (2021) Mankovich, C. R., & Fuller, J. 2021, Nature Astronomy, 5, 1103, doi: 10.1038/s41550-021-01448-3
  * Markham et al. (2020) Markham, S., Durante, D., Iess, L., & Stevenson, D. 2020, , 1, 27, doi: 10.3847/PSJ/ab9f21
  * Marley (1991) Marley, M. S. 1991, Icarus, 94, 420, doi: 10.1016/0019-1035(91)90239-P
  * Marley & Porco (1993) Marley, M. S., & Porco, C. C. 1993, Icarus, 106, 508, doi: 10.1006/icar.1993.1189
  * Miguel et al. (2016) Miguel, Y., Guillot, T., & Fayon, L. 2016, A&A, 596, A114, doi: 10.1051/0004-6361/201629732
  * Miguel et al. (2022) Miguel, Y., Bazot, M., Guillot, T., et al. 2022, A&A, 662, A18, doi: 10.1051/0004-6361/202243207
  * Militzer & Hubbard (2013) Militzer, B., & Hubbard, W. B. 2013, ApJ, 774, 148, doi: 10.1088/0004-637X/774/2/148
  * Militzer et al. (2019) Militzer, B., Wahl, S., & Hubbard, W. B. 2019, ApJ, 879, 78, doi: 10.3847/1538-4357/ab23f0
  * Militzer et al. (2022) Militzer, B., Hubbard, W. B., Wahl, S., et al. 2022, , 3, 185, doi: 10.3847/PSJ/ac7ec8
  * Moore et al. (2021) Moore, K. M., Bolton, B., Cao, H., K. Dougherty, M., & Bloxham, J. 2021, , 2, 181, doi: 10.3847/PSJ/ac173c
  * Movshovitz & Fortney (2022) Movshovitz, N., & Fortney, J. J. 2022, , 3, 88, doi: 10.3847/PSJ/ac60ff
  * Movshovitz et al. (2020) Movshovitz, N., Fortney, J. J., Mankovich, C., Thorngren, D., & Helled, R. 2020, ApJ, 891, 109, doi: 10.3847/1538-4357/ab71ff
  * Nettelmann et al. (2021) Nettelmann, N., Movshovitz, N., Ni, D., et al. 2021, , 2, 241, doi: 10.3847/PSJ/ac390a
  * Ouazzani et al. (2012) Ouazzani, R. M., Dupret, M. A., & Reese, D. R. 2012, A&A, 547, A75, doi: 10.1051/0004-6361/201219548
  * Read et al. (2009) Read, P. L., Dowling, T. E., & Schubert, G. 2009, Nature, 460, 608, doi: 10.1038/nature08194
  * Reese et al. (2006) Reese, D., Lignières, F., & Rieutord, M. 2006, A&A, 455, 621, doi: 10.1051/0004-6361:20065269
  * Rieutord et al. (2016) Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, Journal of Computational Physics, 318, 277, doi: 10.1016/j.jcp.2016.05.011
  * Rosen et al. (1991) Rosen, P. A., Tyler, G. L., Marouf, E. A., & Lissauer, J. J. 1991, Icarus, 93, 25, doi: 10.1016/0019-1035(91)90161-L
  * Saillenfest et al. (2021) Saillenfest, M., Lari, G., & Boué, G. 2021, Nature Astronomy, 5, 345, doi: 10.1038/s41550-020-01284-x
  * Stanley (2010) Stanley, S. 2010, Geophys. Res. Lett., 37, L05201, doi: 10.1029/2009GL041752
  * Stevenson (1982) Stevenson, D. J. 1982, Geophysical and Astrophysical Fluid Dynamics, 21, 113, doi: 10.1080/03091928208209008
  * Stevenson (2020) —. 2020, Annual Review of Earth and Planetary Sciences, 48, 465, doi: 10.1146/annurev-earth-081619-052855
  * Takata & Saio (2013) Takata, M., & Saio, H. 2013, PASJ, 65, 68, doi: 10.1093/pasj/65.3.68
  * Thompson (1990) Thompson, S. L. 1990, ANEOS–—Analytic Equations of State for Shock Physics Codes, Sandia Natl. Lab. Doc. SAND89-2951, http://prod.sandia.gov/techlib/access-control.cgi/1989/892951.pdf
  * Wahl et al. (2017) Wahl, S. M., Hubbard, W. B., Militzer, B., et al. 2017, Geophys. Res. Lett., 44, 4649, doi: 10.1002/2017GL073160
  * Ward & Hamilton (2004) Ward, W. R., & Hamilton, D. P. 2004, AJ, 128, 2501, doi: 10.1086/424533
  * Wisdom et al. (2022) Wisdom, J., Dbouk, R., Militzer, B., et al. 2022, Science, 377, 1285, doi: 10.1126/science.abn1234
  * Wisdom & Hubbard (2016) Wisdom, J., & Hubbard, W. B. 2016, Icarus, 267, 315, doi: 10.1016/j.icarus.2015.12.030
  * Xu & Lai (2017) Xu, W., & Lai, D. 2017, Phys. Rev. D, 96, 083005, doi: 10.1103/PhysRevD.96.083005
  * Yan & Stanley (2021) Yan, C., & Stanley, S. 2021, AGU Advances, 2, e00318, doi: 10.1029/2020AV000318
  * Ye et al. (2018) Ye, S. Y., Fischer, G., Kurth, W. S., Menietti, J. D., & Gurnett, D. A. 2018, Geophys. Res. Lett., 45, 7297, doi: 10.1029/2018GL077976
  * Zhang et al. (2015) Zhang, K., Kong, D., & Schubert, G. 2015, ApJ, 806, 270, doi: 10.1088/0004-637X/806/2/270
  * Zharkov & Trubitsyn (1978) Zharkov, V. N., & Trubitsyn, V. P. 1978, Physics of Planetary Interiors (Tucson, AZ: Pachart)

##  Appendix A Numerical method for calculating $\delta J_{n}$

Consider a rigidly rotating equilibrium state characterized by a density
$\rho_{0}$, pressure $P_{0},$ velocity field ${\bf
u}_{0}={\bf\Omega}_{0}\times{\bf r}$ (${\bf\Omega}_{0}=\Omega_{0}\hat{\bf
z}$), and gravity field $\Phi_{0}$ that satisfy the steady-state Euler and
Poisson equations. The linearized equation of motion for time-independent,
axisymmetric perturbations $\rho^{\prime}(r,\theta)$, $P^{\prime}(r,\theta)$,
${\bf u}=r\sin\theta\delta\Omega(r,\theta)\hat{\boldsymbol{\phi}}$,
$\Phi^{\prime}(r,\theta)$ to this rigidly rotating equilibrium is

${\bf G}\rho^{\prime}-\nabla
P^{\prime}-\rho_{0}\nabla\Phi^{\prime}=2\boldsymbol{\Omega}_{0}\times(\rho_{0}{\bf
u}),$ (A1)

where ${\bf G}=\rho_{0}^{-1}\nabla P_{0}$ is the “effective gravity”
associated with the rigidly rotating state. The barotropic nature of the
rigidly rotating equilibrium implies that $\nabla\times{\bf G}=0,$ so taking
the curl of Equation A1 produces

${\bf
G}\times\nabla\rho^{\prime}+\nabla\rho_{0}\times\nabla\Phi^{\prime}=2{\bf\Omega}_{0}\cdot\nabla(\rho_{0}{\bf
u}).$ (A2)

$\rho^{\prime}$ and $\Phi^{\prime}$ are additionally related by
$\nabla^{2}\Phi^{\prime}=4\pi G\rho^{\prime}$. Equation A2 can be expressed on
the covariant basis of an arbitrary curvilinear coordinate system as

$\epsilon^{ijk}[G_{j}\partial_{k}\rho^{\prime}+\partial_{j}\rho_{0}\partial_{k}\Phi^{\prime}]=2\Omega_{0}^{j}[\partial_{j}(\rho_{0}u^{i})+\rho_{0}\Gamma^{i}_{kj}u^{k}],$
(A3)

where paired upper and lower indices denote summation, $\partial_{i}$ denotes
partial differentiation with respect to coordinate $x^{i}$,
$\epsilon^{ijk}=J^{-1}[i,j,k]$ for $J$ the Jacobian of the coordinate system,
and $\Gamma^{i}_{kj}$ are Christoffel symbols.

As in our mode calculations (see subsection 2.2), we adopt the non-orthogonal
coordinate system $(\zeta,\theta,\phi)$ described in D21 (and originally
proposed by Bonazzola et al., 1998), which involves a quasi-radial coordinate
$\zeta$ constructed so that $\zeta=1$ on the pre-computed surface of the
rigidly rotating planet. Retaining nonzero geometric factors (summarized in
Rieutord et al., 2016), we find

$(G_{\zeta}\partial_{\theta}-G_{\theta}\partial_{\zeta})\rho^{\prime}+(\partial_{\zeta}\rho_{0}\partial_{\theta}-\partial_{\theta}\rho_{0}\partial_{\zeta})\Phi^{\prime}=2\Omega_{0}r\sin\theta\left[(r\cos\theta+r_{\theta}\sin\theta)\partial_{\zeta}-r_{\zeta}\sin\theta\partial_{\theta}\right](\rho_{0}\delta\Omega):=f(r,\theta),$
(A4)

where $G_{i}=\rho_{0}^{-1}\partial_{i}P_{0},$ and $r_{i}=\partial_{i}r$ (note
that $r=r(\zeta,\theta)$).

Inserting expansions
$\rho^{\prime}=\sum_{n}\rho^{n}(\zeta)Y_{n}^{0}(\theta,\phi)$,
$\Phi^{\prime}=\sum_{n}\Phi^{n}(\zeta)Y_{n}^{0}(\theta,\phi)$ in zonal
spherical harmonics, we project onto an arbitrary degree $n=\ell$ by taking
the inner product with $\partial_{\theta}Y_{\ell}^{0}$. This leads to an
infinite set of coupled ordinary differential equations in $\zeta$:

$\sum_{n}\left[\langle\partial_{\theta}Y_{\ell}^{0},\partial_{\theta}Y_{n}^{0}G_{\zeta}\rangle-\langle\partial_{\theta}Y_{\ell}^{0},Y_{n}^{0}G_{\theta}\rangle\partial_{\zeta}\right]\rho^{n}+\sum_{n}\left[\langle\partial_{\theta}Y_{\ell}^{0},\partial_{\theta}Y_{n}^{0}\partial_{\zeta}\rho_{0}\rangle-\langle\partial_{\theta}Y_{\ell}^{0},Y_{n}^{0}\partial_{\theta}\rho_{0}\rangle\partial_{\zeta}\right]\Phi^{n}=\langle\partial_{\theta}Y_{\ell}^{0},f\rangle,$
(A5)

where
$\langle\partial_{\theta}Y_{\ell}^{0},h\rangle=\int_{0}^{\pi}\int_{0}^{2\pi}h\partial_{\theta}Y_{\ell}^{0}\sin\theta\text{d}\phi\text{d}\theta$
for a given function $h(r,\theta)$.

We solve Equation A2 simultaneously with Poisson’s equation (the projected
form of which is given by eq. A32 in D21), using the same Chebyshev
collocation method applied to mode calculations. For boundary conditions, we
use boundary bordering to enforce a vanishing gravitational perturbation at
infinity. We additionally impose regularity boundary conditions on both
$\rho^{\prime}$ and $\Phi^{\prime}$ at the origin.

Lastly, we note that Equation A5 does not constrain $\rho^{0},$ since the
integrals on the lefthand side vanish for $\ell=0$. In a spherical,
nonrotating planet, this would have no effect on the gravitational moments
(e.g., Zhang et al., 2015); in a centrifugally flattened planet this is no
longer the case, but we find that the higher-degree $\delta J_{n}$ are still
largely insensitive to the prescription for $\rho^{0},$ provided that we
ensure $\delta J_{0}=0.$ Since the gravitational moments are directly related
to the gravitational perturbation in the external vacuum (e.g., Dewberry et
al., 2022), including self-gravity allows for convenient application of this
additional constraint.

The lefthand panels in Figure 14 compare the results of calculations under
various approximations against $\delta J_{n}$ reported by Wisdom & Hubbard
(2016) for their “DR3” differential rotation profile. The top panel shows raw
$\delta J_{n}$ calculations, with negative values indicated by empty points,
while the bottom panel shows residuals from the calculations of Wisdom &
Hubbard (2016). Blue points indicate results from self-consistent polytropic
model calculations of Dewberry et al. (2022), and orange points show the
result of calculations solving the oblate TGWE (Equation A5) in full (the
approach taken in this work). Meanwhile green points show calculations in
which we ignore the gravitational perturbation in Equation A5, solving the TWE
while still including the geometric effects associated with planetary
oblateness. Lastly, the red and purple points respectively indicate solutions
of the TGWE and TWE in a spherical geometry. The righthand panels in Figure 14
show analogous results for one of the differential rotation profiles
considered by Dewberry et al. (2022) (adopting a value $d=0.3$ for their depth
parameter). The calculations in the lefthand (righthand) panels were performed
with $N_{\zeta}=80$ ($90$) collocation points in the quasiradial direction,
and a spherical harmonic truncation of $N_{\ell}=30$ ($60$).

In agreement with Cao & Stevenson (2017a), Figure 14 suggests that the
planetary oblateness is more important to include than self-gravity for the
higher degree $\delta J_{n}$, and that $\delta J_{2}$ and $\delta J_{4}$
computed from the thermal wind equation under any approximation should be
treated with caution. As noted by Cao & Stevenson (2017a), however,
uncertainties in $\delta J_{2}$ and $\delta J_{4}$ due to winds are
insignificant compared with uncertainties resulting from uncertain composition
parameters or equations of state.

Figure 14: Gravitational coefficient shifts $\delta J_{n}$ (top) and
associated residuals (bottom), computed for the “DR3” differential rotation
profile considered by Wisdom & Hubbard (2016) (left) and one of the
differential rotation profiles from Dewberry et al. (2022) (right).
Differently colored points show $\delta J_{n}$ calculated directly from the
Euler and Poisson equations (black, blue), and from the TGWE/TWE under various
approximations (orange, green, red, purple).

##  Appendix B Effect of omitting the $m=3,4,5$ constraints

The largest tension between our favored model and the data is in the pattern
speeds of the $m=3,4,5$ Lindblad resonances (Figure 9) that we attribute to
sectoral f-modes. These deep-seated modes are highly sensitive to the Brunt-
Väisälä frequency, and hence to the composition profile, within the deep
interior. Consequently their large pattern speed residuals ($\approx 1-6\ {\rm
deg\ d}^{-1}$) likely stem from the model’s simplistic composition profile.

Figure 15:  As in Figure 9, but for a sample neglecting the $m=3,4,5$
constraints. Here the central panel includes vertical resonances (crosses) in
addition to Lindblad resonances (circles).

To ensure that these less well understood modes, particularly the outlying
$m=4$ mode, do not bias our results, we create a sample that neglects the
$m=3,4,5$ modes but is otherwise identical to the main “precise seismology”
calculation in Figures 9-13. The best model in the alternate sample is
summarized in Figure 15.

The exclusion of the $\ell=m=4$ mode in particular, which all of our models
strongly underestimate, allows this model to better fit the remaining modes by
adopting a modestly faster deep spin (634.0 min versus Figure 9’s 634.8 min)
and deeper wind decay depth ($0.140\,R_{\rm eq}$ versus $0.132\,R_{\rm eq})$.
These differences are within the posterior distributions for each parameter
from our main calculation (Figure 13). The systematic offsets as a function of
$m$ or $\ell-m$ remain virtually unchanged. We conclude that the poorly fit
$m=3,4,5$ modes do not fundamentally bias our main results, but we urge that
future work prioritize these low-$m$ modes as a window into Saturn’s still
uncertain composition profile.