Constraining Saturn's interior with ring seismology: effects of differential rotation and stable stratification

Normal mode oscillations in Saturn excite density and bending waves in the C Ring, providing a valuable window into the planet's interior. Saturn's fundamental modes (f modes) excite the majority of the observed waves, while gravito-inertial modes (rotationally modified g modes) associated with stable stratification in the deep interior provide a compelling explanation for additional density waves with low azimuthal wavenumbers m. However, multiplets of density waves with nearly degenerate frequencies, including an m=3 triplet, still lack a definitive explanation. We investigate the effects of rapid and differential rotation on Saturn's oscillations, calculating normal modes for independently constrained interior models. We use a non-perturbative treatment of rotation that captures the full effects of the Coriolis and centrifugal forces, and consequently the mixing of sectoral f modes with g modes characterized by very different spherical harmonic degrees. Realistic profiles for differential rotation associated with Saturn's zonal winds can enhance these mode interactions, producing detectable oscillations with frequencies separated by less than 1%. Our calculations demonstrate that a three-mode interaction involving an f mode and two g modes can feasibly explain the finely split m=3 triplet, although the fine-tuning required to produce such an interaction generally worsens agreement with seismological constraints provided by m=2 density waves. Our calculations additionally demonstrate that sectoral f mode frequencies are measurably sensitive to differential rotation in Saturn's convective envelope. Finally, we find that including realistic equatorial antisymmetry in Saturn's differential rotation profile couples modes with even and odd equatorial parity, producing oscillations that could in principle excite both density and bending waves simultaneously.


INTRODUCTION
The interiors of gas giants remain enigmatic, even within our own solar system. In particular, we still know very little about the extent and composition of Saturn's core. Cassini 's Grand Finale facilitated pre-cise measurement of the planet's non-spherical gravitational field (Iess et al. 2019), but even knowledge of highorder harmonic coefficients (through J 12 ) leaves room for significant degeneracy in deep-interior density profiles (Movshovitz et al. 2020). Uncertainties surrounding the equation of state for hydrogen and helium mixtures at pressures beyond the reach of laboratory experiments, not to mention the concentration and dissolution of heavier elements, complicate the process of bridging this observational gap (e.g., Helled 2018).
Fortunately, Saturn's rings present a rare opportunity to place seismological constraints on the planet's interior. Analyzing occultation data obtained by Cassini, Hedman & Nicholson (2013) validated a decades-old prediction (Marley 1991;Marley & Porco 1993) that Saturn's "fundamental" oscillation modes (f modes) excite density waves at Lindblad resonances in the C Ring. Continued analysis of Cassini observations has revealed a wealth of additional density (and bending) waves, thought to be excited by planetary oscillations with even (odd) equatorial symmetry (Hedman & Nicholson 2014;French et al. 2016French et al. , 2019Hedman et al. 2019).
Saturn ring seismology has already provided insight into the planet's internal structure. Recently, Mankovich et al. (2019) used f mode identifications to measure Saturn's mean internal rotation rate to a high degree of accuracy (a task previously made difficult by the axisymmetry and nearly exact polar alignment of Saturn's magnetic field; Cao et al. 2020). One of the earliest surprises to come out of Saturn ring seismology was a strong indication that the planet is not, as historically thought, composed simply of a solid core and a convective envelope (but see, e.g., Leconte & Chabrier 2013): Fuller (2014) showed that observations of multiple density waves with the same azimuthal wavenumber m and nearly degenerate frequencies could be explained by the rotational mixing of f modes with gravito-inertial modes (hereafter g modes) associated with stable stratification in the deep interior.
Based on this analysis, Fuller (2014) predicted the appearance of additional density waves driven by lowradial order g modes. French et al. (2016French et al. ( , 2019 subsequently observed density waves identifiable with such g modes, which Mankovich & Fuller (2021) recently fit jointly with Saturn's gravity field to find that Saturn's core is likely very diffuse, and the region of stable stratification associated with composition gradients very extensive. However, the calculations of Fuller (2014) and Mankovich & Fuller (2021) encounter difficulty in reproducing the finest frequency splittings of less than 1% observed for a triplet of m = 3 density waves (Hedman & Nicholson 2013). This difficulty could be related to the authors' method for including the effects of rotation, which they approximated as rigid-body and treated using second-order perturbation theory.
In this paper, we revisit both of these approximations. We first calculate oscillations from the rigidly rotating Saturn models of Mankovich & Fuller (2021) without making any approximation with respect to the modes' modification by either the Coriolis force or the centrifugal distortion of the equilibrium planetary structure. We then investigate the additional coupling of these pulsa-tions by differential rotation associated with Saturn's zonal winds.
Along with conventional (but rotationally modified) g modes, these calculations reveal sequences of so-called "rosette" modes (Ballot et al. 2012;Takata & Saio 2013) that can only be described with higher-order treatments of rotation. Like the relatively low-degree g modes involved in the rotational mixing considered by Fuller (2014), both rosette modes and high-degree g modes gain enhanced surface gravitational perturbations from avoided crossings with f modes. The density of the high-g mode spectrum in frequency-space presents numerous opportunities for such mixing, but the avoided crossings are proportionately narrower. Rigidly rotating models therefore require fine-tuning to produce multiple detectable modes with finely split frequencies.
Limited to the outer envelope, differential rotation from Saturn's zonal winds primarily affects the f modes in our calculations. The zonal winds produce small, but significant shifts in the frequencies of high-degree sectoral ( ∼ m) f modes. Varying the depth of wind decay, we find that differential rotation can enhance rotational mixing between the low-sectoral f modes and the underlying dense spectrum of high-degree g modes. We present examples of three-mode interactions that could in principle explain the observed m = 3 triplet of density waves with finely split frequencies, although producing these interactions still requires fine-tuning (which generally worsens the agreement of our calculations of other low-order modes with observations).
Finally, we find that realistic equatorial antisymmetry in the differential rotation profile couples oscillation modes with even and odd equatorial parity. This coupling is weak, but sufficient to produce equatorially asymmetric oscillations that could simultaneously excite density and bending waves at separate locations in the rings. Observational confirmation of such simultaneous density/bending wave excitation would provide a very useful constraint on Saturn's interior, since the weak antisymmetric component of Saturn's zonal winds is likely only capable of mixing f modes with a limited subset of opposite-parity g modes of low degree and low order. This paper takes the following structure: Section 2 and Section 3 introduce our models of Saturn and our method for calculating their associated oscillation modes, respectively. Section 4 and Section 5 then present the results of our calculations with purely rigid and differential rotation (resp.), Section 6 provides discussion, and Section 7 summarizes our conclusions. The Appendices provide further details about our numerical methods.

Rigidly rotating models
We start from the rigidly rotating, oblate models of Mankovich & Fuller (2021). Calculated using a fourthorder theory of figures (Nettelmann 2017), these parameterize Saturn's interior structure via a smooth variation in heavy element and Helium mass fractions over a core-envelope transition region. The models are thermally adiabatically stratified and adopt a temperature of 135 K at P = 1 bar. We assume an equatorial radius R S = 60.268×10 3 km, GM S = 37, 931, 207.7km 3 s −2 (Jacobson et al. 2006), and a 10.561hr period (Mankovich et al. 2019) corresponding to a bulk rotation rate Ω S 0.397Ω dyn , where Ω dyn = GM S /R 3 S is the equatorial Keplerian rotation rate. Figure 1 plots profiles of density (top) and Brunt-Väisälä frequency (bottom) for the models considered in this paper, which have stably stratified regions extending from r = 0 to r stab = 0.65R S , 0.70R S , and 0.72R S at the equator. We concentrate on these three because they illustrate different aspects of normal mode coupling by differential rotation, and take the model with r stab = 0.70R S as fiducial because it agrees most closely with the best-fit found by Mankovich & Fuller (2021) for smooth composition gradients. This work focuses on exploring the potential effects of differential rotation on oscillations in Saturn, not providing an exact description for the planet's interior. However, we do match the rigid-body component of the observed zonal gravity harmonic J 2 to within 1%, and J 4 , J 6 , J 8 to within 5% (Iess et al. 2019;Galanti & Kaspi 2021).

Zonal winds
Our oscillation calculations directly incorporate Saturn's observed wind profiles, as reported by García-Melendo et al. (2011) and adjusted for the rotation period of 10.561hr constrained by Mankovich et al. (2019). Figure 2 (top left) plots the zonal wind speed with µ = cos θ (here θ is the colatitude), which we smooth and relax to zero at the poles by expanding the surface azimuthal velocity as U = w 60 l=0 c l P l (µ). Here c l are coefficients calculated by projecting onto Legendre polynomials P l , and w is a window function that transitions smoothly to zero for |µ| > 0.975.
Informed by recent investigations combining magnetic and gravity measurements (Galanti & Kaspi 2021), we assume that the zonal wind penetrates deeply as barotropic "rotation on cylinders," producing a perturbation to the angular velocity with the form RΩ D = ηU. Here R = r sin θ is the distance to the rotation axis (cylindrical radius), and η defines a decay function. We adopt the simplified relation where r s (µ) tracks the oblate planetary surface, while ∆ and d parameterize the width and depth of decay to rigid rotation in the interior. We take ∆ = 0.02R s , and vary the decay depth. Figure 2 (bottom) plots an example decay function with a fiducial value of d = 0.125R S 7.5 × 10 3 km, which produces an angular velocity perturbation Ω D with the meridional profile shown in Figure 2 (top right), and a maximum deviation from the bulk rotation rate of 3% at the equator. As demonstrated by the orange dotted line in Figure 2  odd − m), which would otherwise be unable to interact (see Section 5). Unless otherwise stated, we therefore include odd-degree Legendre polynomials in our smoothing expansion, but we note that antisymmetric surfacelevel winds introduce discontinuities when allowed to extend barotropically to equatorial radii r(µ = 0) < R S . We eliminate these discontinuities in Ω D (and more importantly its gradient) by smoothly transitioning the antisymmetric component of the surface wind to zero within |µ| < 0.5, using a window function composed of tanh functions with a scale length set to half the polar grid size. The antisymmetric component of the wind is already small for |µ| 0.5, and our results are insensitive to this choice of window function.

OSCILLATION CALCULATIONS
This section outlines our methods for calculating oscillation modes. We use a combination of non-perturbative (Section 3.1) and perturbative (Section 3.2) techniques for including the effects of rigid and differential rotation (resp.), since Saturn's bulk rotation rate of Ω S 0.397Ω dyn is rapid, but the differential rotation associated with zonal winds constitutes a relatively small perturbation. Throughout this section, we refer to both spherical and cylindrical polar coordinate systems (r, θ, φ) and (R, φ, z).

Non-perturbative calculations with rigid rotation
In a frame rotating with constant angular velocity Ω S = Ω Sẑ , the equations governing a self-gravitating, inviscid gas characterized by velocity u, pressure P , den-sity ρ and gravitational field Φ are (e.g., Lynden-Bell & Ostriker 1967) where D t = ∂ t + u · ∇ denotes the convective derivative, and Γ 1 is the first adiabatic exponent. Here we have neglected non-adiabatic heating in the thermal energy equation (Equation 4). Perturbing these equations around a rigidly rotating, barotropic equilibrium state with pressure P 0 and density ρ 0 at rest in the rotating frame, the linearized equations governing adiabatic oscillations with the harmonic time dependence exp[−iωt], where ω is a (rotating frame) frequency, can be written as Here G = ρ −1 0 ∇P 0 is the effective gravity (including both gravitational and centrifugal accelerations), c 2 A = Γ 1 P 0 /ρ 0 is the square of the adiabatic sound speed, and v is an Eulerian velocity perturbation. Lastly, with primes denoting Eulerian perturbations to the pressure P , density ρ and gravitational field Φ , we have defined the variables h = P /ρ 0 and b = ρ /ρ 0 .
Saturn's rapid bulk rotation breaks the spherical symmetry of the oscillation equations, through both the introduction of the Coriolis force 2Ω S × v, and a centrifugal flattening of the equilibrium state. To account for these effects, both of which can be important for the low-order f-modes and g-modes of interest, we use a non-perturbative numerical method (e.g., Lignières et al. 2006;Reese et al. 2006Reese et al. , 2009Reese et al. , 2013Ouazzani et al. 2012) that treats the eigenvalue problem defined by Equations (6)-(9) as fundamentally two-dimensional.
In this approach, the perturbations are expanded in a series of spherical harmonics and their angular derivatives. For example, we expand a given scalar perturbation X as where ζ is a quasi-radial coordinate that matches the oblate surface of the planetary model at ζ = 1, and transitions to spherical radius at ζ = 0 and ζ = 2. We only solve Laplace's equation in the vacuum region ζ ∈ [1, 2], which is important to include in the computational domain because it permits the application of boundary conditions for the gravitational potential on a spherical surface ). Substituting such expansions into the linearized equations and projecting onto an arbitrary spherical harmonic Y m (this amounts to a Galerkin spectral treatment of the partial differential equations in the angular direction), the dependence of the background state on both radius and colatitude prevents the clean separation of variables that occurs for non-rotating stellar models. Instead, we are left with an infinite series of coupled ordinary differential equations in ζ that must be truncated and solved simultaneously. We employ pseudospectral 'collocation' (e.g., Boyd 2001) for these calculations (see Appendix A for further details). The new code we have developed is similar to that of Xu & Lai (2017), but solves the perturbation equations in a form better-suited to our (non-polytropic) models.
Following Dahlen & Tromp (1998), we normalize the eigenmodes of our rigidly rotating models according to Here ξ = (−iω) −1 v is the Lagrangian displacement assuming a rigidly rotating background state, and the inner product is defined as ξ i , ξ j := V ρ 0 ξ * i · ξ j dV . Under this normalization, rotating-frame mode energies are given simply by ε i = 2ω 2 i (e.g., Schenk et al. 2002).

Perturbative treatment of differential rotation
Defining the Lagrangian displacement more generally as ξ = ∆r, we note that Eulerian and Lagrangian perturbations X and ∆X are related by X = ∆X −ξ·∇X. Taking the Lagrangian variation of Equation 2 then produces (e.g., Lynden-Bell & Ostriker 1967) whereD t = ∂ t +u 0 ·∇, with u 0 a velocity field satisfying the mechanical equilibrium Hereρ 0 ,G = ∇P 0 /ρ 0 andΦ 0 denote a new density, effective gravity, and gravitational field associated with steady flow in the rotating frame; via Equation 13, a non-zero velocity field u 0 will modify the equilibrium pressure, density and gravity field from those of a rigidly rotating fluid body. For our application to Saturn, we consider a velocity field with the form u 0 = Ω D × r, where Ω D = Ω D (r, θ)ẑ is axially symmetric. Then u 0 · ∇u 0 = −RΩ 2 DR , and 2Ω S × u 0 = −RΩ S Ω DR . Figure 2 demonstrates that Ω D is small-amplitude, but varies rapidly over relatively small spatial scales. We therefore ignore the effect of differential rotation on the background pressure, density and gravitational field in Equation 13, but retain terms involving u 0 and its gradient in Equation 12.
With such a velocity field and the modal dependence ξ ∝ exp[i(mφ−ωt)] (for which prograde modes have positive azimuthal wavenumber m and frequency ω), Equation 12 can be written in the operator form where T is the identity (represented with this notation because ξ, T [ξ] relates to kinetic energy), and Given the eigenmodes {ξ j } of a rigidly rotating planetary model, we make the ansatz that normal modes of the differentially rotating model can be expressed as a series expansion ξ = j c j ξ j , with coefficients c j to be determined. Inserting this expansion, we write the momentum equation as (18) where superscripts R and D distinguish operators associated with rigid and differential rotation: is not necessarily zero. Taking the inner product with an arbitrary ξ * i then produces a matrix equation A · c = 0 for coefficient vectors c. Writing ξ i , L[ξ j ] := L ij for a given operator L, A has matrix elements Note that V ij cancels because it does not involve the (in general unequal) eigenfrequencies ω j and ω of the rigidly and differentially rotating models, and because we have ignored modification of the background pressure, density and gravitational field by the differential rotation.
If the effects of differential rotation are small enough that the off-diagonal elements of A are negligible, then the equation A · c = j A ij c j = 0 reduces to A ii c i = 0, from which predicted frequency shifts ∆ω i = ω − ω i can be calculated for a given ξ i as Such shifts only remain accurate in the absence of strong mode mixing. In general, we write A = ω 2 B + ωC + D and solve the quadratic generalized eigenvalue problem Once calculated, the coefficient vectors c can be used to reconstruct approximations to the eigenfunctions of the differentially rotating system from the eigenmodes of the rigidly rotating system. We stress that this treatment of differential rotation is only approximate, not least because it ignores baroclinic modification of the equilibrium state. Our perturbative treatment also relies on the assumption that the oscillations of the differentially rotating planet can be represented by a truncated expansion in the modes of a rigidly rotating model. However, given that Saturn's zonal winds constitute such a small perturbation to the planet's bulk rotation, in this preliminary work we operate under the assumption that they only facilitate the interaction of different oscillations, without drastically altering the rigidly rotating mode spectrum.

RESULTS WITH RIGID ROTATION
In this section we provide a qualitative description of the oscillation modes calculated non-perturbatively for our rigidly rotating models of Saturn, before moving on to a quantitative discussion of their modification by differential rotation in Section 5.

Mode characterization
We focus on fundamental modes, and gravito-inertial modes with super-inertial frequencies |ω| > 2Ω S . Acoustic "p modes" may have been detected in Saturn's gravity field (Markham et al. 2020), but possess frequencies large enough that the associated resonances lie interior to Saturn's rings. On the other hand, although the tidal excitation and damping of inertial waves with subinertial frequencies |ω| < 2Ω S may be relevant to the rapid migration of Saturn's moons (Fuller et al. 2016;Lainey et al. 2020), such low-frequency gravito-inertial modes would produce resonances exterior to the C Ring.
The color-plots in Figure 3 show the meridional (r, θ) structure of the gravitational perturbation associated with a selection of oscillation modes calculated for our model with stable stratification extending to r stab = 0.70R S at the equator. The left-most color-plot illustrates the eigenfunction of the m = 3 sectoral f mode. In our models, the f modes with 4 reside in both the convective envelope and the stably stratified interior, where they take on a gravito-inertial character. Their wavefunctions strongly overlap with those of the lowestorder g modes (independent of frequency degeneracy or rotational mixing), to the point of eliminating any meaningful physical distinction; we simply identify the oscillations with the largest surface gravitational perturbations under Equation 11 as the f modes. The f modes' confinement to the envelope becomes stricter at higher spherical harmonic degrees, resolving their distinction from the g modes and making the oscillations more sensitive to Saturn's zonal winds (see Section 5).
Although the stably stratified interiors of our planetary models are simple, Saturn's rapid rotation complicates the g mode spectrum. As found by Ballot et al. (2010Ballot et al. ( , 2012, in the super-inertial frequency range of interest our calculations produce both conventional g modes ( Figure 3, middle left and middle right), and more exotic rosette modes (Figure 3, right; Figure 4; Figure 5). The former oscillations retain a well-defined number of nodes in the quasi-radial and polar directions, and can be clearly identified with non-rotating counterparts characterized by a single spherical harmonic (thoughout, we use ∼ to denote this correspondence). For these modes, rotation acts primarily to concentrate the oscillations toward the equator (Townsend 2003).

Sectoral f mode
Low-g mode High-g mode Rosette mode . Meridional slices (cross-sections along the rotational axis) illustrating the gravitational perturbations for typical oscillations with azimuthal wavenumber m = 3, for which Φ (r, θ, φ = 0) = −Φ (r, θ, φ = π). These modes were calculated for our model with r stab = 0.70Rs, but generally appear in all of our calculations (with slightly different frequencies). Many f modes and g modes retain eigenfunctions dominated by a single spherical harmonic (dominant = 3, 4, 21 for the left-hand panels), but rotation also produces rosette modes with more exotic angular structures (the right-most panel shows an example).   Figure 4. The heatmaps show, at each spherical harmonic degree (y-axis) and quasi-radial coordinate ζ (x-axis), the (arbitrarily normalized) amplitude of the spectral coefficients ψ (ζ) in an expansion of the modes' gravitational perturbations in spherical harmonics (see Appendix A for variable definitions; for a given mode, all the perturbed fluid variables share similar spectral compositions). The panels demonstrate that the eigenfunctions of these oscillations involve many spherical harmonic degrees. Note that the equatorial symmetry of purely rigid rotation enforces strictly odd (even) equatorial parity for the modes with spectral compositions shown in the outer (inner) panels, so the y-axes include = 4, 6, 8, ... ( = 3, 5, 7, ...).
In contrast, rosette modes take on an altogether different angular structure, characterized by rosette patterns traced in the distributions of kinetic energy (Figure 4 shows examples for a variety of m = 3 rosettes) that inherently involve a coupling across many spherical harmonics (demonstrated by Figure 5). Non-perturbative, or at least degenerate perturbative treatments of the partial differential equations are therefore required to calculate these oscillations. Rosette modes can be identified with the rotational mixing of sequences of g modes that in the non-rotating regime have nearly degenerate frequencies, successive spherical harmonic degrees separated by 2, and different radial quantum numbers n (see Fig. 1 in Takata & Saio 2013). They also appear as stable periodic orbits in ray-dynamic calculations (Prat et al. 2016(Prat et al. , 2018. While the appearance and properties of families of rosette modes depend non-trivially on near-degeneracies in the non-rotating g mode spectrum, asymptotic analyses do predict regularities (Takata 2014), some of which we observe in our calculations. For example, the m = 3 modes with kinetic energies illustrated by the two middle panels of Figure 4 always appear at either end of a sequence with rosette patterns that transition from vertical to horizontal alignment as (closely spaced) frequencies increase. This corresponds to an increase in the real parameter "q" considered by Takata (2014) in their JWKB analysis (see their Fig. 2).
For the purposes of our investigation the rosette modes play an effectively similar role to more conventional high-degree g modes (indeed, the distinction becomes blurry for prograde, non-axisymmetric modes; Saio & Takata 2014). We note, though, that those with vertically rather than horizontally aligned rosette patterns generally couple more strongly with the sectoral f modes, likely due to spectral compositions involving lower-degree spherical harmonics (compare the two middle panels in Figure 5). . The orange points denote relatively low-degree, low-radial order ( , n 12), conventional g modes like that shown in Figure 3  , the sectoral f mode (blue point; see also Figure 3, left) has the largest surface Φ , closely followed by g modes with eigenfunctions dominated by relatively low spherical harmonic degrees and radial orders n (orange points; Figure 3, middle). Rosette modes (green points; Figure 3, right and Figure 4) have smaller surface perturbations, but manifest as sequences with densely spaced frequencies (in this case close to that of the sectoral f mode).

Surface gravitational perturbations
points) and other high-oscillations (black points) carry the majority of their kinetic energy in the stably stratified trapping cavity (i.e., they possess higher mode inertias). Their surface gravitational perturbations are consequently smaller relative to mode energies. Fuller (2014) showed that avoided crossings between f modes and g modes with nearly equal frequencies can enhance the surface amplitudes of the latter. This is particularly true for the conventional g modes with relatively low and n, which in our calculations already possess large surface Φ (our models have much wider trapping cavities than those of Fuller 2014, so the g modes face a smaller distance over which to evanesce in the envelope). However, as demonstrated by the orange points in Figure 6, these low-degree g modes are sparsely distributed in frequency. Additionally, the strength of frequency "repulsion" between the sectoral f mode and low-degree g modes near avoided crossings makes a close frequency splitting between an f mode and multiple lowdegree g modes of the same equatorial parity unlikely.
Rotational mixing with high-g modes presents a possible alternative explanation for the observed triplet of m = 3 density waves, since (as also demonstrated by the green and black points in Figure 6) the high-degree g modes are very densely packed in frequency space; the rosette modes, in particular, appear in sequences with frequencies separated by a few degrees per day (comparable to the finest frequency splitting observed for m = 3). For models similar to our fiducial case with r stab = 0.70R S , the two rosette sequences shown by green points in Figure 6 also appear suggestively close in frequency to the sectoral m = 3 f mode (replacing the conventional high-g modes that mix together to form them). The significance of this proximity should not be over-interpreted, however; rosette modes appear elsewhere in the dense spectrum of high-modes shown in Figure 6, and the two highlighted sequences shift in frequency with changes to the background stratification of the model. Further, as noted by Fuller (2014) rotational mixing between the f modes and high-degree g modes is weak. Although Figure 6 indicates an enhancement for oscillations with frequencies close to ω ∼ 1.35Ω dyn that is common to our models, the surface gravitational perturbations of high-, m = 3 g modes and rosette modes generally fall several orders of magnitude below that of the sectoral f mode for all our calculations with purely rigid rotation (at least under the normalization of Equation 11).

RESULTS WITH DIFFERENTIAL ROTATION
Saturn's zonal winds are strong by terrestrial standards ( 300m/s), but modest in comparison to both the planet's bulk rotation rate and its dynamical frequency (Ω D 0.03Ω S ∼ 0.01Ω dyn ). As a result, we find that the associated differential rotation modifies Saturn's oscillations subtly.

Frequency shifts
The zonal winds directly affect f modes more strongly than g modes. Figure 7 plots the per cent changes in frequency produced by a zonal wind with a decay depth of d = 0.125R S 7.5 × 10 3 km, for f modes with even (left) and odd (right) equatorial parity calculated from our fiducial model with r stab = 0.70R s . The points and plus signs respectively denote frequency shifts calculated by solving Equation 25, and from the approximate frequency shift formula given by Equation 24. Shifts in frequency generally increase with azimuthal wavenumber, a trend explained by the increasingly strict confinement of eigenfunctions to the envelope with increases in the dominant spherical harmonic degree. The lefthand panel in Figure 8 shows the equatorial profiles of Φ (r, θ = π/2) for the sectoral f modes, illustrating this stricter confinement for larger ∼ m. The faint gray line plots the decay profile enforced for the winds (scaled to match the amplitude of the eigenfunctions), demonstrat-ing the domain of influence for the differential rotation profile. Figure 8 (right) shows surface profiles of Φ (r = r s , θ) for m = 15 f modes, demonstrating an additional geometric consideration: Saturn's rapid rotation enhances equatorial confinement, especially for the sectoral f modes. As a result, the high-∼ m sectoral f modes' localization coincides almost exactly with Saturn's rapid equatorial jet [the scaled surface profile of Ω D (r = r s ) is plotted in gray], while the eigenfunctions of the higher-degree oscillations extend to latitudes dominated by alternating regions of under-rotation and super-rotation. These differences in equatorial confinement provide a natural explanation for the monotonic relationship between frequency shift and m for the sectoral oscillations, vis-á-vis the non-monotonic variation for the other f modes.
These geometrical considerations extend to the lowdegree sectoral f modes, even though they are less affected by the zonal winds. Figure 9 plots the per cent frequency changes for the m = 3 sectoral f mode calculated for all of the models shown in Figure 1, for decay depths increasing from zero to an exaggerated d = 30 × 10 3 km. The calculations shown in this plot alone were performed with the antisymmetric component of the wind suppressed, in order to avoid discontinuities deep in the interior (see Section 2.2). The nonmonotonic variation in frequencies with increasing depth can be directly associated with alternating regions of under-rotation and super-rotation at higher latitudes extending deeply enough to overlap with the oscillations' eigenfunctions.
In particular, the relatively large decrease in the frequency shifts plotted in Figure 9 (between d ∼ 5×10 3 km and d ∼ 15 × 10 3 km for all three models) is due to an increased overlap of the m = 3 f mode eigenfunctions with under-rotation at |µ| ∼ 0.5 − 0.6 on the surface (see Figure 2, left) as it extends to meet the equator at r(µ = 0) ∼ 0.8R S . This reduction of the frequency enhancement by zonal winds is relevant given the inference of a sub-corotating layer based on measurements of Saturn's gravitational moments (Iess et al. 2019), but does not significantly alter the frequency shifts of the highdegree sectoral modes shown in Figure 7 (since their eigenfunctions are concentrated toward the surface).

Mode mixing
The zonal wind profiles used in our calculations have little-to-no direct impact on most of the g modes of our models. For all but the lowest degrees and radial orders n, the frequencies of isolated g modes remain largely unaltered except for very deep decay depths . Left: equatorial profiles of Φ for the sectoral ( − m ∼ 0) f modes, calculated from our r stab = 0.70Rs model with/without the effects of differential rotation (dashed/solid lines), and superimposed on the (scaled) decay profile assumed for the zonal winds (thick gray line). Right: surface profiles of Φ for m = 15 f modes with predominantly even equatorial parity (even − m), plotted over the (scaled) surface differential rotation profile ΩD. The left-hand panel demonstrates increasing confinement toward the envelope, likely responsible for larger frequency shifts at larger ∼ m. Meanwhile, the right-hand panel illustrates the overlap between the confinement of the sectoral modes near µ = 0 and Saturn's equatorial jet, which explains the sectoral modes' larger increases in frequency (see Figure 7, left) . d 20 × 10 3 km. This is not surprising, since all but the lowest order oscillations are confined strictly to the stably stratified cavity in the deep interior, where we impose negligible differential rotation. However, highdegree g modes and rosette modes can be affected indirectly through rotational mixing with f modes whose frequencies are altered by the zonal winds. In the following subsections we describe a few examples of such interactions, before discussing observational implications.

Mixing with high-degree g modes and rosettes
Figure 10 illustrates avoided crossings encountered by the m = 3 f mode as the zonal wind decay depth increases for our model with a Brunt-Väisälä profile extending to r stab = 0.65R S at the equator. The left-hand panels plot the frequencies (top) and maximum surface gravitational perturbations (bottom) for the m = 3 sectoral f mode (black circles) and nearby g modes (blue, orange and green plus signs). Since the high-degree g modes and rosette modes are distributed densely in frequency-space, we find that for any given model, at least one decay depth 5 × 10 3 km d 10 4 km produces an avoided crossing with the sectoral m = 3 f mode. During these avoided crossings, the f mode and g mode frequencies nearly overlap, and the oscillations trade character. Close to an avoided crossing the f mode imparts a larger gravitational perturbation to the g mode. This is illustrated by the color-plots in the right-hand panels of Figure 10 . Per cent frequency changes for the sectoral ∼ m = 3 f mode as a function of the wind decay depth, calculated (with equatorially symmetric wind profiles) for the three models considered in this paper. Triangles denote predictions from Equation 24, which can lose accuracy when modes are strongly mixed (this occurs for r stab = 0.65RS near d ∼ 5 × 10 3 km). The non-monotonic variation occurs as Saturn's high-latitude winds descend deeply enough to interact with the oscillation.
both with (right) and without (left) differential rotation, for the relatively high-degree ( ∼ 17) g mode with frequency and surface Φ indicated by the red square. Note that unlike Figure 5, the spectral decompositions shown in Figure 10 include even-components, which remain negligible for this m = 3 mode with even equatorial parity (even − m). Both the meridional slices and the spectral decompositions show an enhancement of the = 3 component of the mode's eigenfunction.
Such avoided crossings provide an attractive explanation for observed m = 3 density waves with frequencies split by less than one per cent. A major caveat is that, as shown in Figure 10, close matches in frequencies (with differences 0.01 − 0.1 per cent) are still required to endow the high-g modes and most rosette modes with surface gravitational perturbations comparable to that of the m = 3 f mode (under the normalization of Equation 11, at least). Consequently, while our calculations for rigidly rotating models typically possess a spectrum of high-g modes and rosette modes that is dense enough near the f mode's frequency to produce at least one avoided crossing at a zonal wind decay depth 10 4 km, explaining the observed m = 3 triplet likely requires at least one other flavor of mode mixing.
For the model with r stab = 0.65R S considered in Figure 10, the m = 3 f mode additionally falls close in frequency to a lower-order g mode, which in the nonrotating regime can be identified with spherical har-monic degree = 5 and radial order n = 3. Because this mode's solution comprises lower-spherical harmonics, the oscillation does not need to be as close in frequency to the f mode to gain a comparable surface amplitude. In fact, Figure 10 (top left) demonstrates clear frequency repulsion between this lower-degree g mode and the sectoral f mode, which prohibits the two oscillations from coming closer than ∼ 0.1 per cent in frequency separation. Figure 9 also illustrates this repulsion, which causes a deviation between frequency shifts calculated in full (blue circles) and from the approximate Equation 24 (blue triangles) near d ∼ 5 × 10 3 km.
The meridional slices in Figure 11 illustrate the eigenfunctions for the three modes from Figure 10 with the largest surface Φ at d = 7.5 × 10 3 km. This decay depth happens to agree with current estimates (Galanti & Kaspi 2021), but in our calculations d is primarily significant for determining the modes' proximity in frequency space; the enhancement of mixing between f modes and g modes by differential rotation does not seem to depend heavily on decay depths between 5 × 10 3 − 10 4 km.
As we discuss in Section 6.2, this type of interaction between i) the m = 3 sectoral f mode, ii) a relatively low-degree g mode mode, and iii) a higher-degree g mode or rosette mode could conceivably explain the observed m = 3 triplet of density waves. However, in general we find that engineering this type of three-mode interaction still requires fine-tuning of models (in order to place a lower-g mode close in frequency to the f mode), which in this case worsens agreement with the observed m = 2 density waves (see Section 6). Figure 12 illustrates another very interesting possible explanation for the three finely split m = 3 density waves. The figure shows similar calculations to those in Figure 10, but for our fiducial model with a profile for the buoyancy frequency that extends to r stab = 0.70R s . For this model (and many others like it), the m = 3 f mode falls very close in frequency to the ∼ 4, m = 3, n = 2 g mode (see the orange plus sign with the largest surface Φ near ω ∼ 1.35Ω dyn in Figure 6). We find that the antisymmetric part of a differential rotation profile like that shown in Figure 2 provides a weak coupling between these two oscillations with differing equatorial parity. The result is an avoided crossing between modes with asymmetric equatorial structure, their eigenfunctions composed of both even and odd-degree spherical harmonics.

Equatorial parity-mixing
This coupling is weaker than that between f modes and g modes with the same equatorial parity, and unlikely to cause significant mixing between a sectoral f  (top) and equatorial surface gravitational perturbation (bottom) with zonal wind decay depth for the m = 3 f mode (black points) and a few nearby g modes (plus signs) calculated from our model with r stab = 0.65RS. The zonal winds have little direct effect on the frequencies of the high-degree g modes, but they gain enhanced surface Φ from avoided crossings encountered as the f mode's frequency changes with decay depth (e.g., the red square at d = 7.5 × 10 3 km). Right: meridional slices (top) and spectral decompositions (bottom; note that unlike the decompositions shown in Figure 5, these include both even and odd ) showing the gravitational perturbation of the ∼ 17 mode denoted by the green plus signs and red square in the left-hand panels, both in the absence of differential rotation (left) and during the avoided crossing at d = 7.5 × 10 3 km (right). These color-plots illustrate the enhancement of the = 3 component of the ∼ 17 mode due to mixing with the f mode. mode and an odd-parity g mode of high-degree. Our model with r stab = 0.72R S does exhibit parity-mixing between the sectoral m = 3 f mode and an odd-parity, ∼ 18 g mode (along with a simultaneous coupling with an even-parity, ∼ 11 g mode). However, the eigenfunction of this higher-degree, odd-parity g mode only gains a significant = 3 component when its frequency differs from that of the f mode by 10 −4 per cent. Aside from requiring a significant amount of fine-tuning, such a small frequency separation so enhances the concentration of kinetic energy in the interior for the f mode (i.e., increases its inertia) that it no longer has the largest surface Φ under the normalization of Equation 11. Low-degree, odd-parity g modes, on the other hand, do not suffer the same frequency repulsion from the sectoral f modes as sectoral g modes. At the same time, they are not as disadvantaged by a high concentration of kinetic energy in the interior as high-degree g modes and rosette modes. Therefore, while substantial mixing between a sectoral f mode and a high-degree, odd-parity g mode is unlikely, a three-mode interaction involving i) the m = 3 f mode, ii) an odd-parity, low-degree g mode, and iii) an even-parity g mode may provide a viable explanation for the finely split triplet. We discuss the observational implications of such an interaction, and how it might be confirmed, in Section 6.2.  Figure 10, but for our fiducial model with r stab = 0.70RS, for which the frequency of the m = 3, ∼ 4, n = 2 g mode (with odd equatorial parity) falls close to that of the sectoral m = 3 f mode. The antisymmetric part of Saturn's zonal winds is weak, but sufficient to mix this odd-parity g mode with the even-parity f mode, producing two oscillations with eigenfunctions involving both even and odd harmonics (as shown for the g mode in the right-hand spectral decomposition). Asymmetric oscillations like that shown in the top-right panel may have a unique observational signature (see Section 6.3), namely the simultaneous excitation of both density and bending waves at different locations in the C Ring.

DISCUSSION
In this section we discuss the implications of the results described in Section 4 and Section 5, and their relation to the waves observed in Saturn's C Ring.

Comparison with observed frequencies
The frequencies of the f modes in our calculations agree broadly with both the perturbative calculations of Mankovich & Fuller (2021), and the majority of the waves observed in Saturn's C Ring. Figure 13 plots azimuthal wavenumbers m against the resonant radii at which f modes (blue points) and low-order g modes (orange/green points) with predominantly even (left) and odd (right) equatorial parity would be expected to excite density waves and bending waves, respectively. We plot the observed waves as black diamonds. Filled circles denote oscillations calculated with our purely rigidly rotating model with positive Brunt-Väisälä frequency extending to r stab = 0.70R s , while plus signs show changes in the resonant radii due to frequency shifts from differential rotation (again assuming a radial decay depth of d 7.5 × 10 3 km). We scale the sizes of all the points by 1 + log 10 (Φ 3 /Φ s ), where Φ s and Φ 3 are the maximum surface gravitational perturbations for the mode in question and the m = 3 sectoral f mode, respectively. Figure 13 indicates nearly universal agreement with the observations for m ≥ 5, which is improved for the sectoral f modes (right-most blue track in the left-hand panel) by increases in mode frequency from Saturn's zonal winds. At m = 4, an avoided crossing with the lowest-order (n = 1) sectoral g mode makes the choice of which mode to call the f mode arbitrary. As discussed in Section 4, we identify the oscillation with the largest surface gravitational perturbation as the sectoral f mode. Regardless, in this case neither oscillation quite matches the observed density wave excited at r 81 × 10 3 km.
This discrepancy is generic to our models. In the absence of additional observations of more m = 4 density waves, it points toward needed refinement of our picture of Saturn's deep interior, since the low-m oscillations are the most affected by interactions with g modes trapped in the stably stratified cavity. Minor discrepancies between the observed m = 2 density waves and our calculations for the m = 2 f mode and ∼ m = 2, n = 2 g mode further emphasize this need, although Mankovich & Fuller (2021) found that a sharp cut-off in buoyancy frequency can help by lessening the eigenfunction overlap between these two modes, allowing their frequencies to come closer. We have focused on models with smooth variation in composition for numerical convenience, but models with discontinuous or rapid variation may provide a more appropriate description.

Optical depth variations and detectability
To identify gravito-inertial modes potentially involved in the excitation of density waves with finely split frequencies, we estimate optical depth variations induced in the C Ring by the oscillations' gravitational perturbations, following the approach of Fuller et al. (2014)  and Fuller (2014). Appendix B reviews the method of these calculations.
6.2.1. Comparison between rigid and differential rotation Figure 14 plots optical depth variations as a function of pattern speed Ω p = σ/m (where σ = ω + mΩ S is the inertial-frame frequency) for m = 2 (left) and m = 3 (right) oscillations calculated for our fiducial (r stab = 0.70R s ) model, both with (orange) and without (black) the effects of Saturn's zonal winds assuming a radial decay depth of d 7.5 × 10 3 km. Points falling in the un-shaded region predict visible density waves, while those falling in the gray-shaded regions would be expected to produce optical depth variations too small to be observable, or resonances lying outside the C Ring. The black diamonds again correspond to the observed density waves.
The majority of modes show only minor changes in predicted |δτ | due to differential rotation. On the other hand, mixing between the even-parity sectoral f modes and odd-parity g modes produces a host of predicted optical depth values |δτ | 10 −5 with no identifiable counterparts from the calculation with rigid rotation, since in the absence of an antisymmetric wind, oddparity modes calculated with purely rigid rotation have Φ = 0 at the equator (and therefore cannot drive density waves). Additionally, the asymmetric avoided crossing shown in Figure 12 produces two finely split, ostensibly visible waves close to the observed m = 3 triplet with Ω p 1736.7, 1735.0, and 1730.3deg d −1 (Hedman & Nicholson 2013).
Very high-, even-parity modes with frequencies close to the m = 3 f mode (such as the rosette mode shown in the middle-right panels of Figure 4 and Figure 5) do show |δτ | enhancements in some of our calculations with differential rotation (depending on the model and wind decay depth). However, we find that modes with solutions involving m + 30 require unreasonably small frequency separations from the sectoral f mode to produce observable optical depth variations. We note that more detailed treatments of differential rotation involving a fully self-consistent modification of the background pressure and density by the zonal winds may produce stronger coupling between the f modes and high-g modes. Also, our predictions of optical depth variations are sensitive to an assumption of energy equipartition between modes; any internal process that preferentially excites g modes and/or rosette modes to larger energies than f modes would be misrepresented by this assumption, imbuing the former with larger predicted |δτ |.
Barring chance over-excitation by stochastic processes, such preferential excitation of g modes in the deep interior would likely require a different mechanism than those already considered for Saturn (Markham & Stevenson 2018;Wu & Lithwick 2019), which should operate most efficiently near the surface (where mode excitation is also most efficient in the Sun; Goldreich et al. 1994). The possibility that double-diffusive convection (e.g., Moll et al. 2017) in Saturn's interior might produce alternating layers of convection and stable stratification opens the door to potentially exotic excitation . Predicted optical depth variations for m = 2 (left) and m = 3 (right) modes calculated from our model with a stably stratified region extending to r stab = 0.70Rs, and plotted against pattern speeds Ωp = σ/m = ω/m + ΩS (here σ is the inertialframe frequency). The black diamonds show the observed density waves, and the orange (black) points show calculations with (without) the inclusion of differential rotation from Saturn's zonal winds (assuming a decay depth of d = 0.125RS 7.5×10 3 km). Differential rotation only appreciably shifts the frequencies of the f modes (by ∼ 0.1 per cent), but can alter the optical depth perturbations predicted for g modes. . Meridional distributions of kinetic energy for a sequence of m = 2 rosette modes with closely split frequencies that produce optical depth variations |δτ | ∼ 10 −2 near Ωp ∼ 1750deg d −1 for our model with r stab = 0.70Rs. This sequence produces nearly detectable optical depth variations in most of our models, and could be associated with the observed m = 2 density waves with finely split frequencies.
mechanisms, although the formation and persistence of high-degree g modes and rosette modes with relatively short wavelengths may be more dubious in a segmented cavity (Belyaev et al. 2015;Pontin et al. 2020). The optical depth predictions in the left-hand panel of Figure 14 show rough agreement with the observed m = 2 density waves: the m = 2 sectoral f mode produces a potentially detectable |δτ | close to the observed density wave with Ω p 1860.8deg d −1 (Hedman & Nicholson 2013), while the sectoral n = 1 g mode matches well with the observed density wave with Ω p 2169.3deg d −1 ). The ∼ m = n = 2 g mode produces a potentially detectable wave with a pattern speed close to, but somewhat less than the observed pair of waves at Ω p 1779.5deg d −1 (Hedman & Nicholson 2013) and Ω p 1769.2deg d −1 (French et al. 2016). As already mentioned, profiles of the Brunt-Väisälä frequency with a discontinuous or sharply varying cut-off in the envelope can reduce the g mode frequency spacing, and therefore shift the n = 2 g mode closer to observations. We also note that our calculations with most models produce a sequence of relatively low-degree rosette modes in this frequency range. For our fiducial model, these rosette modes manifest in the left-hand panel of Figure 14 as sequences of both black and orange points with optical variations |δτ | ∼ 10 −2 that monotonically increase with pattern speeds Ω p 1750deg d −1 . The color-plots in Figure 15 show meridional distributions of kinetic energy for example rosette modes in this se-  Figure 14 but for all three of our models with r stab = 0.65RS (blue), 0.70RS (orange) and 0.72RS (green), all including differential rotation with a decay depth d 7.5 × 10 3 km.
quence. While the optical depth variations calculated for these modes lie just below the detectable regime, their consistent appearance (in many of our calculations with different models) in this frequency range is suggestive. With sufficient mixing with the the n = 2 sectoral g mode, or a departure from our assumption of energy equipartition between modes (i.e., with more energy), such rosettes could be involved in exciting the finely split m = 2 density waves.

Comparison between models
The panels in Figure 16 show |δτ | predictions for the three models with r stab = 0.65R S , 0.70R S and 0.72R S considered in this paper, calculated with decay depth d 7.5 × 10 3 km. As illustrated by the blue points associated with our model with r stab = 0.65R S , the threemode interaction shown in Figure 10 produces a potentially detectable triplet of density waves that lie close in frequency to the observations. However, the frequency splitting of less than ∼ 0.1% between all three modes in this calculation is in fact too close, since the observed pattern speeds of (Ω p 1736.7, 1735.0, 1730.3deg d −1 ; Hedman & Nicholson 2013) are split by ∼ 0.1 − 0.35%.
Additionally, while pairs of detectable waves associated with the mixing of the m = 3 f mode with a single high-degree g mode are not difficult to produce in our calculations (due to denser spacing in frequency space for high-modes), engineering the overlap of such a pair with a third, relatively low-(∼ 5 in this case) g mode requires significant fine-tuning. For the model with r stab = 0.65R S , shifting the overall g mode spectrum to higher frequencies also shifts the m = 2, sectoral g mode with n = 1 away from the observed pattern speed (compare the points near ∼ 2200deg d −1 in Figure 16, left). On the other hand, our model with r stab = 0.72R S significantly under-predicts the frequency of both the n = 1 and n = 2 sectoral g modes, the frequency for the latter falling into the sub-inertial range (which we do not consider here). Figure 16 (right) illustrates another characteristic shared by all of our calculations using the models of Mankovich & Fuller (2021). Invariably, they predict that the n = 1, m = 3 sectoral g mode with pattern speed between Ω p ∼ 1800deg d −1 and 1950deg d −1 (depending on the model) should excite a potentially detectable density wave. This oscillation even possesses the largest predicted |δτ | in our model with r stab = 0.72R S , due to an enhancement of the original f mode's mode inertia via mixing with two g modes with ∼ 11, 18 and nearly degenerate frequencies (see Section 5.2.2). Strong frequency repulsion from the sectoral f mode likely prohibits the n = 1 sectoral g mode from playing a role in the observed m = 3 triplet. The fact that additional m = 3 waves with more widely separated frequencies have not been observed may suggest that our model for Saturn's deep interior still needs refinement, or that energy equipartition may not be good assumption.

Simultaneous density and bending wave excitation
Lastly, we discuss the implications of our finding that even small antisymmetric components of a differential rotation profile can couple oscillations with even and odd equatorial parity. Parity-mixing between a sectoral f mode and an odd-parity, low-degree g mode (see Figure 12) is less likely than avoided crossings involving high-degree g modes of the same parity (see Figure 10); both require a close frequency degeneracy ( 0.1%) with the sectoral f mode, but low-g modes are much more sparsely spaced in frequency. Notably, though, we do find that the sectoral m = 3 f mode and the ∼ 4, m = 3, n = 2 g mode have similar frequencies in many of the best-fitting models of Mankovich & Fuller (2021).
More importantly, such a coupling would have a very unique observational signature: modes with asymmetric equatorial parity could conceivably excite both density and bending waves in Saturn's C Ring, at slightly different radii due to differences between the horizontal and vertical epicyclic frequencies associated with orbital motion. The interaction of the m = 3 f mode with the ∼ 4, m = 3, n = 2 g mode in our calculations could, for example, be confirmed by the additional detection of one or more m = 3 bending waves at the resonant radius r V 83 × 10 3 km (in addition to the observed density wave excitation at r L 82.2 × 10 3 km). If such paritymixing were confirmed, the sparseness of the low-degree g mode frequency spectrum would serve as a powerful asset for seismic inference, since it would significantly reduce the number of mode interactions potentially capable of producing the observed fine splitting.

CONCLUSIONS
We have investigated the effects of Saturn's rapid and differential rotation on the planet's normal mode oscillations, focusing on the mixing of fundamental modes (f modes) and gravito-inertial modes (g modes) likely responsible for exciting density waves with nearly degenerate frequencies in the C Ring (Hedman & Nicholson 2013;Fuller 2014). Using a combination of nonperturbative and perturbative methods to account for Saturn's rapid bulk rotation and deep zonal winds (resp.), we have computed oscillation modes for interior models featuring wide regions of stable stratification constrained by a joint fit of Saturn's low-order seismology and gravity field (Mankovich & Fuller 2021).
The stably stratified interiors in our models produce a rich spectrum of gravito-inertial oscillations (see Figure 6), some adhering to the conventional understanding of g modes in non-rotating planets and stars, and others resembling so-called "rosette modes" (Figure 4; Figure 5; Figure 15) encountered only with higher-order treatments of rapid rotation (Ballot et al. 2012;Takata & Saio 2013).
The high-degree g modes and rosette modes appear densely spaced in frequency, making their interaction with the low-m, sectoral f modes commonplace (Figure 10; Figure 11). We have demonstrated that such interactions can endow the otherwise-undetectable highdegree g modes with large external gravitational perturbations (relative to mode energies) ostensibly capable of exciting visible density waves with frequencies nearly degenerate to those driven by the f modes (Figure 14; Figure 16). The avoided crossings are narrow in frequency width, however, and we infer that additional mixing with at least one lower-degree g mode may be required to reproduce an observed triplet of m = 3 density waves.
With our perturbative treatment of differential rotation, we find that Saturn's zonal winds impact its f modes subtly but measurably. Saturn's equatorial jet most significantly alters the frequencies of high-degree sectoral modes, owing to their strict confinement to both the envelope and equatorial latitudes (Figure 7; Figure 8; Figure 13). Differential rotation in Saturn's envelope also enhances the avoided crossings between f modes and g modes at low azimuthal wavenumbers m.
Interestingly, we find that including a realistically small antisymmetric component of the differential rotation profile leads to weak coupling between modes with different equatorial parity ( Figure 12). The weakness of this coupling makes parity-mixing a relatively unlikely explanation for observed fine-frequency-splitting between density waves in the C Ring. If confirmed, however, the rotational mixing between a sectoral f mode and an odd-parity g mode would provide a powerful constraint on Saturn's deep interior, since such mixing would likely be limited to a sparse spectrum of lowdegree g modes. Such an interaction could be confirmed by observations of simultaneous density and bending wave excitation at Lindblad and vertical resonances with identical pattern speeds.
where Y m Lastly, the solution for the gravitational potential that vanishes at infinity is ψ out ∝ r −( +1) , implying the boundary condition 1 r ζ dψ l out dζ + + 1 ζ ψ l out = 0, which we apply on the spherical surface ζ = r = 2.

B. OPTICAL DEPTH VARIATIONS
Consider the gravitational potential produced by each mode in an inertial frame external to the planet: where σ = ω + mΩ S is the mode's inertial-frame frequency, and A its amplitude. In an equatorial ring with orbital frequency Ω r and surface density Σ r , the effective forcing potential for a density wave with azimuthal wavenumber m is evaluated at the Lindblad radius r L (Goldreich & Tremaine 1979). This Lindblad radius can be computed by solving the equation (mΩ r − σ) = −κ r , where κ r is the horizontal epicyclic frequency. Following Mankovich et al. (2019), we use a multipole expansion for Ω r and κ r , and solve for r L numerically. The forcing potential Ψ excites spiral density waves that, in the (thoroughly justified) "tight-winding" limit and WKB approximation, induce gravitational perturbations with the form where x = (r − r L )/r L , and D 3(m − 1)Ω 2 r | r L near a Lindblad resonance (Cuzzi et al. 1984). In turn, Φ r can be related to a surface density perturbation δΣ and optical depth variation δτ = κ m δΣ (here κ m is the local mass extinction coefficient) via δΣ = δτ κ m = i 2πGr 1/2 d(r 1/2 Φ r ) dr xΨ 2π such that at a given location x in the ring, |δτ | A 3(m − 1)κ 2 m Ω 2 r |Ψ/A| 2 4π 2 G 3 Σ r r L 1/2 |x|.
(B48) Like Fuller (2014), we use this direct relationship between δτ and mode amplitude to calculate the amplitude |A 3 | of the ∼ |m| = 3 f-mode that would be required to produce the optical depth variation δτ 0.21 observed for the largest amplitude |m| = 3 density wave in Saturn's rings. The assumption of energy equipartition with ω 2 |A| 2 = ω 2 3 |A 3 | 2 then provides amplitudes for the rest of the modes (all of which are normalized ahead of time according to Equation 11). We set |r − r L | = 5km, and adopt the nominal values κ m = 0.02cm 2 g −1 and Σ r = 5gcm −2 . These are roughly appropriate to the location of the m = 3 triplet (Hedman & Nicholson 2013), but we note that conditions vary at the locations of many of the other observed waves.