Possible Evidence of p-modes in Cassini Measurements of Saturn's Gravity Field

Gravity field measurements allow us to probe the interior structure of a planet by measuring its deviation from spherical symmetry. For giant planets, the planet's response to its own rotation breaks its spherical symmetry. The deviation away from spherical symmetry depends on the planet's internal density distribution (Stevenson 2019). Therefore a detailed mapping of a planet's gravity field can corroborate or refute interior models. Saturn's non-spherical gravity field was first inferred from spacecraft tracking data of Pioneer 11 (Hubbard et al. 1980; Null et al. 1981), and was later improved using Voyager data (Campbell & Anderson 1989). The arrival of the Cassini spacecraft in the Saturnian system yielded more accurate determination of the gravity field of the gas giant by first looking at the orbits of its satellites. Now the Grand Finale of the Cassini mission has produced exquisite gravity field data for Saturn, providing the first concrete constraints for Saturn's ring mass, zonal wind depths, and evidence for internal differential rotation by offering gravity field measurements up to J₁₂ (Galanti et al. 2017; Iess et al. 2019). But behind these spectacular new findings lurks a dark side: a small component of Saturn's gravity field that cannot be explained with the canonical static, zonally symmetric gravity field expected of gas giants.

Cassini's radioscience experiment is carried out by measuring the Doppler shift of a microwave signal in a two-way configuration: the signal is sent from a ground station to the spacecraft, which retransmits it back to the station preserving phase coherency. The Doppler shift is, to first order, proportional to the relative velocity of the spacecraft with respect to the station. These measurements are compared with predictions based on dynamical and observation models to obtain data residuals. The data we use in this study are two-way Doppler residuals, converted in a radial velocity time series, obtained by removing the effect of empirical acceleration from the reference solution given by Iess et al. (2019).

This additional and unknown source of gravity can be fit with a variety of models. A static tesseral gravity field is possible, but there is no convincing low-order fit (Iess et al. 2019). A low-order tesseral field does not provide a predictive solution with the available data, and also depends on the assumed rotation rate of Saturn. That is, a given gravity harmonic solution for a subset of flybys will not accurately predict the next flyby and requires additional harmonic terms. The nominal method that was employed for the published gravity harmonic results was an agnostic "empirical acceleration" model that, due to the unknown origin of the source of the additional gravity, included random acceleration vectors that changed on a 10 minute timescale. In this context "random" means that each acceleration vector is allowed to have any direction with an a priori amplitude of ±4 × 10⁻¹⁰ km s⁻¹ (Iess et al. 2019). They could be correlated (non-random) even when the process used to create them allows for randomness. This timescale between changing acceleration vectors was determined empirically as the longest timescale that can successfully reduce range-rate residuals to the noise level.

A time-dependent signal does not necessarily require normal modes. For example, there may be a time-dependent or non-symmetric signature from large-scale convection (Kong & Schuber 2016). Additionally, Saturn's envelope is differentially rotating (Galanti et al. 2017; Chachan & Stevenson 2019). If mass anomalies were embedded at different depths or latitudes, then a spacecraft could encounter measurably different quasi-static tesseral gravity fields during each flyby (Iess et al. 2019). However, differentially rotating tesseral structure in Saturn's gravity field has been shown to produce structures in the rings (El Moutamid et al. 2017), and the magnitude of the potential perturbation inferred from observation is orders of magnitude too small to explain the anomalous signal. Because of Saturn's expected internal differential rotation rate (about 5%) (Chachan & Stevenson 2019; Galanti et al. 2017), it is unlikely that such structure could measurably affect the spacecraft trajectory without showing clear structure from resonances in the rings.

This work will specifically explore the hypothesis that Saturn's residual gravity is a consequence of normal-mode oscillations. It has already been demonstrated that normal modes are capable of eliminating the range-rate residuals to the noise level (Iess et al. 2019). This has been done by computationally optimizing for individual mode amplitudes using a large number of free parameters. One possible solution involves only zonal f-modes. This solution, however, is affected by model assumptions such as maximum modeled spherical degree, whether to permit p-modes or g-modes, whether to permit non-zonal normal modes, etc. These uncertainties occur because, when optimizing with a large number of free parameters, there is a risk of over-fitting the data using a too-complex model. These issues are not important in the context of constraining Saturn's zonal gravity harmonics and ring mass because the uncertainty can simply be absorbed in the error ellipses for these values. However, in this work we revisit the residuals data with a different purpose: to try to extract a preferred normal-mode spectrum that is predictive for further flybys, robust to changes in model assumptions, and as simple as possible to capture the qualitative behavior of the spectrum without over-fitting the data. Bearing this in mind, although we find a statistical preference using our simple model for signals dominated by low-order p-modes, readers should remember that our findings are not conclusive proof of such a spectrum on Saturn.

Our investigation has at least two important applications: first, normal modes are themselves a promising method by which to probe the interior structure of giant planets, and this analysis provides some evidence of their power spectrum. Second, any gravitational signal from normal modes above the noise level contaminates spacecraft tracking data and may be aliased into the static model. As we will see, the behavior of the modes are partially degenerate and the solution is non-unique. However, the solutions are clustered in parameter space and predict a high probability of reproducing the observed unexplained gravity signal. The most successful models indicate the signal is likely to be dominated by p-modes between 500 and 700 μHz (see Figure 4).

In the Section 2, we outline some fundamentals of giant planet seismology, spacecraft tracking, and our forward model. In Section 3 we discuss our data reduction method including a novel data-stacking technique, as well as error sources, and a fitting procedure. In Section 4 we present the results of our analysis, finding a simple two-parameter model that has a high probability of producing a good fit to the spacecraft signal. In Section 5 we discuss the implications of our findings.

In order to accurately model seismic effects on Cassini's gravity signal, we must determine the mode's eigenfrequencies, and the scaling relationship between a mode's displacement amplitude and its effect on Saturn's gravity field. These issues are addressed in Section 2.1. Next we must model how a given gravity potential perturbation affects the spacecraft tracking signal, which is done in Section 2.2. Next we need an agnostic parametric model for the modal energy spectrum, discussed in Section 2.3. Finally we account for the intrinsic stochasticity of the problem; a given mode cannot be modeled deterministically, because we have no way of knowing what the temporal phase of each mode was when Cassini was at periapse. This is partially circumvented with our stacking technique, discussed in Section 3 with further technical information in the Appendix.

2.1. Background

In this paper we approximate Saturn as an adiabatic, spherical, uniformly rotating planet. We neglect rotation to compute the eigenfunctions and potential perturbations, but account for rotation when considering Coriolis force frequency splitting and the rotating gravity potential encountered in an inertial frame. In this case, giant planet oscillations can be decomposed into a discrete set of orthogonal normal modes with quantum numbers (n, l, m); n corresponds to the number of radial nodes in the displacement eigenfunction, l to the spherical harmonic degree, and m = [−l..l] to the azimuthal degree. Each mode has a unique displacement eigenfunction

$$\begin{eqnarray}\begin{array}{rcl}{\xi }_{{nlm}}({\boldsymbol{r}}) & = & \left({\xi }_{r,{nlm}}(r)\hat{{\boldsymbol{r}}}+{\xi }_{h,{nlm}}(r)\hat{\theta }\displaystyle \frac{\partial }{\partial \theta }\right.\\ & & +\,\left.\displaystyle \frac{{\xi }_{h,{nlm}}(r)}{\sin \theta }\hat{\phi }\displaystyle \frac{\partial }{\partial \phi }\right){Y}_{{lm}}(\theta ,\phi )\end{array}\end{eqnarray} \tag{ 1 }$$

and a characteristic eigenfrequency ω_nlm so that the total displacement as a function of time is ${\xi }_{{nlm}}\cos ({\omega }_{{nlm}}t+{\alpha }_{{nlm}})$. Because ω_nlm is not precisely determined, the phase α_nlm cannot be coherently specified between flybys and is assumed random for each mode for each flyby. ξ_r and ξ_h correspond to the radial and horizontal eigenfunctions respectively, which together specify the fluid displacement at any point within the planetary sphere. For our purposes the eigenfunctions were obtained using the GYRE stellar oscillation code (Townsend & Teitler 2013), with an n = 1 nonrotating polytrope model for Saturn's interior. Our goal here is independent of accurate interior modeling; we are interested in the relative gravity signal between modes and their order of magnitude, which is not strongly sensitive to small changes in the interior model. However, using an adiabatic interior model precludes g-modes, so we account for contributions from g-modes separately.

Because Saturn's interior structure is not precisely determined, we performed our full analysis on a variety of interior model assumptions to demonstrate that the results are not sensitive to small errors in modal eigenfrequencies. We tested eigenfrequencies produced by this same polytrope model generated with GYRE (Townsend & Teitler 2013), as well as a sampling computed using a more sophisticated Saturn interior model (Sa8; Gudkova & Zharkov 2006). The nominal model uses the eigenfrequencies from Sa8. In addition, we accounted for mode-splitting due to Coriolis forces (Christensen-Dalsgaard 2014). These split according to

$$\begin{eqnarray}&&\delta {\omega }_{{nlm}}=m{\beta }_{\mathrm{nl}}{\rm{\Omega }},\end{eqnarray} \tag{ 2 }$$

where Ω is Saturn's spin rate and

$$\begin{eqnarray}&&{\beta }_{\mathrm{nl}}\equiv \displaystyle \frac{{\int }_{0}^{R}({\xi }_{r}^{2}+l(l+1){\xi }_{h}^{2}-2{\xi }_{r}{\xi }_{h}-{\xi }_{h}^{2}){r}^{2}\rho {dr}}{{\int }_{0}^{R}({\xi }_{r}^{2}+l(l+1){\xi }_{h}^{2}){r}^{2}\rho {dr}}.\end{eqnarray} \tag{ 3 }$$

The nominal frequencies for this paper are plotted in Figure 1. For computational reasons, we consider a finite subset of modes in our model. In the nominal model we consider f-modes and p-modes up to l = 8, n = 7. We found equivalent results when using instead l = 10, n = 5 as bounds on parameter space. We also specially tested f-modes only up to l = 20. We do not expect g-modes to dominate the signal for at least two reasons: first, because the stable layer where they resonate is so deep, its effect on the gravity field would be very weak unless its amplitude were extremely large. Second, we do not expect its amplitude to be extremely large, because its eigenfunction is evanescent near the surface where mode excitation is expected to be most efficient. Nevertheless, for the sake of completeness we tested g-modes using published eigenfrequencies (Gudkova & Zharkov 2006).

Zoom In Zoom Out Reset image size

After choosing eigenfrequencies, we compute the scaling between displacement eigenfunctions and gravity potential perturbations. The gravity field perturbation associated with displacement eigenfunction ξ can be obtained by integrating over the material sphere and accounting for fluid point displacements according to

$$\begin{eqnarray}&&\delta {\rm{\Phi }}=G\int \displaystyle \frac{\rho ({\boldsymbol{r}}^{\prime} )}{| {\boldsymbol{r}}-({\boldsymbol{r}}^{\prime} +\xi ({\boldsymbol{r}}^{\prime} ))| }{d}^{3}r^{\prime} -\displaystyle \frac{{GM}}{r}.\end{eqnarray} \tag{ 4 }$$

For linear perturbations this is equivalent to the Eulerian density perturbation from continuity $\delta \rho ={\rm{\nabla }}\cdot (\rho \xi )$ so that

$$\begin{eqnarray}&&\delta {\rm{\Phi }}=G\int \displaystyle \frac{\delta \rho ({\boldsymbol{r}}^{\prime} )}{| {\boldsymbol{r}}-{\boldsymbol{r}}^{\prime} | }{d}^{3}r^{\prime} .\end{eqnarray} \tag{ 5 }$$

Decomposing this potential perturbation into

$$\begin{eqnarray}\begin{array}{rcl}\delta {\rm{\Phi }} & = & \displaystyle \frac{{GM}}{r}\displaystyle \sum _{n=0}^{\infty }\left(-\displaystyle \sum _{l=2}^{\infty }\displaystyle \sum _{m=-l}^{l}{\left(\displaystyle \frac{R}{r}\right)}^{l}{P}_{l}^{m}(\cos \theta )\right.\\ & & \times \,\left.[\delta {C}_{{nlm}}\cos m\phi +\delta {S}_{{nlm}}\sin m\phi ]\right)\end{array}\end{eqnarray} \tag{ 6 }$$

one can show that the gravity harmonic coefficient perturbation associated with the normal mode is (Marley & Porco 1993; Stevenson 2019)

$$\begin{eqnarray}&&\delta {C}_{{nlm}}(t)=\displaystyle \frac{1}{{{MR}}^{l}}\displaystyle \frac{2l+1}{4\pi }\cos ({\omega }_{{nlm}}t){\int }_{0}^{R}{r}^{l+2}\delta {\rho }_{\mathrm{nlm},r}(r){dr},\end{eqnarray} \tag{ 7 }$$

where $\delta {\rho }_{{nlm}}({\boldsymbol{r}})\equiv \delta {\rho }_{\mathrm{nlm},r}(r){P}_{l}^{m}(\theta )\cos (m\phi )$ is normalized such that the mode surface displacement is 1 cm at the planet's surface. With appropriate choice of coordinates, $\delta {S}_{{nlm}}\to 0$. This leads to Figure 2 that illustrates why f-modes are a priori favored as sources of gravity perturbations. Higher-order p-modes have nodes in their eigenfunction, leading to destructive interference of the gravitational signature. Therefore in order for p-modes to dominate the signal, they must have more than an order of magnitude larger energy than f-modes (two orders of magnitude larger amplitude).

Zoom In Zoom Out Reset image size

A test particle outside of the planet on a prescribed trajectory r(t) encounters the potential perturbation

$$\begin{eqnarray}\begin{array}{rcl}\delta {\rm{\Phi }}(t) & = & \displaystyle \frac{{GM}}{r}\displaystyle \sum _{{nlm}}{a}_{{nlm}}\delta {C}_{{nlm}}{\left(\displaystyle \frac{R}{r(t)}\right)}^{l+1}{P}_{{lm}}(\cos \theta (t))\\ & & \times \cos [m(\phi (t)-{\rm{\Omega }}t-{\phi }_{0})-{\omega }_{{nlm}}t-{\alpha }_{{nlm}}],\end{array}\end{eqnarray} \tag{ 8 }$$

where a_nlm is the maximum surface displacement of mode (n, l, m) in centimeters, and Ω is Saturn's spin rate. Notice that this equation includes two random variables: ϕ₀ is the initial longitudinal orientation of the modes with respect to our coordinates, and α_nlm is the initial temporal phase of the mode. ϕ₀ is shared between all modes, but is random for each flyby. α_nlm is a random variable for each mode and for each flyby. Although the phase difference between flybys can in principle be determined from the mode's eigenfrequency, in practice this is impossible. Eigenfrequencies depend on Saturn's interior structure and cannot be predicted with perfect precision. Because the time between encounters is much longer than the period of a mode, in practice the phase of each mode must be regarded as randomized for each flyby. This stochasticity introduces a complication for modeling the flyby—we do not know the initial phase and orientation of the modes when each flyby occurred. This issue will be addressed in Section 3.

2.2. Numerical Integration and Model Reproduction

The gravity experiments were conducted with an edge-on geometry from Earth's perspective for maximum signal to noise. The orbits were highly inclined and highly eccentric. The closest approach (C/A) of the spacecraft is about 5% of Saturn's radius from its cloud tops, approaching and receding from the planet very quickly during ingress and egress. Because of this orbital geometry, we only expect a significant signal from the planet within about an hour of C/A. Therefore we use the spacecraft orbital elements two hours before C/A to compute the initial conditions. We then numerically integrate the equation of motion for Saturn plus the potential perturbation associated with normal modes, neglecting oblateness. Including the measured static zonal gravity (Iess et al. 2019) affects the simulated normal mode range-rate signal by less than 1%. We integrate the equation of motion using Mathematica's built-in integrator to generate a three-dimensional velocity time series. We then subtract the Keplerian solution from the numerical solution to isolate the signal from the spectrum of modes we are modeling. Finally we project the three-dimensional velocity vector onto the line-of-sight vector connecting Saturn to the Earth. We verify the accuracy of this method by reproducing the signature from values of static zonal J (Iess et al. 2019) using the method from this paper, finding good agreement. This method is fully general for any potential perturbation, and we will use it to inspect the behavior of normal modes.

We verified empirically the approximate linearity of combining the velocity perturbation from various sources. That is,

$$\begin{eqnarray}&&\delta {\boldsymbol{v}}(t)=\displaystyle \sum _{q}{{\boldsymbol{v}}}_{q}(t)\end{eqnarray} \tag{ 9 }$$

within <0.1% for the perturbation magnitudes in question. Strictly speaking this linearity does not hold absolutely; although gravity potential perturbations are exactly linear, a test particle encountering these perturbations may be perturbed from its trajectory—if this perturbation is sufficiently large the linearity breaks down. But for the small perturbation of interest, this nonlinearity is not important.

It is important to discuss at this point a fundamental ambiguity in probing for normal modes from the spacecraft's perspective. The spacecraft is observing two sources of variation of the gravity signal: the intrinsic geometric variation, and the temporal variation. The geometric variation is the physical shape of the mode, which attenuates with distance and varies with the spacecraft's latitude and longitude relative to Saturn. As the spacecraft approaches and recedes from the planet, traveling from north to south and west to east, even a static gravity perturbation would have a time-dependent signal from the spacecraft's frame of reference. On the other hand, the potential perturbation itself varies with time. The convolution of these effects makes it difficult to have a simple intuition for Cassini's response to each mode. See the Appendix for further discussion.

2.3. Spectral Model

For the spectral model, we aim to be as agnostic as possible. We do not know with certainty by what mechanism seismic activity is excited on Saturn, although meteor impacts (Wu & Lithwick 2019) or exotic meteorological phenomena (Markham & Stevenson 2018) have been suggested. Therefore we only use simple parametric models that scale the energy of each mode as a function of parameters. We tested a variety of scaling relationships, including power-law dependence in frequency, as well as power-law dependence on quantum numbers n, l, and $\tfrac{l-| m| }{l}$. We also tested equipartition. None of these models provided convincing fits to the data.

The nominal model used is a Gaussian frequency-dependent model, although here too we aimed to be as agnostic as possible. The Gaussian is convenient because it has the power to probe for a diverse variety of frequency dependences by only varying two parameters. Using this assumption, the mode energy is a function of its eigenfrequency according to

$$\begin{eqnarray}&&E({\omega }_{{nlm}})=\displaystyle \frac{1}{\sqrt{2\pi }\sigma }\exp \left[-\displaystyle \frac{{\left({\omega }_{{nlm}}-{\omega }_{0}\right)}^{2}}{2{\sigma }^{2}}\right].\end{eqnarray} \tag{ 10 }$$

We tested extreme parameters, varying the peak of the Gaussian between 0 and 5 mHz, well above the acoustic cutoff frequency. We also tested widths between an extremely narrow distribution of 10 μHz and an extremely wide distribution of 5 mHz (see Figure 3). By varying the parameters so widely, we can capture a wide variety of possible frequency-dependent behavior. A Gaussian with a faraway peak behaves approximately like an exponential relationship. A Gaussian with an extremely narrow peak behaves approximately like a δ-function, and one with an extremely wide peak approximates equipartition. As discussed in Section 4, we find narrowly peaked distributions centered on low-order p-modes to be the most likely to fit the data, although there is considerable degeneracy within that region of parameter space. We settled on the Gaussian dependence on frequency after trying a variety of parametric models because it provided the best fit to the data, and was flexible in qualitatively approximating many diverse behaviors. We do not claim that the real power spectrum behaves in exactly this way.

Zoom In Zoom Out Reset image size

This frequency dependence on energy has a straightforward connection to excitation and dissipation efficiency, if the energy spectrum is in steady state. In this case,

$$\begin{eqnarray}&&E(\omega )=\displaystyle \frac{{\dot{E}}_{\mathrm{in}}(\omega )Q(\omega )}{\omega },\end{eqnarray} \tag{ 11 }$$

where Q(ω) is the frequency-dependent quality factor and ${\dot{E}}_{\mathrm{in}}(\omega )$ is the frequency-dependent excitation rate. The mode energy is ${E}_{{nlm}}={a}_{{nlm}}^{2}{\omega }_{{nlm}}^{2}{M}_{{nlm}}$, where M_nlm is the modal inertia uniquely defined for each normal mode according to ${M}_{{nlm}}=\int \rho | {\xi }_{{nlm}}{| }^{2}{d}^{3}r$, where ξ_nlm is the mode eigenfunction normalized such that the surface displacement is 1 cm. Therefore the energy scales as the square of the amplitude.

There is a fundamental ambiguity when modeling normal modes that does not exist for a static gravity field: the phase of the mode in question. If we only had one flyby, breaking this ambiguity would be hopeless; since we have multiple, we can do better. By combining multiple flybys, we can average out the effect of initial phases. This can be done perfectly if there is a large number of identical flybys. Indeed, one can show (see the Appendix) for a particle on a prescribed trajectory ${\boldsymbol{r}}(t)$ encountering a potential perturbation of the form of Equation (8), where ϕ₀ and α_nlm are random variables for each flyby, that the summed squared potential obeys the asymptotic relationship

$$\begin{eqnarray}&&\displaystyle \sum _{i}^{N}\delta {{\rm{\Phi }}}_{i}^{2}(t)\sim \displaystyle \frac{N}{2}\displaystyle \sum _{{nlm}}{A}_{{nlm}}^{2}{f}_{{nlm}}^{2}(t),\end{eqnarray} \tag{ 12 }$$

where ${A}_{{nlm}}\equiv {a}_{{nlm}}\delta {C}_{{nlm}}$, and ${f}_{{nlm}}(t)$ is a deterministic function time independent of ϕ₀ and α_nlm. This approximation is valid at large N. It is possible to derive a similar expression for acceleration perturbations, simply the gradient of the potential perturbations, and for velocity perturbations (see the Appendix). In fact these derivations depend on assuming a prescribed trajectory, but in reality the potential perturbations perturb the trajectory itself. We verify empirically that for a sufficiently large number of flybys

$$\begin{eqnarray}&&\displaystyle \sum _{i}^{N}\delta {v}_{i}^{2}=\displaystyle \frac{N}{2}\displaystyle \sum _{{nlm}}{a}_{{nlm}}^{2}\displaystyle \sum _{i}^{N}\delta {v}_{{nlm},i}^{2},\end{eqnarray} \tag{ 13 }$$

where $\delta {v}_{{nlm},i}^{2}$ is a randomly generated squared time series of velocity perturbation associated with the spacecraft encountering the potential perturbation due to a 1 cm displacement amplitude mode with quantum numbers (n, l, m).

We note that the stochastic behavior of the modes approaches deterministic behavior when summing over a large number of flybys to demonstrate why such an exercise is useful: it reduces the stochastic component of the signal and amplifies the deterministic component. In the real experiment, however, there were only five flybys, not enough to simply stack the data and compare them against the asymptotic average. Therefore, we ran a Monte Carlo simulation, leaving the initial phase as a free random variable, and combined the signal from five randomly selected flybys with an input spectrum scaled according to Section 2.3. In this case, we use the relationship from Equation (9) to obtain the expression

$$\begin{eqnarray}&&\displaystyle \sum _{i}^{N}\delta {v}_{i}^{2}=\displaystyle \sum _{i}^{N}{\left(\displaystyle \sum _{{nlm}}{a}_{{nlm}}\delta {v}_{{nlm}}\right)}^{2},\end{eqnarray} \tag{ 14 }$$

which is equivalent to Equation (13) in the limit of large N, but for finite N has stochastic components that should be accounted for.

With N = 5 we can eliminate significant ambiguity. The raw data are shown as the red scattered points in Figure 6 for each flyby (the black curves are model fits to the data). The raw data were obtained by subtracting the observed spacecraft signal from a model excluding stochastic acceleration (Iess et al. 2019). We take these points and bin them into 150 s windows so that most points in time will have contributions from all flybys (see Section 3.1 for why this is important). We then average the square value of the corresponding data point across the five flybys to obtain an average value. After accounting for various quantifiable sources of error, we produce Figure 5 that shows the stacked data with error bars in red, with a black curve as a good-fit forward model.

3.1. Error Sources

Now that we have added error bars to account for the straightforwardly quantifiable sources of error, we can attempt to fit them. We do this with a reduced χ² test according to

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{D-f}\displaystyle \sum _{i}^{D}\displaystyle \frac{{\left({x}_{i}-{m}_{i}\right)}^{2}}{{e}_{i}^{2}},\end{eqnarray} \tag{ 15 }$$

where D is the number of data points, f is the number of degrees of freedom in the model (three in our case: the two Gaussian parameters and the scaling coefficient), x_i is data point i, m_i is its corresponding modeled value, and e_i is the error. Choosing an appropriate model is subtle. One choice is to use the asymptotic average of an infinite number of flybys for a given spectral model. As demonstrated in Section 3, we expect the data to converge toward this average. But given the finite number of flybys, there will be variation from this asymptotic mean. Therefore in order to evaluate the likelihood of a given model, we conducted 2000 tests of five simulated flybys for each modeled spectrum.

To produce our forward model, we ran 10⁴ simulated signals from individual modes, for the subset of considered modes (recall our results are not sensitive to the specific choices of considered modes or computed eigenfrequencies. For more discussion see Section 2). After simulating a large number of range-rate signals for each individual mode, we chose five to combine their squared signal using Equation (14), where a_nlm for a given model is computed according to Equation (10) with ${E}_{\mathrm{nlm}}={a}_{\mathrm{nlm}}^{2}{\omega }_{\mathrm{nlm}}^{2}{M}_{\mathrm{nlm}}$. This is the forward model we use to try to fit Figure 5. Our tolerance threshold for goodness of fit is χ² = 50. This is a large value, but we consider it sufficient to qualitatively reproduce the essential shape of the data (for further discussion as to why this is appropriate, see Section 3.1). Choosing a different threshold does not significantly affect the results, but reduces the probability of fitting within the tolerance threshold for all models. We then use the large number of experiments to assign a probability of reproducing the data within tolerance for a given input spectrum. The results show a degenerate set of distributions that can reproduce the signal. We found a strongly favored region of parameter space after coarse sampling using extreme parameters, then followed up with a finer sampling in that region.

The probability plots are shown in Figure 3, which illustrates a clearly preferred region of parameter space, but degeneracy within that region. Each grid cell of Figure 3 represents a particular spectral model. The shading indicates the probability of reproducing the data within our tolerance threshold, if that spectrum were Saturn's normal-mode spectrum. The degeneracy can be partly understood by considering the contributions to the gravity signal a moving spacecraft encounters. The degeneracy can be understood in two ways. First, as shown in Figure 4, there is a great deal of overlap between favored models. Second, because the spacecraft is moving through space, a static field would have a time-dependent signature. Because normal modes oscillate in time, there is another source of time dependence that would be experienced even by a stationary test particle. The synthesis of these two contributions allows gravity perturbations with different properties to produce a similar signal along the Saturn–Earth line-of-sight axis from the spacecraft's frame of reference. We elaborate on this second degeneracy source in Section 2.2 and in the Appendix.

Zoom In Zoom Out Reset image size

Although we cannot identify a single conclusive power spectrum, we can exclude a wide variety of simple spectra, and find the highest-probability models favor a relatively narrowly peaked distribution. The location of the peak is also constrained, with the most likely models having a peak between 500 and 700 μHz.

Although the exact width and frequency peak cannot be precisely determined, we can exclude a wide variety of models as implausible, and note a clear clustering of models that have a high probability of reproducing the observed signal. Low-order f-modes lie generally below 200 μHz (see Figure 1). No good-fit models favor significant contributions from f-modes. We demonstrate an example of what we consider to be a "plausible" fit near the cutoff threshold in Figure 5. This particular run has χ² = 4, among the better fits we were able to obtain.

Zoom In Zoom Out Reset image size

We also assessed the frequency content of the residual data. This proved to be less diagnostic than fitting the time series. This is perhaps unsurprising, because in order to fit the frequency content, we only need a model that varies on the correct timescales, driven largely by the geometric ${\left(r/R\right)}^{l}$ effect. By contrast, in order to fit the time series data, we have to match much more specific behavior. Looking at Figure 5, the time series model must fit several specific phenomena. From left to right on the figure, the best models the steepness of the "ramp up" before C/A, the timing of the peak after C/A, the width of the main curve, the timing of the plateau/turnover, and other features. We find that the plausible time series fits are also compatible with the data's Fourier transform. However, the Fourier transform is much more degenerate and possible to fit with a wide variety of models, and it is difficult to obtain any new information.

To verify that our results were not excessively biased by our assumptions, we ran a variety of tests and alternatives, in addition to trying various parametric models as described in Section 2.3. We also explored the possibility that the signal may be dominated by a single mode, by testing that hypothesis against each mode in our sample. This possibility seems plausible based on our results given the narrowness of the peak in many best-fit cases. We found some modes within the preferred region of parameter space had a finite probability of reproducing the data, but the probability was lower than our preferred spectral models. Consistent with our spectral method, f-modes were not favored. No f-mode had a probability higher than 2% of producing the observed signal.

We separately tested all f-modes up to spherical order 20 for completeness, because of the a priori expectation that they should be the most gravitationally important modes, and because some of their frequencies overlap with the degenerate region of parameter space that can give some probability of providing a tolerable fit to the data. High degree f-modes are discussed in more detail in the Appendix. Even allowing for higher-order f-modes, we did not find any simple combinations that satisfactorily reproduce the data. Perhaps f-modes are inefficiently excited for reasons beyond their frequency; for example, some excitation models depend on compressibility (Goldreich et al. 1994), and in the Sun mode power declines for increasing l even at fixed eigenfrequency.

We also tested g-modes. Although our interior model assumptions did not produce g-modes, we tested them specially using published eigenfrequencies (Gudkova & Zharkov 2006), and separately testing asymptotic approximations for their eigenfrequencies (Tassoul 1980). Without the eigenfunctions, we tested the spectral model by varying the gravity harmonic coefficient perturbations directly with frequency dependence. We did not find any solution that could satisfactorily reproduce the data with g-modes. We also explored the possibility of using a given mode's eigenfrequency as a tunable parameter, varying our expectation for each f-mode's eigenfrequency between half and three times its theoretically predicted value. Although some frequencies fit better than others, none came close to the goodness of fit we obtain with our spectral model. We also tested the full pipeline omitting one flyby, testing each subset of four flybys to ensure the results were consistent, and not a spurious peculiarity of these five particular flybys. That is, we wanted to ensure that if one of the flybys had had a problem such that it did not successfully transmit data, it would not have altered our conclusion. As expected, with fewer flybys the preferred region of parameter space could not be as tightly constrained, but the results were consistent and favored the same region shown in Figure 3. If future missions can perform the same experiment with a larger number of flybys, we may be able to make stronger conclusions.

All plausible spectral models predict large peak mode amplitudes on the order of several kilometers for a small number of modes (of order 5–10) near the peak frequency. Mode amplitudes inferred from velocity map time series of Jupiter are of order 100 m (Gaulme et al. 2011), so in order to explain our findings we require the peak amplitudes to be at least an order of magnitude larger on Saturn than have been observed on Jupiter.

This method also allows us to fit the range-rate residuals from each individual flyby. We begin with a sample amplitude spectrum with a high probability of reproducing the data (see Figure 3). We then run a suite of simulations with random initial phases of each mode and show the best-fit results for each flyby in Figure 6. Each fit optimizes for the best-fit scaling coefficient, and all are in agreement within a factor of two.

Zoom In Zoom Out Reset image size

We can reproduce the behavior of the non-zonal and/or non-static component of Saturn's gravity field using a simple three-parameter forward model for mode amplitudes (Gaussian peak, width, and scaling). Our model uses a simple interior model and is not sensitive to detailed assumptions about Saturn's interior structure, spin rate, or rotation profile. For the amplitude spectrum of the modes, we rely only on a general understanding of Saturn's eigenfunctions and the equation of continuity in order to compute its gravitational effect for a given amplitude. We can also compute the mode inertia to scale the gravity signal for a given mode energy.

The best way to think of these results is in the Bayesian sense; we do not claim incontrovertible proof of p-modes on Saturn. Rather, given this particular data set, we present the probability that different models will reproduce these data in Figure 3. These probabilities should be used to update prior assumptions about Saturn's normal mode spectrum, bearing in mind that more complex models with more degrees of freedom (e.g., from Iess et al. 2019) are not captured by our analysis.

We find a moderate preference for models that have a frequency peak between about 500 and 700 μHz with a narrow width, although models with peaks as low as 250 μHz or as high as 1000 μHz also have nonzero probability. Intriguingly, the inferred narrowness of the peak is analogous to the narrowly peaked five minute modes observed on the Sun. Observing the power spectrum of the solar modes, one finds a peak frequency near 3000 μHz with a full width half maximum (FWHM) of order 300 μHz (Frohlich et al. 1997), indicating a ratio of the FWHM to peak frequency of order 1/10. Converting the standard deviation of our Gaussian frequency-dependent functions to the corresponding FWHM, we likewise obtain a solution of order 1/10 (for example in the case of the fit shown in Figure 5). This may indicate some similarity between the two systems, for example that the peak frequency is set by some dynamical process with a characteristic timescale. On the Sun this timescale is the eddy turnover time in the top scale height of the convective zone. Although the same mechanism cannot excite the observed amplitudes on Saturn, less frequent moist convective events with a characteristic turnover timescale (see, e.g., Markham & Stevenson 2018) could produce similar strong frequency dependence. On Jupiter, new theories to explain the ammonia distribution require updrafts that traverse 100 km in 1000 s (Guillot et al. 2020). If similar dynamics occur on Saturn, the timescale is roughly consistent with the peak frequencies inferred by this work. Others have suggested a large impact as a source of Saturn's oscillations (Wu & Lithwick 2019). Although a 150 km impactor could in principle excite kilometer-scale oscillations in p-modes, the scaling suggests the gravity signal from f-modes should always dominate. Therefore, if this is indeed the dominant excitation mechanism on Saturn, there must be some other reason to preferentially dissipate f-modes or preferentially amplify p-modes.

Most notably, we cannot reproduce the time series data with f-modes, either with a single f-mode dominating the signal or with a straightforward superposition of f-modes. This finding is consistent with inferred amplitudes of f-modes that have been measured using Saturn's rings (kronoseismology; Hedman & Nicholson 2013; Fuller 2014; Wu & Lithwick 2019), which are determined to be on the order of a meter in amplitude and should not produce this large a signal in Cassini's gravity experiment. Assuming the amplitudes inferred from ring data, the detectability should have been marginal (δv ∼ 0.05 mm s⁻¹ instead of the observed ∼2 mm s⁻¹). The p-modes that are required to produce the observed signal would need surface amplitude on the order of kilometers, implying radial velocities of meters per second. These p-modes, despite their large amplitudes, are not expected to show structure in the rings, because the relevant resonant radius for these frequencies is well inside Saturn's C ring.

The required p-mode amplitudes are at least an order of magnitude larger than were observed on Jupiter (Gaulme et al. 2011). Interestingly, early analysis of the Juno mission indicates a similar unexplained gravity signal on Jupiter that is approximately 20 times weaker than the signal observed on Saturn (Durante et al. 2020). This is interesting, because this analysis indicates that if the relevant amplitudes were those observed by Gaulme et al. (2011), then we should expect a similar time-dependent signal diminished in scale by about an order of magnitude. Replicating our analysis of Cassini's gravity data for Juno, which has many more planned gravity orbits than Cassini, may be a promising future application of the method outlined in this paper. We can test to see if the inferred normal mode spectrum from gravity measurements on Jupiter is consistent with the corresponding power spectrum obtained with Earth-based observations.

Because we predict large peak amplitudes, we must consider whether these are plausible and consistent with existing data. Voyager radio occultation measurements of Saturn have error estimates between 6 and 10 km, and the measurements found incompatible radii between the northern and southern hemisphere on the order of 10 km (Lindal et al. 1985). These uncertainties are compatible with the time-dependent shape variations our analysis predicts. Our analysis here would predict that future measurements of Saturn's shape cannot obtain better accuracy than around a few kilometers. A series of highly accurate measurements of Saturn's shape should have a time-dependent component on the order of kilometers.

We must also consider how our findings can constrain Q. Using a set of N modes excited to 10 km amplitude a (a high estimate, see Table 1), powered by the full luminosity of Saturn ${ \mathcal L }$ as a lower bound on Q, we compute

$$\begin{eqnarray}&&Q\gt \displaystyle \frac{{a}^{2}{\omega }^{3}{MN}}{{ \mathcal L }}\sim {10}^{7}\end{eqnarray} \tag{ 16 }$$

for the relevant peak modes we have identified, where M is the modal mass and N ∼ 5–10 corresponds to the number of modes with significant amplitude (for example, greater than a kilometer). This Q is compatible with estimates so far based on theory (Markham & Stevenson 2018; Wu & Lithwick 2019).

We must also check that the modes can still be approximated as linear perturbations, i.e., ${\boldsymbol{u}}\cdot {\rm{\nabla }}{\boldsymbol{u}}\ll \tfrac{d{\boldsymbol{u}}}{{dt}}$, where u is the velocity vector. For the peak modes we identified, the frequency is sufficiently high that the motion is almost purely vertical. $\tfrac{{du}}{{dz}}$ is most significant near the surface when the atmospheric properties vary quickly. We know in this region $\tfrac{{du}}{{dz}}\sim \left(1-\sqrt{1-{\left(\tfrac{\omega }{{\omega }_{c}}\right)}^{2}}\right)u/H$ (Christensen-Dalsgaard 2014), where H is the scale height and ω_c is the acoustic cutoff frequency. Therefore the condition that the system can be treated linearly is

$$\begin{eqnarray}&&\left(1-\sqrt{1-{\left(\displaystyle \frac{\omega }{{\omega }_{c}}\right)}^{2}}\right)\ll \omega H/u\end{eqnarray} \tag{ 17 }$$

the ratio of the left-hand side to the right-hand side using ω ∼ 5 × 10⁻³ for 1 km amplitude modes is about 10⁻³. This is the maximum value near the surface; the value is much smaller in the interior where most of the mode inertia is. If the inferred amplitudes are correct, nonlinear effects are important for p-modes on both Saturn and Jupiter.

Many questions remain, and it is clear that the field of giant planet seismology—both observational and theoretical—is in its infancy. Here we demonstrate that the unexpected and unexplained components of Saturn's "dark" gravity field can be straightforwardly modeled as simple frequency-dependent seismic activity. This provides one more piece of plausible evidence that the giant planets are seismically active, and should motivate further observations and theoretical study.

A.1. Potential Perturbation

In the following derivation, we will assume the spacecraft is on a prescribed Keplerian orbit, and assess the gravity potential field it encounters as a function of time. In reality, the data we have are a velocity time series; we cannot measure the gravity potential directly. However, we demonstrate how a gravity potential with stochastic elements can be averaged toward deterministic behavior that isolates information about the amplitude spectrum. This same basic procedure will be employed to forward model the velocity time series, although the details in that case are considerably more complicated.

We begin with Equation (8). The full potential of all modes for a given flyby i is

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{i}(t) & = & \displaystyle \sum _{{nlm}}{A}_{{nlm}}{\left(\displaystyle \frac{R}{r(t)}\right)}^{l+1}{P}_{{lm}}(\cos (\theta (t)))\\ & & \times \,\cos [m(\phi (t)-{\rm{\Omega }}t-{\phi }_{i})-{\omega }_{{nlm}}t-{\alpha }_{{nlm},i}]\end{array}\end{eqnarray} \tag{ A1 }$$

with ${A}_{{nlm}}\equiv {a}_{{nlm}}\delta {C}_{{nlm}}$. ϕ_i is the random initial longitudinal orientation of Saturn (the same for each mode, but random for each orbit because Saturn's spin rate is not precisely known), and ${\alpha }_{{nlm},i}$ is the temporal phase of each mode—random for each mode and for each flyby because we do not know the eigenfrequency with sufficient precision to impose phase coherency between subsequent close encounters. The uncertainty in ϕ₀ can be absorbed into α; for zonal modes, ϕ₀ does not matter, and for tesseral/sectoral modes it can be added into α so that there is only one relevant random variable. We can rewrite this expression using the harmonic addition theorem:

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{i}(t) & = & \displaystyle \sum _{{nlm}}{A}_{{nlm}}{f}_{\mathrm{lm}}(t)\\ & & \times \,\left(\cos ({g}_{m}(t))\cos {\alpha }_{{nlm},i}-\sin ({g}_{m}(t))\sin {\alpha }_{{nlm},i}\right),\end{array}\end{eqnarray} \tag{ A2 }$$

where ${f}_{\mathrm{lm}}(t)={\left(\tfrac{R}{r(t)}\right)}^{l+1}{P}_{{lm}}(\cos (\theta (t)))$ and ${g}_{m}(t)=m(\phi (t)\,-{\rm{\Omega }}t)-\omega t$ is the same for every orbit and does not depend on the random variable ${\alpha }_{{nlm},i}$. Now we square this expression

$$\begin{eqnarray}\begin{array}{rcl}{{\rm{\Phi }}}_{i}{\left(t\right)}^{2} & = & \displaystyle \sum _{q}\displaystyle \sum _{q^{\prime} }{A}_{q}{A}_{q^{\prime} }{f}_{q}(t){f}_{q^{\prime} }(t)\\ & & \times \,(\cos ({g}_{q}(t))\cos ({g}_{q^{\prime} }(t))\cos {\alpha }_{q,i}\cos {\alpha }_{q^{\prime} ,i}\\ & & -\,\cos ({g}_{q}(t))\sin ({g}_{q^{\prime} }(t))\cos {\alpha }_{q,i}\sin {\alpha }_{q^{\prime} ,i}\\ & & -\,\sin ({g}_{q}(t))\cos ({g}_{q^{\prime} }(t))\sin {\alpha }_{q,i}\cos {\alpha }_{q^{\prime} ,i}\\ & & +\,\sin ({g}_{q}(t))\sin ({g}_{q^{\prime} }(t))\sin {\alpha }_{q,i}\sin {\alpha }_{q^{\prime} ,i}),\end{array}\end{eqnarray} \tag{ A3 }$$

where we substitute a single index q to refer to a given mode (n, l, m) for notation convenience.

From here, we perform the crucial step of summing over many such flybys. This total can be expressed as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \sum _{i=1}^{N}{{\rm{\Phi }}}_{i}{\left(t\right)}^{2} & = & \displaystyle \sum _{i=1}^{N}\displaystyle \sum _{q}\displaystyle \sum _{q^{\prime} }{A}_{q}{A}_{q^{\prime} }{f}_{q}(t){f}_{q^{\prime} }(t)\\ & & \times \,(\cos ({g}_{q}(t))\cos ({g}_{q^{\prime} }(t))\cos {\alpha }_{q,i}\cos {\alpha }_{q^{\prime} ,i}\\ & & -\,\cos ({g}_{q}(t))\sin ({g}_{q^{\prime} }(t))\cos {\alpha }_{q,i}\sin {\alpha }_{q^{\prime} ,i}\\ & & -\,\sin ({g}_{q}(t))\cos ({g}_{q^{\prime} }(t))\sin {\alpha }_{q,i}\cos {\alpha }_{q^{\prime} ,i}\\ & & +\,\sin ({g}_{q}(t))\sin ({g}_{q^{\prime} }(t))\sin {\alpha }_{q,i}\sin {\alpha }_{q^{\prime} ,i}).\end{array}\end{eqnarray} \tag{ A4 }$$

In this case, we can assume random variables behave as true statistical averages. We make use of the fact that ${\sum }_{q,q^{\prime} ,i,j}\cos {\alpha }_{q,i}\cos {\alpha }_{q^{\prime} ,j}\to {\delta }_{{qq}^{\prime} }\delta {ij}$. Similarly ${\sum }_{q,q^{\prime} ,i,j}\cos {\alpha }_{q,i}\sin {\alpha }_{q^{\prime} ,j}\to 0$. Using this asymptotic behavior, we can write for a sufficiently large number of N flybys,

$$\begin{eqnarray}&&\displaystyle \sum _{i}^{N}{{\rm{\Phi }}}_{i}^{2}(t)\sim \displaystyle \frac{N}{2}\displaystyle \sum _{{nlm}}{A}_{{nlm}}^{2}{f}_{{nlm}}^{2}(t).\end{eqnarray} \tag{ A5 }$$

Thus a particle on a prescribed trajectory encountering potential perturbations due to normal modes will, when averaging over many such encounters, approach a deterministic curve that does not depend on the initial phase of each mode.

A.2. Acceleration and Velocity

The gravity experiment does not directly measure the gravity potential. The data we have are the velocity perturbation along a single axis. Our data-stacking method is most straightforwardly derived for gravity potential perturbations. To calculate gravitational acceleration we apply $\delta {g}_{{nlm}}=-{\rm{\nabla }}\delta {{\rm{\Phi }}}_{{nlm}}$ to obtain

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\delta {{\boldsymbol{g}}}_{{nlm}}}{{A}_{{nlm}}}=\displaystyle \frac{{GM}}{{r}^{2}}{\left(\displaystyle \frac{R}{r}\right)}^{l}\\ \ \times \,\left[-(l+1)\cos [m(\phi -{\rm{\Omega }}t)-\omega t-\alpha ]{P}_{l}^{m}(\cos \theta )\hat{{\boldsymbol{r}}}\right.\\ \ +\,\cos [m(\phi -{\rm{\Omega }}t)-\omega t-\alpha ]\displaystyle \frac{\partial {P}_{l}^{m}(\cos \theta )}{\partial \theta }\hat{\theta }\\ \ -\,\left.\displaystyle \frac{m}{\sin \theta }\sin [m(\phi -{\rm{\Omega }}t)-\omega t-\alpha ]{P}_{l}^{m}(\cos \theta )\hat{\phi }\right].\end{array}\end{eqnarray} \tag{ A6 }$$

This gives us the components we need to project onto $\hat{\oplus }$, the unit vector pointing toward Earth. The total gravitational acceleration perturbation from normal modes can be written

$$\begin{eqnarray}&&\delta {g}_{\oplus }(t)=\displaystyle \sum _{{nlm}}{A}_{{nlm}}{g}_{\oplus ,{nlm}}(t),\end{eqnarray} \tag{ A7 }$$

where we can express

$$\begin{eqnarray}&&{g}_{\oplus ,{nlm}}(t)={A}_{{nlm}}\left({f}_{1}(t)+{f}_{2}(t)\right)({f}_{3}(t)\cos \alpha +{f}_{4}(t)\sin \alpha ),\end{eqnarray} \tag{ A8 }$$

where f_i(t) are tedious but nevertheless well defined functions of time independent of α. Just as above, if we have an expression in this form, we can express the asymptotic average of the sum of a large number of flybys as

$$\begin{eqnarray}&&\displaystyle \sum _{i=1}^{N}{g}_{\oplus ,i}^{2}(t)=\displaystyle \frac{N}{2}\displaystyle \sum _{{nlm}}{A}_{{nlm}}^{2}{f}_{{nlm}}{\left(t\right)}^{2},\end{eqnarray} \tag{ A9 }$$

where f(t)² is a function of f₁(t) and f₂(t).

A similar procedure can be followed for velocity perturbations, and this is the source of frequency dependence (absent in averaged potential perturbations) ${\omega }_{{nlm}}$ by temporally integrating the gravitational perturbations. In fact it is not this simple; a perturbed spacecraft will deviate from its Keplerian trajectory, which leads to errors from the exact solution of order $({\rm{\Delta }}r/r){\rm{\Delta }}t$. A better approximation is to dynamically solve the equation of motion, accounting for displacement from the initial Keplerian trajectory to linear order, although this method produces higher-order errors from the exact solution. The method we eventually used, solving the exact equation of motion explicitly then subtracting the Keplerian solution, is not conducive to an analytic expression. Because of this, we verified these results by testing them numerically against simulations to be sure they did approach a deterministic curve when stacked.

To do this, we performed 10⁴ Monte Carlo simulations of random flybys for each mode. We then averaged random subsets of these samples to verify that they asymptotically approach an asymptotic curve. We also directly simulated various examples of a spectral superposition of different modes, running 10⁴ flybys. We then compared these stacked results against the superposition of a spectrum of averaged squared modes, finding excellent agreement with Equation (13). We therefore have an analytic approximation that motivates the stacking procedure, as well as more exact numerical tests to make sure the results from the analytic approximation are robust to this application. These tests verify the validity of our data-stacking procedure. In practice, we have a finite number of flybys, N = 5. A random set of five flybys will deviate somewhat from the asymptotic behavior of 10⁴ flybys, so for the actual probabilistic fitting routine we used Equation (14). We used the statistical method outlined in Section 4 to account for this, by directly solving for the probability that a given spectrum will produce a good fit to the data.

A.3. Spacecraft Sensitivity and Model Intuition

We must consider whether the particular orbit of Cassini biases its ability to detect certain modes. If the convolution of a particular mode's space-dependent eigenfunction and time-dependent eigenfrequency appears stationary in Cassini's frame during close approach, then Cassini will preferentially detect signals from this mode. For example, because Cassini's orbit is nearly polar, there may be an intrinsic preference for sectoral (m = ±l) modes, as seems to be true for f-modes (see Figure A1). A different example is if a mode's half-period is near the timescale of Cassini's motion from the northern to the southern hemisphere, which appears to be the case for tesseral modes with m = l − 1 in Figure A1. Moreover, because Cassini's orbit moves from west to east, it is plausible that it would preferentially detect m > 0 modes whose pattern rotates in the prograde direction, an intuition also supported by Figure A1. According to Figure A1, this quasi-resonant effect does indeed make Cassini more sensitive to certain modes. All these effects are implicitly accounted for in the forward model in the main text.

Zoom In Zoom Out Reset image size

Another important consideration is high-degree f-modes. If the strong frequency dependence we predict applies equally to f-modes, then we must consider whether Cassini would detect them. There are two sources of attenuation: the intrinsic geometric attenuation with distance of high-degree gravity harmonics (R/r)^l, and the monotonically decreasing gravity potential coefficient response to each mode that obeys ${(l(l+1))}^{-1/2}$ (Wu & Lithwick 2019). We use a characteristic radius for Cassini's close encounter ${r}_{c}=v/R{\int }_{-R/2v}^{R/2v}r(t){dt}\sim 1.15R$, where v is the periapse velocity. High-degree f-mode eigenfrequencies obey ${\omega }^{2}\sim \tfrac{{GM}}{{R}^{3}}{\left(l(l+1)\right)}^{1/2}$, so the modes in the most likely peak frequency region have l ∼ 50–100. Using our attenuation estimates, these modes would be reduced by at least five orders of magnitude compared to low-order f-modes, significantly more severe than the two–three orders of magnitude reduction in low-order p-modes. Therefore, we do not expect Cassini to detect high-order f-modes, even if they likewise had kilometer-scale amplitudes.