<figure><img src="image_1.png" /></figure>

# Long-Term Evolution of the Saturnian System

Matija Ćuk¹ ID • Maryame El Moutamid² ID • Giacomo Lari³ ID • Marc Neveu⁴,⁵ ID • Francis Nimmo⁶ ID • Benoît Noyelles⁷ ID • Alyssa Rhoden⁸ ID • Melaine Saillenfest⁹ ID

Received: 1 May 2023 / Accepted: 18 January 2024
© The Author(s) 2024

## Abstract

Here we present the current state of knowledge on the long-term evolution of Saturn's moon system due to tides within Saturn. First we provide some background on tidal evolution, orbital resonances and satellite tides. Then we address in detail some of the present and past orbital resonances between Saturn's moons (including the Enceladus-Dione and Titan-Hyperion resonances) and what they can tell us about the evolution of the system. We also present the current state of knowledge on the spin-axis dynamics of Saturn: we discuss arguments for a (past or current) secular resonance of Saturn's spin precession with planetary orbits, and explain the links of this resonance to the tidal evolution of Titan and a possible recent cataclysm in the Saturnian system. We also address how the moons' orbital evolution, including resonances, affects the evolution of their interiors. Finally, we summarize the state of knowledge about the Saturnian system's long-term evolution and discuss prospects for future progress.

**Keywords** Satellites of Saturn • Tidal evolution • Orbital resonances • Obliquity of Saturn

## 1 Tidal Dissipation Within the Planet and Orbital Evolution

Satellites raise a tidal bulge on the planet they orbit. If there is any kind of lag in the planet's response, the tidal bulge will not point directly towards the satellite and as a result, the latter will experience a torque. Generally, planets spin faster than satellites orbit and, as a result, the torque pushes the satellites outwards. In the limit of low-eccentricity and inclination, a spin-synchronous satellite, and static solid tidal dissipation, for a tidal bulge averaged over an orbital cycle and apsidal precession, the semi-major axis increases as (Murray and Dermott 1999)

$$ \dot{a} = 3 \frac{k_{2,p}}{Q_p} \sqrt{\frac{G}{M_p}} R_p^5 m_s a^{-11/2}, \quad (1) $$

where G is the gravitational constant, $M_p$, $R_p$ and $k_{2,p}$ are the planet's mass, radius and tidal Love number, respectively, while $m_s$ and $a$ are the mass and semi-major axis of the satellite that raises the tides. Traditionally, the lag has been assumed to be proportional to $Q_p^{-1}$, where $Q_p$ is the so-called dissipation factor of the planet and is often taken to be constant.

Extended author information available on the last page of the article
In this theory of “equilibrium tides”, the semi-major axis evolves according to $t^{2/13}$, where $t$ is time (Murray and Dermott 1999) and more distant satellites tend to evolve outwards more slowly. The result is converging orbits, which raises the possibility of capture into mean-motion resonance (see Sect. 2). Tides on the planet also cause the eccentricity of satellites to increase, even though this effect is generally smaller than the eccentricity damping due to the tidal dissipation within the moons (see e.g. Goldreich and Gold 1963).

More recently, a different mechanism, termed “resonance locking”, has been proposed (Fuller et al. 2016). This mechanism is described in more detail in Fuller et al. (2024, this collection). The key results, however, are as follows: 1) orbital migration can be significantly more rapid than previously thought; 2) outward evolution is exponential or power-law in time, rather than going as $t^{2/13}$, meaning that the rate of outwards motion accelerates; and 3) more distant satellites tend to evolve outwards faster. All three statements are true only if a satellite is in a resonance lock, which is not assured and depends on the internal structure of the planet.

Point (1) has received some observational confirmation from astrometry and radio science (Lainey et al. (2020) but cf. Jacobson (2022)). Point (2) means that an observed rapid present-day outward migration does *not* necessarily imply that the satellite formed recently. Point (3) means that the opportunities for capture into resonance are more limited than in the traditional picture, though the results are sensitive to very poorly known parameters.

Another important consequence of the resonance-locking hypothesis is that for satellites in a mean-motion resonance, the tidal heat production rate corresponding to constant (equilibrium) orbital eccentricities may be much higher than previously thought, because the effective $Q_p$ of the planet is lower (e.g. Meyer and Wisdom 2007). This issue is discussed in more detail in Nimmo et al. (2023, this collection), and in Sects. 3 and 4.2 below.

## 2 Mean-Motion Resonances

A mean-motion resonance (MMR) arises when two orbiting bodies (in this case satellites of a planet) have a ratio of orbital mean motions that is close to a simple fraction, i.e. $1/2$, $2/3$, $3/5$ etc. MMRs are usually locations of strong gravitational interaction between the moons, and depending on many factors, this interaction could be temporary or the bodies may be permanently *captured into the resonance*, meaning that their period ratio stays fixed as their orbits evolve.

More rigorously, two bodies are in resonance if a relevant angle called the *resonant argument* is librating rather than circulating. For a MMR, the general form of the resonant argument is $j_1\lambda_1 + j_2\lambda_2 + j_3\omega_1 + j_4\omega_2 + j_5\Omega_1 + j_6\Omega_2$, where $\lambda, \omega$ and $\Omega$ respectively designate the mean longitude, the longitude of pericenter and the longitude of the ascending node (Murray and Dermott 1999), and the subscripts 1 and 2 refer to the inner and the outer satellite. The coefficients $j_{1-6}$ are integers, with the condition that $j_1 + j_2 + j_3 + j_4 + j_5 + j_6 = 0$ (D'Alembert's rule, arising because the geometry of MMR does not depend on the origin of longitude; Murray and Dermott 1999). The order of the resonant argument is $|j_3| + |j_4| + |j_5| + |j_6|$, and the corresponding resonant term in the Hamiltonian (i.e. interaction potential) will be multiplied by $e^{j_3|j_3|}e^{j_4|j_4|}s_1^{|j_1|s_1^{|j_1|}}s_2^{|j_2|s_2^{|j_2|}}$, where $e$ and $s = \sin(i/2)$ refer to orbital eccentricity and inclination. As $e < 1$ and $s < 1$, and especially because the orbits of the major satellites of the giant planets tend to have small eccentricities and inclinations, lower-order resonances are typically stronger than higher-order ones. Note that each period commensurability ($j_2 : j_1$) contains many formal resonances (often
called “sub-resonances”) with different resonant arguments, which may affect either inclination or eccentricity. For example, within the 2:1 commensurability between Mimas (“M”) and Tethys (“Θ”), we can have 1st order resonance (e.g. $2\lambda_{\Theta} - \lambda_{M} - \omega_{M}$), 2nd order resonance ($4\lambda_{\Theta} - 2\lambda_{M} - \Omega_{M} - \Omega_{\Theta}$), which is currently librating) as well as 3rd order resonance (e. g. $2\lambda_{\Theta} - \lambda_{M} + \omega_{M} - \Omega_{M} - \Omega_{\Theta}$). In general, $|j_1 + j_2|$ is the lowest order that a resonance within the commensurability can have, but higher order sub-resonances are always also present. Note that $j_5 + j_6$ has to be an even number (due to inclinations being relative, rather than absolute, and the geometry of the MMR not depending on the choice of reference plane; Murray and Dermott 1999), with the implication that resonances affecting inclination have to be at least of second order, while eccentricity-affecting resonances can be of first order.

Satellite systems of the giant planets introduce two elements in the resonant dynamics that are not as prominent for heliocentric orbits. Major satellites experience significant tidal torques from the primary, which makes moons migrate, potentially encountering resonances with each other (see Sect. 1). Additionally, the highly oblate shapes of the giant planets strongly perturb the orbits of the satellites, inducing fast orbital precession. Different sub-resonances of the same commensurability are separated in frequency space by secular (i.e. precessional) frequencies, so fast orbital precession leads to well-separated sub-resonances that do not overlap in frequency (e.g. Meyer and Wisdom 2008). This non-overlap is important in order for a stable resonance capture to occur, as resonance overlap produces dynamical chaos (Wisdom 1980), which sometimes breaks resonances and sometimes leads to a much more complex evolution (as seen in the Uranian system; Dermott et al. 1988; Ćuk et al. 2020).

One fundamental rule of resonant encounters is that convergent orbital evolution between two bodies can lead to resonance capture, while divergent crossings of resonances always lead to “jumps”, i.e. passages through resonance during which relevant orbital elements experience discontinuous “kicks” but no long-term resonant evolution (e.g. Dermott et al. 1988). As the resonances within the same commensurability are separated by precession frequencies, it is always the case that inclination-type resonances (i.e. those with resonant arguments including longitudes of the node) are encountered first during convergent migration, as the precession of orbital nodes is generally retrograde (against the direction of orbital motion). Eccentricity-type resonances (which include the longitudes of pericenter in the resonant argument) tend to be crossed by converging satellites later as the pericenters precess in the direction of orbital motion. Whether convergent migration results in capture into a specific resonance depends on the parameters of the bodies and the resonance (lower-order resonances and higher masses make capture more likely), as well as migration rates (slower migration enhances capture probability) and pre-resonance orbital elements (low eccentricity or inclination are required for capture; Murray and Dermott 1999). While dissipation within satellites does not directly affect the capture into resonances, it can indirectly enable capture by damping previously existing eccentricity (or inclination).

The above discussion assumes that a MMR is only affecting the two satellites whose orbits are in near-commensurability. Close proximity of several sizeable satellites in the Saturnian system means that three-body resonances are sometimes present. Most relevant three-body resonances are semi-secular ones involving Titan, in which the longitudes of the apse or node of Titan (the most massive moon and the most important perturber) appear in a resonant argument of a MMR between two moons other than Titan (Ćuk and El Moutamid 2022). Three-body resonances independent of two-body resonances are also possible, but tend to be relatively weak and any capture tends to be temporary (Ćuk and El Moutamid 2022).
Once two moons are captured into a resonance, their semi-major axes evolve in unison in order to preserve the resonant relationship. In a stable MMR, eccentricity or inclination of at least one of the moons has to increase over time. This eccentricity or inclination is determined by the resonant argument, with the increasing orbital element being a conjugate of the secular angle in the resonant argument (Murray and Dermott 1999). For example, in the current Titan-Hyperion resonance with the resonant argument $4\lambda_H - 3\lambda_T - \omega_H$ the secularly increasing element is $e_H$, conjugate to $\omega_H$. Higher-order resonances can have orbital elements of both moons increase over time, as is the case with the inclinations of both Mimas and Tethys in the current Mimas-Tethys 4:2 resonance. The eccentricities and/or inclinations of bodies in MMR are expected to grow indefinitely, unless balanced by dissipation in the satellite or more usually until some other dynamical mechanism arrests further evolution or even breaks the resonance. Often this is due to the onset of resonance overlap. Even if closely-spaced sub-resonances are separate at the moment of resonant capture, their widths are proportional to eccentricities and inclinations and become wider over time. Eventually this can produce resonance overlap, resulting in dynamical chaos and likely breaking of the resonance relation.

One special case in which the orbital eccentricity can be constant indefinitely in the resonance is when the tidal dissipation within the satellite (described in the next subsection) acts to damp the eccentricity. Depending on the masses of the moons, their rates of tidal evolution and the strength of dissipation within the affected moon, the system can settle in an equilibrium in which the increase of eccentricity within the resonance is exactly counteracted by the damping by satellite tides (Meyer and Wisdom 2007). It is widely thought that Enceladus is currently in such equilibrium, in which eccentricity pumping by 2:1 MMR with Dione is counteracted by tidal dissipation within Enceladus, keeping eccentricity constant and providing energy for the geothermal activity on Enceladus (Lainey et al. 2012). It is important to note that this situation can be stable only for some combinations of tidal dissipation parameters that damp the librations of the resonant argument, as in other cases the equilibrium may be temporary as librations grow and the moons evolve out of the resonance (Meyer and Wisdom 2008).

## 3 Tidal Dissipation and Heating Within Satellites

A synchronous satellite (i.e. one that has identical rotational and orbital periods) will experience a permanent tidal bulge as a consequence of its proximity to its parent planet. However, if the satellite's orbit is eccentric, the size and orientation of this tidal bulge will fluctuate; the latter is also true if the satellite has an axial tilt (obliquity). If there is any friction in the system, some of the mechanical energy deforming the satellite will be converted to thermal energy. This is the origin of tidal heating.

Quantitatively, the heat production rate $\dot{E}$ is given by (e.g. Wisdom 2008)

$$ \dot{E} = \frac{3}{2} \frac{n^5 R^5 k_2}{G Q} (7e^2 + \sin^2 \theta). \qquad (2) $$

Here $n$ is the satellite mean motion, $R$ is its radius, $e$ its eccentricity and $\theta$ its obliquity. This expression is correct in the limit of small $e$ and $\theta$. The quantity $k_2$ is the satellite's tidal Love number, which describes the response of its gravity field to the perturbing tidal potential; the maximum value (for a uniform, fluid body) is 1.5 and it is reduced if the body has appreciable rigidity or central concentration of mass. The quantity $Q$ describes the friction
in the satellite; $Q^{-1}$ is proportional to the lag angle between the applied potential and the response, so that a low $Q$ implies a large lag angle and correspondingly higher dissipation.

Calculation of $k_2$ and (especially) $Q$ for satellites<fnref n="1" /> is not straightforward. In the simplest case, a solid viscoelastic (Maxwellian) body, the density, rigidity and viscosity at each point in the satellite need to be specified to derive $k_2/Q$ (e.g. Moore and Schubert 2000). Since viscosity is strongly temperature-dependent, and the heat production rate depends on $k_2/Q$, even this very simple case raises the possibility of feedbacks between heat production and internal structure (Ojakangas and Stevenson 1986; Hussmann and Spohn 2004).

The Maxwell description of viscoelasticity, though often used, fails to adequately describe how real materials respond to forcing at different frequencies, so more complicated models (e.g. Andrade, Burgers) are needed (Renaud and Henning 2018; Renaud et al. 2021). If oceans are present, other modes of tidal heating are possible, such as turbulent dissipation in the oceans themselves (Chen et al. 2014), or tidal flushing of water through a permeable crust (Rovira-Navarro et al. 2022). However, the contribution of these mechanisms to the overall tidal heating budget appears in general to be small.

Because tidal heating represents an energy loss, the semi-major axis of the satellite must contract. If it is not experiencing external torques, for angular momentum to be conserved the eccentricity must also decrease. Thus, for an isolated satellite any eccentricity will be damped relatively rapidly. The damping rate for small $e$ is given by (Murray and Dermott 1999)

$$ \dot{e} = -\frac{a\dot{E}}{eGMm}, \qquad (3) $$

where $a$ is the semi-major axis, $m$ and $M$ the mass of the satellite and primary, respectively, and $\dot{E}$ is the tidal heating rate (equation (2)).

If eccentricity tides dominate then Eq. (3) can be used to derive the eccentricity damping timescale $\tau_e = e/\dot{e}$:

$$ \tau_e = \frac{GMme^2}{a\dot{E}}. \qquad (4) $$

Since $\dot{E}$ depends on $e^2$ (Eq. (2)), the damping timescale is independent of eccentricity.

If dissipation is dominated by obliquity rather than eccentricity tides, the inclination ($i$) damping timescale $\tau_i$ for small inclination is given by Chyba et al. (1989)

$$ \tau_i = \frac{GMm \sin^2 i}{a\dot{E}}. \qquad (5) $$

Because the satellites are expected to be in a Cassini state,<fnref n="2" /> the obliquity $\theta$ will scale with inclination and thus $\tau_i$ is independent of inclination. Other things being equal, inclination damping is slower than eccentricity damping because of the factor of $7(\sin i/\sin\theta)^2$ from Eqs. (2), (4) and (5). Note that the factor $(\sin i/\sin\theta)$ is typically $\gg 1$ for relatively close-in satellites of giant planets that are in the Cassini State 1 (Chen et al. 2014). However, for satellites with oceans, obliquity tidal heating can be much more effective than eccentricity

<footnote n="1">¹Not to be confused with the analogous parameter $Q_p$ for the planet; see Sect. 1.</footnote>

<footnote n="2">²Cassini states are the end-points of energy damping within a rotating satellite, in which the obliquity of the satellite is determined solely by its shape and orbital inclination, and any initial free obliquity is damped (Ward 1975).</footnote>
tidal heating (Tyler 2008; Chen et al. 2014), and so inclination damping may be more rapid than eccentricity damping.

Equation (3) is the reason that orbital resonances (see Sect. 2) are so important; although the eccentricity of an isolated satellite (like Triton) will swiftly damp to zero, if the eccentricity is being excited by an orbital resonance, then long-lasting tidal heating can arise. In this case the ultimate source of energy is the rotational kinetic energy of the planet. Because the damping rate depends on the tidal heating rate, and the tidal heating rate depends on the thermal structure, complicated and non-monotonic orbital and thermal histories can result (e.g. Ojakangas and Stevenson 1986; Hussmann and Spohn 2004).

Because torques from the planet tend to increase the eccentricity and semi-major axis (Sect. 1), while dissipation in the satellite decreases both, for satellites in a MMR an equilibrium can result in which $e$ is constant (see Sect. 2). Meyer and Wisdom (2007) show that the power of tidal heating of the inner moons of a resonant pair, assuming no tidal torque on the outer moon and all eccentricities being in equilibrium, is:

$$H = n_E T_E - \frac{T_E}{L_E + L_D} \left( \frac{Gm_S m_E}{a_E} + \frac{Gm_S m_D}{a_D} \right), \quad (6)$$

where $T_E$ is the tidal torque on Enceladus, while $L$ is angular momentum (= $ma^2n \times \sqrt{1-e^2}$), $n$ is mean motion, $a$ is semi-major axis, $m$ is mass, and $G$ is gravitational constant (subscripts E, D and S refer to Enceladus, Dione and Saturn, respectively). If we ignore terms of order $e^2$ (Meyer and Wisdom 2007):

$$H = n_E T_E \left( 1 - \frac{1 + m_D a_E / (m_E a_D)}{1 + (m_D / m_E) \sqrt{a_D / a_E}} \right). \quad (7)$$

The torque on Enceladus is given by:

$$T_E = \frac{3 G m_E^2 R_p^5 k_{2p}}{2 a_E^6 Q_p}. \quad (8)$$

In this special case, the dissipation rate depends on the $Q_p$ of the planet but *not* $Q$ of the satellite, which affords a very significant simplification. Since the effective $Q_p$ of Saturn has been measured at various frequencies (Lainey et al. 2020), the equilibrium tidal heat production rate in satellites in MMRs (either current or ancient) can in principle be calculated (e.g. Meyer and Wisdom 2007). Thus, for instance, equilibrium tidal heating rates of 11 GW and 4.8 GW are obtained for the 2:1 present-day Enceladus-Dione resonance and a putative earlier 3:2 Mimas-Enceladus resonance, assuming a Saturn $Q_p$ of 1800 as suggested by astrometry (Lainey et al. 2020). A higher Saturn $Q_p$ would yield a correspondingly lower heat flow. Note that a tidal equilibrium for Enceladus requires it to be quite dissipative (with $k_{2E}/Q_E \simeq 0.01$ or even larger), likely due to the presence of an internal ocean, making solid-body dissipation estimates (e.g. Murray and Dermott 1999) invalid. Extensive further discussion of this issue is provided in Sect. 4.2 and Nimmo et al. (2023, this collection).

# 4 Important Resonance Passages in Saturn's Moon System

## 4.1 Mimas-Tethys 4:2 MMR

Mimas and Tethys are locked in a 4:2 inclination-type MMR, associated with the resonant argument $2\lambda_M - 4\lambda_\Theta + \Omega_M + \Omega_\Theta$. This argument slowly librates around 0° with a large
amplitude of 95° and a period of 70 years. Such a resonance is expected to raise the incli-
nations of the two involved bodies over time. However, since Tethys is about 16 times more
massive than Mimas, only the inclination of Mimas is significantly affected. Allan (1969)
derived the time evolution of these two inclinations since the trapping and has shown that,
while the current inclination of Mimas (1.6°) could be due to the excitation by this MMR,
the one of Tethys (1.1°) cannot, and was probably nearly as high as it is today when the
capture happened. Ćuk et al. (2016) propose an explanation for this inclination, involving
a past 5:3 MMR between Dione and Rhea, followed immediately by a secular resonance
between Tethys and Dione.

The trapping into this current resonance and its stability require that the orbits of Mi-
mas and Tethys converge. Since the trapping, assuming no other perturbations, the libration
amplitude has decreased from 180° at the onset of resonance to the currently observed 95°,
while the inclination of Mimas increased from the initial value to 1.6°. Using a linear model,
Sinclair (1972) found that when reaching the separatrix delimiting the resonance, the satel-
lites had only a 4% chance to be trapped. This low probability of the capture is due to the
large current libration amplitude of the resonant argument, which implies an initial inclina-
tion of 0.4° for Mimas. Champenois and Vienne (1999a,b) solved this paradox by showing
that a non-linear model, which is more accurate, allows for secondary resonances to affect
the libration amplitude in the main resonance. Allowing for much lower resonant libration
amplitudes post-capture leads to much higher capture probabilities. These probabilities may
almost reach 100% for a null initial free inclination of Mimas, which means that when Mi-
mas and Tethys encountered the resonance during their convergent migration, the trapping
into the currently observed configuration was the only plausible outcome.

Before the current 4:2 inclination resonance was reached, Mimas and Tethys had to cross
the sub-resonance with the argument $2\lambda_M - 4\lambda_\Theta + 2\Omega_M$ ("$i_M^2$ resonance"). In the simplest
case in which the pre-resonance inclination of Mimas was 0.4° the probability of capture
into the $i_M^2$ resonance was only 7% (Luan 2014), explaining why the capture did not happen.
However, scenarios in which the initial inclination of Mimas was small and the libration am-
plitude of the current ("$i_M i_\Theta$") resonance was subsequently changed (such as Champenois
and Vienne 1999b) would also make the capture into the $i_M^2$ resonance highly likely. There-
fore any hypothesis on the origin of the current 4:2 $i_M i_\Theta$ Mimas-Tethys resonance would
also need to account for the prior crossing of the $i_M^2$ resonance not resulting in capture.

None of these analytical and numerical studies actually constrain the dissipation rate inside Saturn, which is modelled as a scale temporal factor. The capture process and the evolution into resonance are assumed to be slow enough, i.e. adiabatic, to not be affected by the dissipation timescale. In assuming a Love number $k_{2,p} = 0.341$ for Saturn and that the dissipation function $Q_S$ is constant over the whole orbital evolution, Champenois and Vienne (1999a) estimate that the capture happened some 17,000 $Q_S$ years ago. This implies that the resonance was assembled only about 30 Myr ago using the $Q_S$ value from Lainey et al. (2012). Note that in equilibrium tidal theory (Eq. (1)) orbital evolution timescales of Mimas and Tethys are relatively similar. On the other hand, if the tidal response of Saturn depends on the frequency (Fuller et al. 2016), the estimated age of the resonance could greatly differ from the value given above.

Recent results, which imply non-equilibrium tides in Saturn, complicate the issue of the
history and age of this resonance. Capture into resonance requires convergent evolution
of Mimas and Tethys, which contradicts the end-member case in which all moon orbits
are expanding with constant semi-major axis ratios (Lainey et al. 2020). While the directly
measured tidal accelerations of Mimas and Tethys still have large uncertainties (Lainey et al.
2020; Jacobson 2022), the relatively fast observed evolution makes it unlikely that the age
of the Mimas-Tethys resonance is greater than a few tens of Myr. Better determination of the moons' tidal accelerations is necessary to put firmer limits on the age of this resonance.

## 4.2 Enceladus-Dione 2:1 MMR

Enceladus and Dione are currently in a 2:1 MMR with an argument $2\lambda_D - \lambda_E - \omega_E$ that acts to increase the eccentricity of Enceladus over time. The libration amplitude of this resonant argument is below one degree (Murray and Dermott 1999), implying long-term damping of the libration amplitude. Given that the eccentricity of Enceladus is modest ($e_E = 0.0047$), it was always considered likely that the tides within Enceladus act to decrease this eccentricity (Sect. 3). This eccentricity damping releases heat within Enceladus, affecting its geophysics and geology. As the eccentricity of Enceladus is modified over time by satellite tides, it is not possible to use the current orbits of the moons to determine the age of this resonance.

After the *Voyager* mission, it was known that parts of Enceladus's surface were very young, but current tidal heating was considered insufficient for resurfacing due to low eccentricity and then-prevailing estimates of the rate of Enceladus's orbital evolution (Squyres et al. 1983). As shown by Meyer and Wisdom (2007), if we assume that the Dione-Enceladus resonance is in a long-term equilibrium, the resonance produces $(18,000/Q_S) \times 1.1$ GW of heat within Enceladus. $Q_S = 18,000$ is the smallest value for which Mimas stays outside Saturn's Roche limit over the age of the Solar System (assuming frequency-independent equilibrium tides; Eq. (1)). As it was widely accepted before the *Cassini* mission that Saturn's major moons were primordial, long-term average tidal heating of Enceladus was thought to be restricted to 1 GW or less.

*Cassini* observations have shown that the heat flux from the south pole of Enceladus is 10-15 GW (Porco et al. 2006; Howett et al. 2011, 2024, this collection), an order of magnitude larger than predicted. Early solutions to this discrepancy suggested non-equilibrium scenarios (e.g. O'Neill and Nimmo 2010), until Lainey et al. (2012) found that fast tidal evolution of Saturn's moons suggested by astrometric observations can naturally explain the tidal heating in Enceladus being in equilibrium. Equilibrium heating of 15 GW implies $Q_E/k_{2E} = 100$ for Enceladus (Sect. 3), giving us an eccentricity-damping timescale (in absence of excitation) of 0.5 Myr. Therefore Enceladus would need only a few Myr to settle into a dynamical steady state, and an assumption of equilibrium does not preclude geophysical evolution on longer timescales (Nimmo et al. 2023, this collection).

More recent findings of non-equilibrium, likely resonant, tides dominating the evolution of the Saturnian system (Lainey et al. 2020) complicate this picture. Just like with the Mimas-Tethys resonance (previous subsection), the existence of the Enceladus-Dione resonance requires that the orbits of those two moons converge. Furthermore, unlike in the case of Mimas-Tethys resonance which only requires past net convergent evolution, tidal heating within the Enceladus-Dione resonance requires ongoing orbital convergence. This is consistent with the small libration amplitude which would be growing if the two moons' orbits were diverging, making the resonance evolve toward dissolution.

Convergent evolution is not compatible with the end-member orbital evolution model in which all moons are locked to resonant modes that all evolve on the same timescale (Lainey et al. 2020), or with a model in which the modes are divergent (Fuller et al. 2016). However, a situation in which Enceladus is in a resonance lock, but Dione is not, can result in a MMR. Using the approach of Meyer and Wisdom (2007) and adding the assumption about Enceladus being much less massive than Dione, we can roughly estimate the (equilibrium) tidal heating of Enceladus in the resonance as:

$$H \approx \frac{Gm_S m_D}{2t_a a_D} = \frac{5 \times 10^{28} \text{ J}}{t_a}, \qquad (9)$$
where $t_a = (a_D/\dot{a}_D)_{MMR}$ is the timescale for the evolution of Dione's orbit forced by resonance with Enceladus. In the case of equilibrium tides, $t_a$ would be the timescale of the two moons' orbital convergence multiplied by the moons' mass ratio $m_D/m_E$, so $t_a \approx 100$ Gyr if $Q_S = 1700$, yielding $H \approx 16$ GW. If Enceladus is evolving through a resonant lock, then $t_a$ is simply the timescale of the resonant mode's evolution, as Dione would be forced to evolve at the same rate. If we assume $t_a = 9$ Gyr, as proposed by Lainey et al. (2020) as being a typical value for orbital evolution due to resonance locking in the Saturnian system, $H \approx 180$ GW, more than an order of magnitude in excess of the observed value (Howett et al. 2011).

On the other hand, if the resonance lock evolves uniformly in a reference frame rotating with Saturn and the satellites' orbits evolve divergently (Fuller et al. 2016), steady-state heating of Enceladus will be much lower. The Saturn-frame resonance lock would produce $a/\dot{a} \approx 200$ Gyr for Enceladus, based on Titan's evolution timescale of 11 Gyr (Lainey et al. 2020) and $a/\dot{a} \propto n/(\Omega_p - n)$, where $\Omega_p$ is the planet's rotation rate (Fuller et al. 2016). This rate of evolution gives us a steady-state heat flow of about 8 GW for Enceladus, relatively close to the measured value, but almost certainly less than the true heat production rate. However, if we use Rhea's observed evolution ($a/\dot{a} = 6$ Gyr; Lainey et al. 2020) to estimate Enceladus's tidal heating in the same model, we get $t_a \approx 27$ Gyr and $H \approx 60$ GW. This last estimate is motivated by the main result of Jacobson (2022), who finds that Titan does not exhibit non-equilibrium tidal evolution but Rhea does. If we use the observational results for the secular acceleration of Rhea from Jacobson (2022), the timescale for Enceladus's evolution through resonant lock (assuming divergent modes) could be in the 30-40 Gyr range and the associated tidal heating in the 40-55 GW range. Note that the total tidal heating rate of Enceladus must exceed the measured value, as the thermal observations are not sensitive to distributed tidal heating outside the South Polar Terrain. Models of Enceladus's ice shell suggest a heat loss rate in the range 25-40 GW (Hemingway and Mittal 2019), so there may not be a discrepancy between the heating rates predicted here and the actual heat production rate (see Nimmo et al. 2023, this collection).

A separate constraint on the possible resonant-lock driven heating of Enceladus comes from the observational limits on the orbital deceleration of Dione. The minimum possible effective $k_{2,p}/Q_p$ for Dione found by Lainey et al. (2020) is $0.6 \times 10^{-4}$, which implies a slowest evolution of about 30 Gyr. This means that 55 GW is the minimum tidal heating in the Enceladus-Dione system allowed by the astrometric observations, if all of Dione's orbital evolution is due to Enceladus. Jacobson (2022) has larger error bars on the observed evolution of Dione, and allows for minimum heating rates that are about twice smaller. Note that these “lower limits” apply only to pure evolution of Enceladus by resonance lock that includes zero tidal torque on Dione. Models in which equilibrium tides are significant are not constrained by these limits.

We conclude that Enceladus's observed heat flux (which is a lower bound on the actual heat production rate) may be consistent with the resonance lock model in which the resonant modes are divergent, but is well short of that expected from the simple picture in which all resonant modes evolve on the same timescale. A picture in which modes affecting the inner satellites evolve more slowly (Fuller et al. 2016) can explain the current Enceladus-Dione resonance, but requires Dione to have avoided direct resonance locking to an internal mode within Saturn. Further refinement of our understanding of Enceladus's tidal heating will require more precise measurement of the current orbital evolution of the Saturnian moons, as well as more complete constraints on the thermal flux from Enceladus (Howett et al. 2024, this collection). Further discussion of Enceladus's tidal heat budget and thermal evolution may be found in Nimmo et al. (2023).
Unlike in the case of the Mimas-Tethys resonance, less attention has been given to the formation of the Dione-Enceladus resonance than to its continuing maintenance. In order for Enceladus and Dione to reach the current sub-resonance of their 2:1 MMR through convergent orbital evolution, numerous other sub-resonances must be crossed first. Past work (Meyer and Wisdom 2008; Zhang and Nimmo 2009) has concentrated on the sub-resonances that can be modelled assuming two bodies on planar orbits, including the second-order $4:2 e_E e_D$ mixed-eccentricity resonance. However, inclination-type resonances (equivalent to the one Mimas and Tethys are currently in) would be encountered first, and their effect on the (currently very low) inclination of Enceladus has yet to be fully explored. Additionally, using numerical integrations featuring all major moons, Ćuk and El Moutamid (2022) have recently shown that two-body resonances also contain additional three-body sub-resonances that can result in capture. Clearly more work is needed to understand the history of the Enceladus-Dione 2:1 resonance.

## 4.3 Past Resonances Between Mid-Sized Moons?

In addition to currently active resonances, researchers have long suspected that Saturn's moons may have experienced additional resonances in the past which have since been broken (Fig. 1). Here we will address several proposed past resonances among the inner moons (out to Rhea), while those involving Titan and the Sun (i.e. the evection resonance) are discussed in separate subsections. For an introduction to the topics of orbital migration, resonance capture and evolution within the resonance, we refer the reader to Sect. 2.

Among the inner moons, past resonances are often considered responsible for excited eccentricities and inclinations that cannot be accounted for with observed interactions. For example, while the inclination of Mimas is likely a product of its current $4:2$ MMR with Tethys, the eccentricity of Mimas is unaffected by this resonance and must predate it. Similarly, the large inclination of Tethys mostly precedes the resonance with Mimas, requiring an origin mechanism (Sect. 4.1). Here it is assumed that the moons formed on coplanar,

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Current location of nominal MMR for each satellite in terms of variation of its semi-major axis, assuming that the other moons do not move. Within a 10% variation in *a*, we reported all first-order resonances up to 5:4, all second-order resonances up to 5:3, the evection resonance for Rhea and the 5:1 resonance between Titan and Iapetus. Overlapping between resonances is due to currently active MMR</figcaption></figure>
circular orbits, as both eccentricity and inclination are damped by expected collisions with satellitesimals (Peale 1999). Apart from current orbital elements, the record of tectonic activity, such as Ithaca Chasma on Tethys, could indicate past tidal heating driven by high eccentricity (Chen and Nimmo 2008). Tethys's orbital eccentricity is presently very low, but past high eccentricity implies that Tethys experienced one or more orbital resonances during its history.

Reconstruction of past resonance crossings among the inner moons is highly dependent on the assumed strength and nature of tidal dissipation within Saturn. The potential past resonances that were studied in past work and are therefore discussed here are expected to have happened in the context of equilibrium tides. Equilibrium tidal theory predicts that all adjacent inner moon pairs are on converging paths, with the sole exception of Enceladus and Tethys, which were thought to be on divergent orbits (Murray and Dermott 1999). Frequency-dependent tidal dissipation within Saturn certainly changes this picture and may rule out some past resonances while introducing new ones. Here we will review past work on these resonances and follow with a brief discussion on how the latest results on Saturn's tidal response (Lainey et al. 2020; Jacobson 2022) change our assessment of these proposed resonant encounters.

**Mimas-Enceladus 3:2 Resonance** In the classical picture of equilibrium tides, one of the most recent resonance passages was that of the Mimas-Enceladus 3:2 MMR. This is expected to have been a convergent encounter, and has been proposed by Meyer and Wisdom (2008) as the source of Mimas's eccentricity. In particular, Meyer and Wisdom (2008) have suggested capture into $e_M$ and $e_M e_E$ (i.e. the 6:4 Mimas-Enceladus MMR) as possible sources of Mimas's eccentricity, while ruling out past capture into the $e_E$ sub-resonance as it would have been impossible to break. Tian and Nimmo (2020) later revisited this work using a wider range of tidal parameters, and found regions of tidal parameter space in which the hypothesis of Meyer and Wisdom (2008) is likely. Both of these studies relied on semi-analytical treatment of the resonant Hamiltonian, and assumed a planar system. The last issue is crucial, as the inclination of Enceladus is presently very low ($i_E \approx 0.01^\circ$), putting strong constraints on past resonance crossings involving Enceladus. El Moutamid et al. (2019) proposed that a long-lived Mimas-Enceladus 3:2 $e_E$ resonance was broken when the pair encountered respective 3:1 and 2:1 resonances with Dione. However, ongoing work by some of this work's authors (El Moutamid and Ćuk) using direct numerical integrations suggests that this triple resonance usually produces dynamical chaos, exciting the inclination of Enceladus beyond the observed value. In principle, a high inclination could have been damped by fluid tides in Enceladus's ocean (Sect. 3) but no quantitative study has been attempted. Clearly more work is needed to explore the dynamical viability of a past Mimas-Enceladus 3:2 MMR crossing, especially as the tidal evolution rates of these two moons (which may also not be constant over time) are hard to constrain by other means.

**Mimas-Dione 3:1 Resonance** The fact that Enceladus and Dione are currently in a 2:1 resonance implies that their orbital evolution (at least recently) has been convergent. Therefore, before this resonance was established Enceladus must have been interior to the resonant location. A clear implication is that if Mimas ever crossed the 3:2 MMR with Enceladus, it must have also crossed the 3:1 resonance with Dione (either at the same time or later). The reverse is not true: depending on the relative tidal evolution rates and the age of the system, Mimas could have in principle crossed the 3:1 resonance with Dione and not the 3:2 MMR with Enceladus. Given the large mass difference between Dione and Mimas, Mimas is likely to be affected by their 3:1 MMR much more than Dione. Meyer and Wisdom (2008) suggested the Mimas-Dione 3:1 MMR as another possible source of the eccentricity of Mimas.
They prefer capture into the $e_M e_D$ sub-resonance as it would also excite the eccentricity of Dione in line with constraints from Dione's later resonance passage with Enceladus. More recently Ćuk and El Moutamid (2022) simulated the Mimas-Dione 3:1 resonance and found that three-body resonances involving the eccentricity of Titan can be important. This three-body sub-resonance, not included in the semi-analytical model of Meyer and Wisdom (2008), usually leads to a short-lived capture that would excite the eccentricity of Mimas, consistent with the present state. Since Mimas's eccentricity is not being excited at present, survival of its eccentricity from an earlier time places limits on how dissipative it could be. At present, the 3:1 Mimas-Dione resonance appears to be a viable candidate for the source of Mimas's eccentricity, but more work is needed to confirm this.

**Tethys-Dione 3:2 Resonance** This resonance should have taken place relatively recently if Tethys and Dione were to follow equilibrium tidal theory in their orbital evolution. In the classical picture in which Saturn's tidal $Q_S = 18,000$ (Murray and Dermott 1999) and is constant, this resonance would have happened about 1 Gyr ago, while $Q_S \approx 1800$ proposed by Lainey et al. (2012) would put this resonance at 100 Myr ago. This resonance was proposed by Chen and Nimmo (2008) as the source of tidal heating of Tethys, and this hypothesis was explored in more detail by Zhang and Nimmo (2012) using numerical integrations. Zhang and Nimmo (2012) found that the Tethys-Dione 3:2 resonance was impossible to break once established, and proposed a large impact on Tethys (that formed the Odysseus basin) as a mechanism of breaking this resonance. However, Zhang and Nimmo (2012) assumed a planar system, therefore ignoring all constraints from the satellites' inclinations, and also underestimating the amount of chaos in the system (as all inclination-type resonances were absent in their simulations). Ćuk et al. (2016) studied this resonance using direct integrations with non-zero inclinations, and found that the inclination of Tethys is excited up to several degrees, after which a chaotic phase excites both the eccentricities and inclinations of both moons. The chaotic interactions ultimately break the resonance, but leave Tethys and Dione with inclinations well in excess of the observed ones. Unless both bodies experienced efficient inclination damping, these results imply that the Tethys-Dione 3:2 MMR never happened and therefore that the present inner satellite system of Saturn must be younger than 100 Myr (Ćuk et al. 2016, assuming $Q_p = 1800$). However, the recent astrometric results implying that the tidal response of Saturn is frequency-dependent suggest that Saturn's response to Tethys is significantly weaker than to other moons (Lainey et al. 2020; Jacobson 2022). A slowly-migrating Tethys would have either converged with Dione much more slowly or not at all (observations are still too uncertain to distinguish between these possibilities), so this resonance clearly cannot be used to put an upper limit on the age of the system. In summary, a passage through the Tethys-Dione 3:2 resonance would have produced inclinations inconsistent with the present values (Ćuk et al. 2016, unless the inclinations of both Tethys and Dione were greatly damped by some process), and its past occurrence is also disfavored by recent observational results on the moon's orbital evolution (Lainey et al. 2020; Jacobson 2022), so this resonance is unlikely to have much bearing on the history of the system.

**Dione-Rhea 5:3 Resonance** In the equilibrium tides-based relative chronology, the Dione-Rhea 5:3 resonance happens soon after the Tethys-Dione 3:2 resonance. Ćuk et al. (2016) studied the former using numerical simulations and unexpectedly found that this resonance between Dione and Rhea is a good candidate for the origin of Tethys's inclination. They find that following a largely chaotic Dione-Rhea 5:3 MMR passage, Dione and Tethys encounter a secular resonance with the argument $(\omega_\Theta + \Omega_\Theta) - (\omega_D - \Omega_D)$ (see Sect. 2 for
definitions). The two expressions within parentheses designate very slow-moving angles, as the apsidal and nodal precession rates for each moon are approximately equal and in opposite directions. However, precession rates are perturbed in proximity of MMRs, and in this case Rhea's perturbations just outside the Dione-Rhea 5:3 MMR modify the sum of Dione's apsidal and nodal precession rates just enough to produce the secular resonance. In the secular resonance, the eccentricity and inclination of Dione are transferred to Tethys in equal amounts (i.e. $sin(i_\Theta) \approx e_\Theta$). This sequence of resonances (Dione-Rhea 5:3 MMR followed by Tethys-Dione secular resonance) can explain the high inclinations of Tethys ($i_\Theta = 1^\circ$ before resonance with Mimas) and Rhea ($i_R = 0.33^\circ$) despite the low inclination of Dione ($i_D = 0.02^\circ$, decreased by the secular resonance). These resonances also produce past large eccentricity of Tethys ($e_\Theta \approx 0.01$) and the current eccentricity of Dione, assuming reasonable tidal dissipation. A past excited eccentricity of Tethys could in principle explain a past heating episode on Tethys inferred from Ithaca Chasma (Chen and Nimmo 2008). The free eccentricity of Rhea, however, is currently just $e_R \approx 2 \times 10^{-4}$ requiring either strong dissipation or some other mechanism for damping.

The biggest challenge to the past passage through the Dione-Rhea 5:3 MMR is the currently observed tidal evolution of Rhea (Lainey et al. 2020), which is five times faster than equilibrium tides would predict and would make the orbits of Dione and Rhea divergent. Two obvious solutions are that either the inclinations of Tethys and Rhea require different explanations, or Rhea's fast evolution is a phenomenon that started relatively recently. More work is urgently needed, both on exploring alternative possible sources of Tethys's inclination, and on the nature and operation of tidal dissipation within Saturn.

## 4.4 Rhea and the Evection Resonance

The evvection resonance is a semi-secular resonance in which a moon's precession period is equal to the parent planet's orbital period (Touma and Wisdom 1998). While no moons in the solar system are currently in the evvection resonance, studies of the early orbital history of Earth's Moon have identified the evvection resonance as an important dynamical mechanism (Touma and Wisdom 1998; Čuk and Stewart 2012; Tian et al. 2017; Rufu and Canup 2020). Evvection is a “semi-secular” resonance in the sense that the (apparent) mean motion of the perturber (the Sun) is in resonance with the apsidal precession of the satellite's orbit. As the moons' orbital precession are overwhelmingly driven by Saturn's oblateness, the dynamics of the evvection resonance does not depend on the relative arrangement of moons, as MMR do. The evvection resonance affects any moon at a certain distance from Saturn, which is given by:

$$a_{eve} = R_S \left( \frac{3}{2} J_2 \frac{\Omega_0}{n_S} \right)^{2/7}, \qquad (10)$$

where $J_2$ is Saturn's oblateness moment, $\Omega_0 = \sqrt{GM/R_S^3}$ is the orbital frequency at $a = R_S$ (Saturn's radius), and $n_S$ is Saturn's heliocentric orbital mean motion. Equation (10) gives $a_{eve} = 8.1R_S$, but Titan's perturbations shift the location of the evvection resonance to $a_{eve} = 8.2R_S$. This distance is somewhat smaller than the current orbital distance of Rhea (which has $a = 8.7R_S$).

Estimates of the past orbital evolution of Rhea based on extrapolation from observations (Lainey et al. 2020; Jacobson 2022) suggest that Rhea should have crossed the evvection resonance about 300-400 Myr ago. The evvection resonance, with the resonant argument $2\lambda_S - 2\sigma$ ($\lambda_S$ is Saturn's mean longitude, $\sigma$ Rhea's longitude of pericenter), is expected to
**Fig. 2** Simulations of Rhea's encounter with the solar ejection resonance using different parameters for the width of Saturn's resonant mode and Rhea's tidal dissipation. From top to bottom, the panels plot Rhea's semi-major axis, eccentricity and inclination. The purple and the green lines had $Q/k_2 = 1000$ for Rhea, and the light blue and the orange had $Q/k_2 = 100$ (i.e. likely value for Enceladus). The purple and the light blue lines had effective tidal $Q = Q_0 + (n-v/vσ)^2$ (where v is the mode's frequency and $σ = 10^{-6}$) due to the resonant mode, while the green and the orange lines used the same Q profile with $σ = 10^{-5}$. The evolution timescale is $a/ā = 6$ Gyr, as found for Rhea by Lainey et al. (2020). The dotted line on the bottom panel shows the present-day inclination of Rhea. This simulation included the Sun, Jupiter, Saturn, Titan, and Rhea. From work in preparation by MČ

<figure><img src="image_3.png" /></figure>

significantly excite Rhea's eccentricity. One possible outcome is that Rhea's satellite tides become strong enough to push the satellite inwards through the resonant mode, leaving it stranded interior to the mode. Another possibility would be Rhea evolving with the mode, while captured in the ejection resonance, producing a continuously increasing eccentricity. Yet another outcome would be a Rhea that continues evolving through the mode, but moves beyond the resonance, with Rhea's excited eccentricity slowly damping. Only the last of these possibilities would be consistent with the present orbital distance, small eccentricity ($\approx 0.001$) and observed evolution rate of Rhea.

Figure 2 shows four preliminary simulations done by author M.Ć. using a modified version of the numerical integrator SIMPL (Ćuk et al. 2016). One of the simulations shows a long-term resonance capture (green line), one “drops out” from the resonant mode (light blue line), while the other two simulations have Rhea staying with the mode but moving beyond the ejection resonance. However, the bottom panel that plots inclination, indicates that the ejection and associated resonances also significantly excite the *inclination* of Rhea, well in excess of the observed value ($i_R = 0.33°$).

The excitation of inclination appears to be chaotic and occurs concurrently with the ejection resonance. This is consistent with results obtained by Ćuk et al. (2016) assuming that Rhea migrated by equilibrium tides. The best candidate for the excitement of the inclination
is the resonance as $\lambda_S - \omega_S + \Omega - \Omega_S$ where $\Omega$ refers to the longitude of the ascending node and the subscript S refers to Saturn's heliocentric orbit. As the reference plane is the Laplace plane<fnref n="3" /> of Rhea's orbit, which is close to Saturn's equatorial plane, $\Omega_S$ is identical to the longitude of Saturn's equinox. Both $\Omega_S$ and $\omega_S$ are relatively slow-precessing (with periods of about 1.8 Myr and 50 kyr), so the resonance happens because $\dot{\lambda}_S = n_S \approx -\dot{\Omega}$. One can think of this resonant term as a combination of the Sun's “annual equation” (a term in classical lunar theory; Brouwer and Clemence 1961) and the Sun's main secular perturbation on the inclination of Rhea. In other words, the Sun induces perturbations in Rhea's inclination once per precession period, and solar perturbations vary in strength during one orbit of Saturn because of the planet's orbital eccentricity. When these two cycles are commensurable, this strong resonant perturbation occurs. This resonance is chaotic, as its strength depends on the eccentricity of Saturn's orbit, which varies in the $0.01 < e_S < 0.09$ range with a period of $5 \times 10^4$ yr due to Jupiter's perturbations. Jumps in inclination also affect the apsidal precession, making the main evection resonance also chaotic. Direct numerical integrations are clearly the only practical and accurate way of studying this problem.

In some cases, the orbital inclination of a moon can be damped by obliquity tides (Chyba et al. 1989), but this is unlikely to be applicable to Rhea. Large-scale inclination damping requires large forced obliquities, either through a Cassini State transition like the one experienced by the Moon (Chen and Nimmo 2016), or spin-orbit resonance, as proposed for Uranian moons (Čuk et al. 2020). The forced obliquity of Rhea is currently $\theta = 0.03^\circ$ (Chen et al. 2014), and the timescale for inclination damping is longer than that for eccentricity damping by a factor of $7(\sin\theta/\sin i)^2$ (Chyba et al. 1989), which in the case of Rhea is on the order of $10^3$. If Rhea is as dissipative as “Enceladus in equilibrium” with $Q/k_2 = 100$, the eccentricity damping timescale is on the order of 10 Myr (Fig. 2), which puts the timescale for inclination damping to many Gyr. Resonant tides excited in satellite oceans offer a different mechanism of inclination damping (Tyler 2008), but Chen et al. (2014) find that ocean obliquity tides would produce less dissipation in Rhea than classic obliquity tides that assume $Q/k_2 \simeq 100$. The above upper limits on dissipation may be overestimates, as Rhea currently does not appear to be in hydrostatic equilibrium (Tortora et al. 2016), and a solid tidal Love number $k_2 = 0.01$ (Chen et al. 2014) may be more applicable. Therefore, it appears that the damping of multi-degree inclination of Rhea within the last 0.5 Gyr is implausible.

It appears that Rhea should have crossed the evection resonance about 300-400 Myr ago and acquired a large inclination that could not have been subsequently damped. This dynamical history is clearly in conflict with Rhea's modest inclination ($i_R = 0.33^\circ$), and this disagreement suggests that one or more assumptions that produced this result are incorrect. While it cannot be completely ruled out that Rhea somehow avoided acquiring high inclination during resonance crossing, numerical results available so far suggest that it is highly unlikely that Rhea crossed the evection resonance, at least in the present dynamical environment. If the obliquity of Saturn is only 1-2 Gyr old (Saillenfest et al. 2021a), it is possible that Rhea did cross the evection resonance in the distant past when the out-of-plane perturbations by the Sun were much weaker. This would still require a separate solution for the timescale problem. One possible explanation is that Rhea became locked to a normal mode (Fuller et al. 2016) more recently than 400 Myr ago, with its orbital evolution before that being much slower. Another explanation would be that Rhea re-accreted outside the evection resonance distance, as a part of a dynamical cataclysm (e.g. Čuk et al. 2016). In order

<footnote n="3">³See Sect. 5.1 for a definition.</footnote>
to distinguish between these possibilities a better understanding of tidal dissipation within Saturn, both past and present, is likely needed.

The only firm conclusions we can draw from the proximity of Rhea to the ejection resonance is that we can exclude the global picture proposed by Lainey et al. (2020) in which the moons evolve by resonance lock practically since the system's formation. Unlike the tidal evolution of Titan where the results of Lainey et al. (2020) and Jacobson (2022) are in disagreement, both studies agree on the fast orbital evolution of Rhea, corresponding to a timescale of $(a/\dot{a}) \approx 6$ Gyr. Large-scale migration of Rhea, if it is a long-term steady state, assures relatively recent (300-400 Myr) crossing of the ejection resonance which would overexcite Rhea's inclination. Therefore, within the last few hundred Myr, Rhea must have experienced either a dramatic speeding up of its tidal evolution, or possibly a disruption and re-accretion.

## 4.5 Titan-Hyperion 4:3 MMR

Titan and Hyperion are currently in the 4:3 MMR, which keeps their orbital periods fixed in a 3:4 ratio. As Hyperion's orbit is relatively eccentric ($e_H = 0.104$), the resonance is important for protecting Hyperion from having close encounters and colliding with Titan. It is established that the eccentricity of Hyperion in the resonance grows as Titan migrates outward (Murray and Dermott 1999); in the context of equilibrium tides, the orbital evolution of Hyperion on its own has a negligible effect on the resonance. Therefore, the eccentricity of Hyperion could be used to constrain the past orbital evolution of Titan.

According to Eq. 8.242 in Murray and Dermott (1999), the increase in the eccentricity of Hyperion is given by:

$$ \frac{\dot{e}_H}{e_H} = \frac{1}{e_H^2} \frac{m_T}{m_S} n_H a_H \frac{F}{3g}, \qquad (11) $$

where $m, a, e$ and $n$ are masses, semi-major axes, eccentricities and mean motions as defined before, and subscripts $S, T$, and $H$ refer to Saturn, Titan and Hyperion. The other variables are defined as $F = 4\dot{n}_H - 3\dot{n}_T \approx -3\dot{n}_T$, and $g = 16Gm_T/a_H^2 + 9Gm_H/a_T^2 \approx 16Gm_T/a_H^2$. Expanding $F$ and $g$ into Eq. (11), and using $Gm_S/a_H^3 = n_H^2$, $n_H = (3/4)n_T$ and $\dot{n}/n = -(2/3)\dot{a}/a$, we get:

$$ 2e_H \dot{e}_H = \frac{d(e^2)}{dt} = \frac{1}{4} \frac{\dot{a}_T}{a_T}. \qquad (12) $$

Assuming that Hyperion started with a low eccentricity, and evolved to $e_H = 0.1$, Titan's semi-major axis must have evolved by about 4%, assuming a constant rate of migration. This is a very rough estimate, but it gives us the correct order of magnitude of how much Titan could have migrated since Hyperion was captured into the resonance.

Before 2012, when Saturn was estimated to have $Q_p > 18,000$ (Murray and Dermott 1999), Titan was thought to have migrated less than 1% of its orbital distance, and therefore the resonance was thought to have been established and largely evolved by non-tidal means, such as gas drag in protosatellite nebula (Peale 1999). Greenberg (1973) proposed an alternative view that the Titan-Hyperion resonance was evolved by tides alone, and implied a much faster global evolution rate of Saturn's moons.

When Lainey et al. (2012) presented evidence for the fast tidal evolution of Saturn's moons, Ćuk et al. (2013) noted that the new value of Saturn's tidal quality factor $Q_p \approx 1700$ would make the Titan-Hyperion resonance about as old as the Solar System. The more
recent direct measurement of Titan's orbital evolution puts the effective tidal $Q_p$ of Saturn at Titan's frequency as $Q_p \approx 120$, equivalent to an evolution timescale of $a/\dot{a} = 11$ Gyr (Lainey et al. 2020). Such a rapid expansion of Titan's orbit would place the origin of the Titan-Hyperion resonance at about 400-500 Myr ago.

The above calculation relies on two assumptions. One is that Hyperion's own rate of tidal evolution is negligible, and the other is that Hyperion did not damp its eccentricity. The first assumption about negligible tidal evolution is straightforward in the theory of equilibrium tides, in which a moon's tidal evolution rate is proportional to its mass. In the context of resonant mode locking (Fuller et al. 2016), the issue shifts from the rate of evolution to the question of whether Hyperion can lock to a mode within Saturn. Some quick calculations assuming $\dot{a}/a = 11$ Gyr (i.e. the observed rate of evolution of Titan) and using standard definition of tidal parameters (Murray and Dermott 1999) yield that Saturn would need to have a tidal quality factor $Q_p \approx 10^{-3}$ at Hyperion's frequency, which is nonphysical. Therefore, Hyperion is too small and too distant to lock to a resonant mode by itself.

The second issue of eccentricity damping is somewhat more difficult to model as Hyperion is in chaotic rotation (Wisdom et al. 1984). One way to estimate the amount of energy dissipation is to treat it like never-ending tidal spindown (Murray and Dermott 1999). This approach gives an eccentricity damping timescale of $10^{11}$ yr (assuming $Q_H = 10$ and $k_{2H} = 0.01$). Alternatively, we can estimate the energy dissipation within chaotic rotation as if Hyperion was a wobbling asteroid (Sharma et al. 2005), and this approach gives us an eccentricity damping timescale of $10^{12}$ yr. Unless another source of dissipation is identified, it appears unlikely that the eccentricity damping within Hyperion contributed significantly to the evolution of its resonance with Titan.

While the age of the Titan-Hyperion resonance is not the same as the age of Hyperion, there are strong reasons to think Hyperion cannot be much older than its 3:4 resonance with Titan. As Titan's and Hyperion's orbits converged, a primordial Hyperion would have likely been captured into other resonances that were crossed before the 4:3 resonance (notably the 3:2 and 7:5 resonances). Capture of Hyperion into those resonances is clearly inconsistent with Hyperion being on a relatively low-$e$, low-$i$ orbit before capture into the 4:3 resonance. While more work is needed to directly confirm the likelihood of capture into these outer resonances, we would expect at least the capture into the first-order 3:2 resonance to be robust. If Hyperion formed just interior to Titan's 3:2 resonance, Hyperion cannot be older than 1.5 Gyr, assuming Titan maintained the migration rate measured by Lainey et al. (2020). If capture into the mutual 7:5 resonance is also found to be a certain outcome, then Hyperion must be younger than about 1 Gyr.

Hyperion is a relatively small moon, but its late formation would very likely imply a wider cataclysm at that time. Hamilton (2013) has proposed that Titan was a late merger between multiple satellites, with Hyperion being an unaccreted fragment. Asphaug and Reufer (2013) also proposed a late formation mechanism for the Saturnian moons that involved a major impact on Titan. Confirming that Hyperion is not primordial would support the hypothesis of a late cataclysm involving Titan. On the other hand, if the slower migration of Titan implied by the orbital solutions of Jacobson (2022) is correct, Hyperion is likely to be primordial or almost primordial.

## 4.6 Past Titan-lapetus 5:1 MMR Crossing?

Iapetus is the third-largest moon of Saturn, as well as the major moon that is the most distant from the planet. Iapetus is notable for its albedo dichotomy (Buratti and Mosher 1995; Porco et al. 2005), oblate shape (Thomas et al. 2007; Castillo-Rogez et al. 2011), and equatorial
ridge (Levison et al. 2011; Dombard et al. 2012; Stickle and Roberts 2018), but here we will restrict ourselves to studying its orbital motion. Like other regular satellites, Iapetus has a relatively low orbital eccentricity ($e_I = 0.028$), but it also has a substantial orbital inclination ($i_I = 8°$ with respect to its Laplace plane<fnref n="4" />), the origin of which has been a long-standing problem (Ward 1981; Nesvorný et al. 2014). As the solar perturbations on Iapetus’s orbit are comparable to those arising from Saturn’s oblateness and the inner moons (chiefly Titan), the Laplace plane of Iapetus is significantly tilted to Saturn’s equator ($i_L = 14°$). As Iapetus’s orbit precesses around its Laplace plane, the instantaneous inclination of Iapetus to Saturn’s equator varies approximately over a $5°–21°$ range over Iapetus’s nodal precession period of about 3400 yr.

While in the classical picture (e.g. Murray and Dermott 1999) Iapetus does not take part in any resonances with other satellites, faster tidal evolution (Lainey et al. 2012, 2017) would make Titan and Iapetus cross their mutual 5:1 MMR in the past. This crossing should have happened about 500 Myr ago if we assume a uniform tidal quality factor $Q = 1500–2000$ for all satellites (Ćuk et al. 2013), or could have happened at a very different epoch if the tidal evolution of Saturn’s moons is driven by resonant modes inside the planet (Fuller et al. 2016).

Ćuk et al. (2013) have modeled the Titan-Iapetus 5:1 MMR crossing and found that the orbits of both bodies are chaotic during the crossing of the resonant region, which consists of numerous sub-resonances of the 5:1 resonance (they assumed $Q_p/k_{2,p} = 4000$ for Saturn). If Titan was as eccentric as it is now, the most likely outcome is eccentricity growth for Iapetus, followed by orbit crossing. If the eccentricity of Titan was low (about 0.005), Iapetus typically survives the resonance, acquiring an eccentricity of a few percent, consistent with its present orbit ($e_I = 0.03$). Ćuk et al. (2013) conclude that if the Titan-Hyperion 5:1 resonance was crossed in the past (which is unavoidable for $Q_p/k_{2,p} = 4000$), the current large eccentricity of Titan must postdate this resonance crossing. This suggestion is consistent with the later excitation of Titan’s eccentricity (and inclination) by an encounter with a large moon, now lost (Wisdom et al. 2022) (see Sect. 5).

Ćuk et al. (2013) also find that the inclination of Iapetus is affected only weakly by the 5:1 resonance crossing, with the typical change being only a degree or so, compared to the current $8°$. Therefore, this resonance cannot constrain the timing or the source of Iapetus’s inclination excitation. Interestingly, changes to the inclination of Titan can be comparable to the free inclination itself (which is $0.3°$), making it possible that Titan’s inclination was significantly modified by this resonance.

Polcarpe et al. (2018) have carried out numerical simulations of the resonance crossing using an N-body code as well as using averaged equations of motion. A large span of migration rates were explored for Titan and Iapetus was started on its local Laplace plane ($14°$ with respect to the equatorial plane) with a circular orbit. Polcarpe et al. (2018) find that the resonance crossing can trigger a chaotic evolution of the eccentricity and the inclination of Iapetus. The outcome of the resonance is highly dependent on the migration rate (or equivalently on $Q_p$). For a quality factor $Q_p$ of over around 2000, the chaotic evolution of Iapetus in the resonance leads in most cases to its ejection, while simulations with a quality factor between 100 and 2000 show a departure from the resonance with post-resonant eccentricities spanning from 0 up to 15%, and free inclinations capable of reaching $11°$. Usually high inclinations come with high eccentricities but some simulations (less than 1%) show elements compatible with Iapetus’ current orbit. Polcarpe et al. (2018) conclude a

<footnote n="4">⁴The instantaneous Laplace plane can be defined for every perturbed orbit as plane normal to the vector around which the orbit normal is precessing. See Sect. 5.1 for a formal definition.</footnote>
quality factor between 100 and 2000 at the frequency of Titan would bring Titan and Iapetus into a 5:1 resonance, which could in some cases perturb Iapetus' eccentricity and inclination to values observed today. Such rapid tidal migration would have avoided Iapetus' ejection around 40-800 Myr ago.

Čuk et al. (2018) find that Iapetus is currently in a secular resonance with an argument $\bar{\omega} - \Omega + \bar{\omega}_J - \Omega_{eq}$ librating around 180°, where $\bar{\omega}$ and $\bar{\omega}_J$ are the longitudes of pericenter of Iapetus and Jupiter, while $\Omega$ and $\Omega_{eq}$ are the longitudes of Iapetus's ascending node and Saturn's vernal equinox. The libration period is several Myr and the libration is likely to persist for several tens of Myr. Longer-term stability of this resonance is tied to the precession of Saturn's spin axis, and more definite predictions await better determinations of Saturn's precession rate (see Sect. 5). Most allowable solutions for Saturn's pole precession lead to eventual breaking of the Iapetus-$g_5$ secular resonance,<fnref n="5" /> but some solutions preserve the secular resonance for at least 100 Myr. The Iapetus-$g_5$ secular resonance was almost certainly established more recently than the proposed 5:1 MMR crossing between Titan and Iapetus (500-50 Myr ago, depending on the Titan's unknown tidal evolution rate). While Čuk et al. (2018) find cases when the secular resonance was established in the aftermath of this MMR (with the more rapidly evolving Titan offering promising results), they did not find a high-probability mechanism for establishing the secular resonance.

There are a number of profound mysteries regarding Iapetus's orbital history. There is no agreed-upon high-probability mechanism of producing Iapetus's very large inclination. The history of its secular resonance is also unclear. The eccentricity of Iapetus is probably the orbital parameter that is most closely coupled to the overall evolution of the system. Survival of Iapetus with a low eccentricity after the 5:1 Titan-Iapetus resonance is consistent with rapid evolution found by Lainey et al. (2020), with no restrictions on Titan's own eccentricity at the time. However, it is also consistent with a slower evolution of Titan (Jacobson 2022), but this requires that Titan acquired its eccentricity in the last 500 Myr (Čuk et al. 2013), requiring some kind of catastrophic event (cf. Asphaug and Reufer 2013; Hamilton 2013; Čuk et al. 2016). It is still not clear whether the low eccentricity of Iapetus can be consistent with a recent loss of a moon exterior to Titan (Wisdom et al. 2022), as such a moon would strongly interact with Iapetus before its disruption.

# 5 Obliquity of Saturn

In the previous sections, specific features of the satellite dynamics were reviewed in the context of strong tidal dissipation within Saturn. On a gigayear timescale, however, satellites cannot be considered in isolation; they are part of a vast coupled system that comprises the spin-axis dynamics of their host planet. The fast migration of Saturn's satellites, and of Titan in particular, has been shown to strongly affect the motion of Saturn's spin axis (see Saillenfest et al. 2021a,b; Saillenfest and Lari 2021; Wisdom et al. 2022). In this section, we review the main implications of these findings.

## 5.1 Basic Effect of Saturn's Satellites

From a satellite's perspective, the orientation of a planet's spin axis has direct dynamical consequences. During their formation in a disc, regular moons are naturally damped towards an equilibrium configuration with near-zero eccentricity. Their orbital inclination at

<footnote n="5">⁵The eigenmode with fundamental frequency $g_5$ dominates the apsidal precession of Jupiter, enabling us to use $\bar{\omega}_J$ as its proxy.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Equilibrium orbital inclination of a small satellite as a function of its semi-major axis. The orbital inclination of the satellite is measured with respect to the planet's equator. $\epsilon$ is the planet's obliquity</figcaption></figure>

equilibrium defines the “Laplace plane”, which results from a balance between the attraction of the equatorial bulge of the planet and the attraction of the star. For a massless satellite, the inclination of the Laplace plane measured from the planet’s equator is:

$$I_L = \frac{\pi}{2} + \frac{1}{2}\operatorname{atan2}[-\sin(2\epsilon), -r_M^5/a^5 - \cos(2\epsilon)] \quad (13)$$

(see Tremaine et al. 2009; Saillenfest and Lari 2021), where $\epsilon$ is the planet's obliquity. As shown in Fig. 3, the characteristic length $r_M$ is the distance at which the satellite's Laplace plane lies exactly halfway between the equator and the orbital plane of the planet (the index M stands for 'midpoint'). It can be written

$$r_M^5 = 2 \frac{M}{m_\odot} J_2 R_{eq}^2 a_P^3 (1 - e_P^2)^{3/2}. \quad (14)$$

Here, $m_\odot$ is the mass of the star, and $M$, $J_2$, and $R_{eq}$ are the mass, second zonal gravity coefficient, and equatorial radius of the planet. The orbit of the planet around the star is described by its semi-major axis $a_P$ and eccentricity $e_P$. The value of $r_M$ for Saturn is about 42 $R_{eq}$, so that the Laplace planes of most satellites are close to Saturn's equator. After their formation, regular moons keep oscillating about their local Laplace planes. Among the regular satellites of Saturn, only Iapetus presents today a substantial deviation from the exact equilibrium, about which it oscillates with an offset of 8° (see Sect. 4.6). From Eq. (13), we deduce that the spin-axis orientation of the planet directly sets the mean orbital inclination of its moons. Saturn's obliquity changes resulting from high dissipation within Saturn (see below) has therefore an indirect effect on Saturn's whole satellite system.

From a planet's perspective, a system of satellites has a long-term contribution to its spin-axis motion. In absence of satellites, the average motion of a planet's spin axis is merely due to the stellar torque applied on its equatorial bulge. Noting $\mathbf{s}$ as the unit spin-axis vector and $\mathbf{n}$ the unit orbit normal, this motion can be written

$$\dot{\mathbf{s}} = -\alpha(\mathbf{n} \cdot \mathbf{s})(\mathbf{n} \times \mathbf{s}) \quad (15)$$

(see e.g. Goldreich 1966; Tremaine 1991), where $\mathbf{n} \cdot \mathbf{s} = \cos \epsilon$. This equation produces a precession of the spin axis at a rate $\alpha \cos \epsilon$ around the instantaneous orbit normal $\mathbf{n}$ of the planet. The timescale of this precession is set by the “precession constant” $\alpha$, which can be expressed as

$$\alpha = \frac{3}{2} \frac{G m_{\odot}}{a_{\mathrm{P}}^{3}(1 - e_{\mathrm{P}}^{2})^{3/2}} \frac{J_{2}}{\lambda \omega}, \quad (16)$$

where $\lambda$ is the normalised polar moment of inertia of the planet, and $\omega$ is its spin rate. Satellites alter this motion in several ways: they contribute to the angular momentum of
the system, they apply a torque on the planet's equatorial bulge, and their orbits are themselves torqued by the stellar attraction. The relative magnitude of these effects depends on the satellite distances. The moon-planet coupling is most obvious for the Moon around the Earth, as the lunar mass is about 1% of the Earth's mass (for comparison, Titan's mass is only $2 \times 10^{-4}$ of Saturn's mass). As such, the first formulas that took into account the coupled motion of a planet's spin axis and the orbit of its moons were developed for the Earth-Moon system (see e.g. Goldreich 1966). From these formulas, we know that satellites enhance the mean spin-axis precession rate of their host planet in a way that is intimately linked to their local Laplace plane. Analytical expressions were obtained by Ward (1975) in the two extreme cases of close-in satellites (whose mean orbital plane is the planet's equator) and far-away satellites (whose mean orbital plane is the planet's own orbital plane; see Fig. 3). The first closed-formed expressions valid for any distance of the satellites were introduced by Tremaine (1991) and French et al. (1993). These expressions can be written as a modified precession constant $\alpha'$ for the planet, obtained by replacing $J_2$ and $\lambda$ in Eq. (16) by their effective values:

$$ 
\begin{aligned}
J'_2 &= J_2 + \frac{1}{2} \sum_k \frac{m_k}{M} \frac{a_k^2}{R_{\text{eq}}^2} \frac{\sin(2\varepsilon - 2I_k)}{\sin(2\varepsilon)}, \\
\lambda' &= \lambda + \sum_k \frac{m_k}{M} \frac{a_k^2}{R_{\text{eq}}^2} \frac{n_k}{\omega} \frac{\sin(\varepsilon - I_k)}{\sin\varepsilon}.
\end{aligned}
\qquad (17) 
$$

In this expression, the sum runs over all regular satellites (assumed to have circular orbits); $m_k$ is the mass of the $k$th satellite, $a_k$ is its semi-major axis, and $I_k$ is the inclination of its local Laplace plane with respect to the planet's equator. Equation (17) can be retrieved as a particular case of the self-consistent expressions of Boué and Laskar (2006): here, the motion of satellites is averaged over their precession around the Laplace plane. This approach is valid as long as they precess much faster than the planet's spin axis. This condition is well verified in practice for any distance of the satellites. As a result, the planet and its satellites rigidly precess as a whole (Goldreich 1965).

Equation (17) implies that the boost imparted by a moon on a planet's spin-axis precession rate strongly depends on the moon's distance. As an example, Fig. 4 shows the enhancement factor of Saturn's spin-axis precession rate as a function of Titan's distance. Despite the small mass of Titan as compared to Saturn's, the spin-axis precession rate of Saturn is currently multiplied by roughly a factor 4 due to the presence of Titan. This precession boost strongly depends on Titan's distance. The high dissipation within Saturn and

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Enhancement factor of Saturn's spin-axis precession rate due to its most massive satellite, Titan. The current semi-major axis of Titan is shown by a vertical line. The obliquity of Saturn is fixed to its current value $\varepsilon \approx 27^\circ$. The curve reaches its maximum at $(a/r_M)^3 \approx (\sqrt{\cos^2(2\varepsilon) + 24 - \cos(2\varepsilon)})/6$, see Saillenfest and Lari (2021)</figcaption></figure>
the resulting fast migration of Titan therefore have strong implications for Saturn's spin-axis dynamics.

## 5.2 Secular Spin-Orbit Resonances

Because of the existence of secular spin-orbit resonances, a variation in a planet's spin-axis precession rate can drastically change its spin-axis dynamics. To understand this mechanism, we must notice that the orbit pole **n** appearing in Eq. (15) is actually not fixed, but it precesses itself as a result of mutual planetary perturbations. As the orbital motion of the planet remains largely unaffected by its spin state, the vector **n** behaves in Eq. (15) as an autonomous forcing term. Building on the early works by Lagrange and Laplace, very accurate theories now exist for the long-term orbital motion of the Solar System planets. As an example, Table 1 gives the dominant terms of the secular motion of Saturn's orbit pole computed by Laskar (1990). The solution is written through the complex variable

$$ \zeta = \sin \frac{I}{2} \exp(i\Omega) = \sum_{k} S_{k} \exp[i\phi_{k}(t)], \quad (18) $$

in which the amplitudes $S_k$ are real constants and the sum runs over all terms with non-negligible amplitude. The angles $\phi_k$ evolves linearly over time as $\phi_k(t) = \nu_k t + \phi_k^{(0)}$, where $\nu_k$ is a fixed frequency and $\phi_k^{(0)}$ is the phase at $t = 0$.

Each term in the motion of the orbit pole **n** in Eq. (15) can generate a resonance for the motion of the spin axis **s**. These resonances are called “secular spin-orbit resonances”, because they involve the precession of the spin axis and the precession of the orbit. The geometry of the lowest-order resonances is given by Colombo’s top Hamiltonian (see Colombo 1966; Henrard and Murigande 1987; Saillenfest et al. 2019; Su and Lai 2022). As a function of the planet’s precession constant $\alpha$, the system near a resonance can have up to four equilibrium points which are called the “Cassini states” (Peale 1969). For a critical value of $\alpha$, a separatrix appears around the Cassini state 2; this separatrix divides the inside from the outside of the resonance. At this point, it becomes clear how even a set of small moons

<table><caption>**Table 1** First ten terms in the series decomposition of Saturn's inclination and longitude of ascending node, from Laskar (1990). The orbital elements of Saturn are measured in the J2000 ecliptic and equinox reference frame. Terms are sorted by decreasing amplitude and identified in terms of the fundamental frequencies $g_j$ and $s_j$ of the Solar System (using standard definitions, see e.g. Murray and Dermott 1999). The period associated with each frequency $\nu_k$ is given in the last column</caption><thead><tr><th>k</th><th>Identification</th><th>ν<sub>k</sub> ("yr<sup>-1</sup>)</th><th>S<sub>k</sub> × 10<sup>6</sup></th><th>φ<sub>k</sub><sup>(0)</sup> (°)</th><th>P<sub>k</sub> (kyr)</th></tr></thead><tbody><tr><td>1</td><td>s<sub>5</sub></td><td>0.00000</td><td>13,774</td><td>107.59</td><td></td></tr><tr><td>2</td><td>s<sub>6</sub></td><td>-26.33023</td><td>7850</td><td>127.29</td><td>49</td></tr><tr><td>3</td><td>s<sub>8</sub></td><td>-0.69189</td><td>560</td><td>23.96</td><td>1873</td></tr><tr><td>4</td><td>s<sub>7</sub></td><td>-3.00557</td><td>391</td><td>140.33</td><td>431</td></tr><tr><td>5</td><td>g<sub>5</sub> - g<sub>6</sub> + s<sub>7</sub></td><td>-26.97744</td><td>59</td><td>43.05</td><td>48</td></tr><tr><td>6</td><td>2g<sub>6</sub> - s<sub>6</sub></td><td>82.77163</td><td>34</td><td>128.95</td><td>16</td></tr><tr><td>7</td><td>g<sub>5</sub> + g<sub>6</sub> - s<sub>6</sub></td><td>58.80017</td><td>20</td><td>212.90</td><td>22</td></tr><tr><td>8</td><td>2g<sub>5</sub> - s<sub>6</sub></td><td>34.82788</td><td>16</td><td>294.12</td><td>37</td></tr><tr><td>9</td><td>s<sub>1</sub></td><td>-5.61755</td><td>14</td><td>168.70</td><td>231</td></tr><tr><td>10</td><td>s<sub>4</sub></td><td>-17.74818</td><td>13</td><td>123.28</td><td>73</td></tr></tbody></table>
can dramatically affect a planet's spin dynamics: because of their effective enhancement of the precession constant $\alpha'$, the presence of moons may result in the planet being inside instead of outside a given resonance. On top of that, the moons' tidal migration may slowly drive the value of $\alpha'$ past the critical value and produce a resonance capture. In case of resonance capture, the spin-axis precession rate of the planet becomes roughly constant, as it remains close to a given frequency $v_k$ appearing in the orbital precession (see Table 1). In other words, $\alpha' \cos\varepsilon \approx -v_k$. As the moons continue to migrate and produce large changes in the value of $\alpha'$ (see Fig. 4), the planet's obliquity $\varepsilon$ is forced to evolve so as to maintain the equality. This phenomenon is likely to happen for Jupiter in the future and produce an increase in its obliquity (Lari et al. 2020; Saillenfest et al. 2020).

## 5.3 Saturn Is Close to a Resonance...

By taking into the account the effect of Saturn's satellites, Ward and Hamilton (2004) and Hamilton and Ward (2004) have shown that Saturn is today very close to a strong secular spin-orbit resonance with $s_8$ (third term in Table 1). It was immediately recognised that this striking match and Saturn's large obliquity ($\varepsilon \approx 27^\circ$) are unlikely to be coincidental. Indeed, due to the process of gas accretion, giant gaseous planets are expected to form with a near-zero obliquity ($\varepsilon \approx 0^\circ$); some mechanism should therefore have tilted Saturn after its formation. As the strong resonance with $s_8$ has necessarily affected Saturn's obliquity in the past, it offers a natural explanation for its large obliquity (see Ward and Hamilton 2004; Hamilton and Ward 2004). Yet, in order to produce a resonance capture and tilting, the two frequencies involved (i.e. Saturn's spin-axis precession rate and Neptune's precession mode $s_8$) must have evolved and crossed the 1 to 1 commensurability. The question is when this crossing has happened in the past, and under which circumstances.

It is known today that shortly after the planetary formation, the orbits of the giant planets have been reshaped, as planets were migrating through a swarm of planetesimals (Malho tra 1993). The last large-scale changes in the orbits of the Solar System planets probably happened during this stage, more than four gigayears ago, when planetary migration trig gered a phase of orbital instability (see e.g. Tsiganis et al. 2005; Nesvorný and Morbidelli 2012). The dominant orbital precession modes of the planets are a direct function of their semi-major axes (see e.g. Murray and Dermott 1999); for this reason, planetary migration and instability will have produced a variation in the frequency $s_8$. By virtue of Eq. (16), Saturn's spin-axis precession rate also varied during this period of time. Hence, planetary migration first appeared to astronomers as a natural explanation for the crossing of Saturn's secular spin-orbit resonance with $s_8$. More than that, the current obliquity of Saturn was viewed as a practical constraint for the planetary migration process itself (Boué et al. 2009; Vokrouhlický and Nesvorný 2015; Brasser and Lee 2015).

However, the high tidal dissipation inside Saturn measured by Lainey et al. (2020) has changed our understanding of Saturn's tilting process. The migration rate measured for Titan is so fast that Titan has most probably migrated over a distance of several radii of Saturn in the past (and possibly even more than that; see Lainey et al. 2024, this collection, but also the contrary opinion of Jacobson 2022). Due to the steepness of its enhancement factor on Saturn's precession constant (see Fig. 4), this migration for Titan must have greatly altered Saturn's precession rate over time. As a result, it appears very unlikely that Saturn crossed the $s_8$ resonance during late planetary migration, more than four gigayears ago, because its spin-axis precession motion was too slow at that time. Instead, Titan's orbital expansion offers a new explanation for Saturn's capture into resonance and tilting to its current $27^\circ$ obliquity (Saillenfest et al. 2021a; Wisdom et al. 2022).
<figure><img src="image_6.png" /><figcaption>**Fig. 5** The three types of possible obliquity evolution for Saturn (adapted from Saillenfest et al. 2021a). Each panel shows an example of simulation for a given value of $\lambda$ (see labels). The explored range of values lies in the horizontal green bar. The first and last point of each trajectory are shown by a green spot. First-order secular spin-orbit resonances are shown in pink: the large area on the right is the resonance with $s_8$ and the thin area on the left is the resonance with $g_7 - g_8 + s_7$ (not shown in Table 1; 19th term). The separatrices of the $s_8$ resonance are highlighted in blue. Panels a and b show resonant capture while panel c does not</figcaption></figure>

The actual dynamical path followed by Saturn depends on its current location with respect to the resonance. As described by Eq. (16) and (17), the precise precession rate of Saturn depends on several physical parameters. Each of these parameters are well known today from direct measurements, except Saturn's normalised polar moment of inertia $\lambda$. Indeed, the gravitational potential measured by spacecraft only provides differences between the moments of inertia ($J_2, J_4$, etc.), but no measure of the individual moments of inertia. In order to obtain the value of Saturn's polar moment of inertia, authors usually use semi-empirical models of Saturn's interior structure that are fitted to Saturn's gravitational moments (see e.g. Hubbard and Marley 1989; Nettelmann et al. 2013; Movshovitz et al. 2020). Such estimates, although formally accurate, are model-dependent and generally do not agree with each other. Moreover, it is not clear how to convert the moment of inertia derived from state-of-the-art interior models to that used in simplified models of rotation dynamics (which use the gyroscopic approximation, a purely rigid rotation, no internal dynamical processes, etc.); the small discrepancy between these two concepts of moment of inertia may lead to a slight shift of Saturn's precession constant in the neighbourhood of the $s_8$ secular spin-orbit resonance. For these reasons, previous authors did not rely on a given value of $\lambda$, but they explored a range of possible values spanned by the estimates obtained through interior modelling. Using this approach, Saillenfest et al. (2021a) showed that Saturn could have followed three possible types of trajectories as a result of Titan's fast migration (see Fig. 5). For a given range of $\lambda$ values, one obtains very small past obliquities for Saturn (panel a) and a resonance capture with 100% certainty without crossing the separatrix. For slightly lower values of $\lambda$, one obtains a recent resonance capture by crossing the separatrix (panel b) or a resonance crossing without capture (panel c). By conducting large-scale Monte Carlo experiments, Saillenfest et al. (2021b) investigated the likelihood for each pathway as a function of Saturn's initial obliquity. The result is shown in Fig. 6 for Titan's nominal migration rate measured by Lainey et al. (2020). Trajectories labelled “b” are by far the less likely, because separatrix crossings are probabilistic events and in this case they rarely lead to a capture (see Saillenfest et al. 2021b). Trajectories labelled “a” are about ten times more likely that those labelled “c”. This higher likelihood is mostly due to the precession phase of Saturn's spin axis, which points today roughly in the direction of the
**Fig. 6** Likelihood of reproducing Saturn's current state as a function of its past obliquity (adapted from Saillenfest et al. 2021b). Each point is made of 240 numerical simulations with random initial precession angle for Saturn, started 4 Gyrs ago and propagated until today. The colour scale measures the relative likelihood of obtaining a final obliquity and precession angle within 0.5° and 5° of Saturn's current values, respectively. A small grey dot means that no successful trajectory was found. The labels a, b, and c refer to the trajectory types illustrated in Fig. 5

<figure><img src="image_7.png" /></figure>

centre of the resonance.<fnref n="6" /> Hence, the tilting of Saturn from a small obliquity – and Saturn still being in resonance today – is by far the most likely scenario from a dynamical point of view.

In this situation, Saturn's obliquity would continue to increase in the future as a result of the still ongoing migration of Titan; it could even reach values as high as $\epsilon \gtrsim 75^\circ$ within 5 Gyrs, at which point Titan would be strongly destabilised (see Saillenfest and Lari 2021).

## 5.4 ...but Saturn Is Just Out of Sync

The dynamical pathway followed by Saturn depends on its current location with respect to the $s_8$ resonance. If Saturn is inside the resonance today, then it would have been captured from a small obliquity and it would continue to follow the resonance as Titan migrates (Fig. 5a). If instead Saturn is outside the resonance today, then it would have recently crossed the resonance without capture, from a higher obliquity (Fig. 5c). In order to determine whether Saturn is currently inside or outside the resonance, Wisdom et al. (2022) developed new models of Saturn's interior structure that were fitted to the latest estimates of Saturn's gravitational moments measured by the Cassini spacecraft. From these models, they deduced refined values of Saturn's polar moment of inertia $\lambda$ (estimated to be equal to $0.2182_{-0.0003}^{+0.0006}$) that they used in long-term numerical integrations of Saturn's spin axis.

According to their results, Saturn is not in resonance today, but just outside its edge, at the level of 1% on the value of $\lambda$. As for the ambiguity and model dependency inherent to this approach (see above), Wisdom et al. (2022) argue that the four models that they consider cover the entirety of possible solutions that one could obtain from realistic models of Saturn's interior. The more recent interior models of Mankovich et al. (2023), which include a fit to ring seismology, also suggest that Saturn is out of resonance at present (see

<footnote n="6">⁶The current value of the resonance angle is $\sigma \approx 30^\circ$, which is relatively close to zero (Ward and Hamilton 2004). If Saturn is not locked in resonance today (as in case c), then $\sigma$ would circulate between $0^\circ$ and $360^\circ$ and have no particular reason to be small today; it would rather be a coincidence. If instead Saturn is locked in resonance today (as in case a), then $\sigma$ would oscillate around $0^\circ$ and always keep a relatively small value.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Variant scenario proposed by Wisdom et al. (2022). The small black dots show the types of trajectories "a" and "c" depicted in Fig. 5. The two portions of pathway followed by Saturn are superimposed in green and blue, respectively. When Saturn's ancient satellite is lost, the curve instantly jumps from the green to the blue curve (dashed segment)</figcaption></figure>

below). On this basis, Saturn should definitely be considered out of the $s_8$ resonance today. Saturn should therefore follow a pathway similar to that in Fig. 5c.

In this case, we are back to the question of the origin of Saturn's large obliquity. As recalled by Saillenfest et al. (2021a), an early giant impact could be a possibility. Instead of a giant impact, Wisdom et al. (2022) propose to connect the two possible types of pathways ("a" and "c") by a timely reorganisation in Saturn's moon system. Indeed, if a moon of Saturn is suddenly ejected or engulfed into Saturn, then Saturn's effective precession constant $\alpha'$ would jump to a lower value (see Eq. (17)). If this moon is massive enough, this would instantly kick Saturn out of the resonance. Wisdom et al. (2022) show that a relatively small ancient moon would be enough to produce this effect, for instance with a mass of the order of Iapetus's. Then, the still ongoing migration of Titan would make Saturn cross the resonance again, but without capture. This scenario is illustrated in Fig. 7. Ejection of the moon is hypothesized to occur when its orbit is destabilized during a 3:1 MMR with Titan (Wisdom et al. 2022).

As Saturn is still very close to the resonance today, and Titan's migration drives Saturn farther and farther away as time goes by, the removal of Saturn's ancient moon must have occurred recently in the history of the Solar System. For a former moon with the mass of Iapetus (that would not kick Saturn very far off the resonance), Wisdom et al. (2022) estimate that this event should have happened about 100 Myrs ago. Such a recent event reminds us of the controversial young age of Saturn's ring (see Crida et al., this collection). Wisdom et al. (2022) thus suggest to link these two phenomena: the lost moon would have been ripped apart within Saturn's Roche limit, which would have supplied the ring material, similarly to what has been proposed by Canup (2010). A recent instability in Saturn's system may also explain some peculiarities in the moons' orbits, such as Titan's eccentricity or Iapetus's inclination (see Sect. 4.6). It is still unclear, however, how such a catastrophic event would relate to the existence of Saturn's mid-sized moons, which all have rather unexcited orbits today and are connected through a fragile web of MMR (see Sect. 4).

The scenario of Wisdom et al. (2022) has the merit of preserving Titan's migration as the explanation for Saturn's large obliquity (see Saillenfest et al. 2021a), while allowing Saturn to be out of the resonance today. However, the likelihood analysis shows that, if one remains somewhat sceptical about the exact value of Saturn's polar moment of inertia, Saturn's current dynamical state is about ten times more likely to result from Saturn still
being in resonance today (mostly due to the current small value of the resonance angle $\sigma$; see footnote 6 above). It is also more probable in a cosmogonic point of view: if Saturn had freshly crossed the resonance, then humanity would be witnessing today a transient event that could have happened any time over the history of the Solar System. That being said, Saturn's ring does seem exceptional in the Solar System, so one may need to accept some exceptional level of chance to account for its existence (see the discussion by Maryame El Moutamid 2022). Furthermore, the Saturn system seems to have possessed a population of planetocentric impactors which are not observed elsewhere (Ferguson et al. 2022) and suggest that something unusual happened.

The debate is revived by the recent article by Jacobson (2022). His self-consistent adjustment of Saturn's spin axis and satellite dynamics to all available data yields a lower migration rate for Titan than previously measured (see Lainey et al. 2024, this collection). If we extrapolate the results of Jacobson (2022) into the past, we find that Titan would have migrated just the right amount to explain Hyperion's eccentricity (Ćuk et al. 2013; see Sect. 4.5) and Saturn's tilt (Saellenfest et al. 2021b; see Sect. 5.3) over the age of the Solar System. Moreover, Jacobson (2022) obtains a value for Saturn's moment of inertia<fnref n="7" /> that would put Saturn near the centre of the most probable region in Fig. 6; this would imply that Saturn's obliquity has evolved as in Fig. 5a, rather than the scenario proposed by Wisdom et al. (2022) and displayed in Fig. 7. If so, then the origin of Saturn's rings, Titan's eccentricity and the population of planetocentric impactors all remain to be explained. In the next few years, we can expect different teams of astronomers to revisit the studies of Wisdom et al. (2022) and Jacobson (2022); the future will tell which scenario will reach scientific consensus. Similarly, future work will hopefully yield agreement on the correct value for Saturn's polar moment of inertia $\lambda$.

A recent step in this direction has been made by Mankovich et al. (2023). The authors fitted interior models not only to the measured gravity harmonics of Saturn, but also on the frequencies of Saturn's internal oscillation modes imprinted in the rings. This approach allows the deep rotation profile of Saturn to be probed. From their models of internal structure, Mankovich et al. (2023) computed the expected value of Saturn's moment of inertia: they note that their results are consistent with the upper end of estimates of Wisdom et al. (2022), and approximately 0.5% smaller than the value required for the $s_8$ resonance to be active today. This supports the notion that Saturn is today outside the resonance, though close to its edge. Mankovich et al. (2023) also stress that large residuals in the frequencies of several modes (mostly for Saturn's deep structure, but also for outer layers) reflect the fact that current models for Saturn's interior are still not complete. According to them, a parallel can be drawn between this mismatch and the persistent difficulty encountered in fitting Jupiter and Saturn interior models to gravity while retaining the supersolar envelope heavy element enrichments indicated by spectroscopic measurements. These discussions reveal the richness and complexity of planetary interior modelling, and new avenues to better match the available observational constraints. Mankovich et al. (2023) advocate giving more freedom to the composition, temperature, and rotation profiles used as working interior models for their adjustments. It will be interesting to know whether such future refined interior models could change estimates for Saturn's moment of inertia by a non-negligible amount, and whether the strongly non-rigid rotation of Saturn that they reveal – still never taken into account in simulations of Saturn's long-term spin-axis dynamics – can alter in some way the precession rate of Saturn's spin axis.

<footnote n="7">⁷When converted to the same normalising radius $R_{eq}$ and spin rate $\omega$ as those used in Sect. 5.3, the values obtained by Jacobson (2022) for Saturn's normalised moment of inertia are $\lambda = 0.2283$ and $0.2313$ (see Eq. (16)), with a $1\sigma$ uncertainty of about 1%. The two values correspond to two different sets of observations.</footnote>
# 6 Interior Evolution of Moons

The interior evolution of moons is driven by heating processes that can affect the state (e.g., solid or liquid), location (e.g., homogeneous or segregated), and composition of their materials. These materials appear to be primarily ices (chiefly water ice) and, as inferred from meteorites and interplanetary dust particles, silicate rock and refractory organic material. The proportion of ice is inferred from the moons' bulk densities (985 kg m⁻³ for Tethys to 1879 kg m⁻³ for Titan) to be at least 50 vol.%, with some uncertainty arising from porosity (for smaller moons such as Mimas) and compression (for Titan). Heating predominantly arises from the long-term decay of radioactive isotopes sourced in silicate rock and from tidal dissipation (Eq. (2)) with additional contributions, likely limited in time, from the radioactivity of short-lived isotopes, release of potential energy during accretion or separation of materials by density (differentiation), and chemical reactions between materials (e.g., rock hydration). Radiogenic heating drives larger, more silicate-rich moons, with a lower ratio of surface area (heat loss) to volume (heat production), to have warmer, more evolved interiors. Tidal heating drives closer-in, more eccentric or inclined moons to be warmer and thus more evolved (Eq. (2)). Tidal dissipation also couples the interior and orbital evolution of a moon and, because of gravitational interactions (e.g., resonances) within the moon system, to those of the other moons, Saturn, and its rings as well. While interior and orbital evolution aspects can and have been studied separately for specific satellites, or a subset of the moon system, system-level studies of the coupled interior-orbital evolution are challenging to carry out owing to the diversity of timescales involved (from orbital periods of days to radioactive decay timescales of billions of years) and the larger space of uncertain or unknown parameters.

A number of observational constraints provide a starting point. In addition to the dynamical constraints on ongoing and past orbital evolution described in previous sections, these include geological observations (see Rhoden et al. 2024, this collection, for an in-depth discussion). The moons' ice-rich surfaces and (assuming hydrostatic equilibrium) moments of inertia indicate some degree of ice-rock separation, but this separation is incomplete inside at least Titan (Iess et al. 2010) and perhaps Enceladus (Hemingway and Mittal 2019) and Rhea (Tortora et al. 2016). The data at other moons are too uncertain to tell (Tajeddine et al. 2014; Beuthe et al. 2016). Enceladus (see Rhoden et al. 2024, this collection) and, perhaps, Titan (Nixon et al. 2018) show evidence of ongoing endogenic geological activity, and Tethys and Dione show evidence of past activity linked to past high heat fluxes (Chen and Nimmo 2008; Hammond et al. 2013; White et al. 2017). The sizes of the moons increase with increasing semi-major axis from Mimas to Titan, but their bulk densities (rock content, including heat-producing radioisotopes) do not.

An especially noteworthy observation, the 'Mimas paradox', is the puzzling contrast in geological activity between inactive Mimas and (cryo)volcanic Enceladus. From Eq. (2), the dissipation of solid tides raised by Saturn, which dominates radiogenic heating at these small, close-in moons, is roughly 40 times stronger at Mimas than at Enceladus, assuming a similar interior structure and material response (i.e., $k_2/Q$). Why, then, is Enceladus so active and Mimas so quiescent? Proposed solutions to this paradox involve different interior material responses between Mimas and Enceladus, i.e., different loci in the positive feedback loop in which increasing interior temperatures lead to decreasing material viscosities, which lead to increasing tidal dissipation and more heating. Malamud and Prialnik (2013) proposed that Enceladus was nudged into this positive feedback loop by a spike of chemical heating during the hydration of the core, possible at Enceladus but not Mimas due to the former's higher radiogenic heat supply (Enceladus is denser and thus has a higher rock
fraction than Mimas). Czechowski and Witek (2015) suggested that this nudging could have also been caused by higher short-lived radiogenic heating at Enceladus. Neveu and Rhoden (2019) invoked more radical differences: a young age for Mimas and earlier encounter of MMR by Enceladus which would not be today in equilibrium heating. Rhoden and Walker (2022) suggested that Mimas may have an ocean just like Enceladus, which could be compatible with observations of its libration (Tajeddine et al. 2014). Whether this suggestion is compatible with the lack of obvious indications of extensive relaxation of Mimas's large impact basins remains to be established. An ocean-bearing Mimas would also have to have avoided damping of Mimas's eccentricity (Sect. 4.3).

The long-term evolution of the moons depends on how and when they formed. This is largely an open question for the Saturn system. While Titan and Iapetus are likely primordial, the age of Rhea and closer-in moons is unclear (Rhoden et al. 2024, this collection). They too could have formed concurrently with Saturn at the dawn of the Solar System (Mosqueira et al. 2010, and references therein), possibly coexisting and interacting with a long-lived massive ring (Nakajima et al. 2019), or from the reaccretion of moon fragments from prior collisions with Sun-orbiting impactors (Movshovitz et al. 2015), or from accretion at the outer edge of Saturn's rings followed by tidally driven orbital expansion (Charnoz et al. 2010). Each formation scenario is compatible only with a subset of possible orbital expansion histories (Sect. 1): roughly, the end-members are accelerating expansion for primordial or old moons and decelerating expansion for young moons. Astrometric observations of fast orbital expansion (Lainey et al. 2020) currently lean toward the former (but see Jacobson 2022). If moons formed from rings, rock (silicates and organics) and ice may have been separated inside moons at the time of formation: the lesser propensity of rock to undergo tidal disruption would have led first to the formation of a rock-rich core, onto which ice would then have accreted as tidally driven orbital expansion led to decreased tidal forcing (Charnoz et al. 2011). This appears to be the case for Mimas, whose predicted fast orbital expansion, dominated by interactions with Saturn's nearby rings, suggests a young age if Mimas postdates the rings (Rhoden et al. 2024, this collection), and whose libration necessitates some degree of internal differentiation (Tajeddine et al. 2014). However, there is no independent evidence that Mimas postdates the rings. The higher propensity of close-in moons (Enceladus, Tethys, Dione) to encounter MMR implies a greater potential for interior evolution via episodes of high heat flow or collisions, especially if the moons are old. This is reflected in the broader diversity of geological features at their surfaces relative to Mimas, which may be younger, and Rhea, which is farther out.

In contrast, the interior evolution of Rhea, Titan, and Iapetus appears dominated by radiogenic heating and thus less tumultuous (see Rhoden et al. 2024, this collection). Dating events in the moons' history, such as past periods of geological activity, can constrain or test hypotheses on their internal evolution (Zahnle et al. 2003); progress on crater-based dating accounting for planetocentric impactor populations (Ferguson et al. 2020, 2022) and future in situ measurements could decrease the currently considerable uncertainties on surface ages.

# 7 Conclusions

The long-term evolution of the Saturnian system has probably never been more contentious. The proposed views of the system have changed over the years, from a primordial one driven by relatively slow tides (Murray and Dermott 1999; Meyer and Wisdom 2008), through a fast-evolving system (Lainey et al. 2012) that may not be primordial (Ćuk et al. 2016),
to a resonance-lock driven one that could either be primordial (Fuller et al. 2016; Lainey et al. 2020) or recently reset (Wisdom et al. 2022). Our models of the system's long-term evolution are closely related to the observational results on the Saturnian rings, satellite astrometry and crater statistics which are covered in other articles of this collection. Some of the past controversies have been resolved, as we now know that the rings have a relatively low mass (Iess et al. 2019), even if their age is still in question (Crida et al. 2019), and the impact craters on the satellites are at least partially caused by planetocentric debris (Ferguson et al. 2020, 2022). New controversies include the current migration rate of Titan (Lainey et al. 2020; Jacobson 2022) and whether Saturn's spin precession is still in a resonance with Neptune's orbit (Saillenfest et al. 2021a; Wisdom et al. 2022).

One firm data point is that the observed tidal evolution of Rhea is well constrained (Lainey et al. 2017, 2020; Jacobson 2022) and is incompatible with equilibrium tides, indicating that highly frequency-dependent (probably resonant) tides must be operating in the Saturnian system (Lainey et al. 2024, this collection). While this is technically just an instantaneous migration rate, the existence of resonant tides greatly influences our understanding of the past evolution of the system. In this article we reviewed past work on the long-term dynamical evolution of Saturn's satellites, and attempted to show how various resonances and other dynamical features constrain the current theories about the system's age and tidal evolution. Our findings can be summarized as follows:

1. The currently observed Mimas-Tethys 4:2 and Enceladus-Dione 2:1 resonances indicate that in each of these two pairs, the satellites are on converging orbits. This is incompatible with all satellites being on diverging (Fuller et al. 2016) or parallelly-evolving orbits (Lainey et al. 2020). If the system is primarily driven by resonance locking, these resonances indicate that not all satellites are currently locked to Saturn's internal modes.

2. The current heat flux of Enceladus (if it is in equilibrium) is consistent both with strong equilibrium tides (Lainey et al. 2012) or diverging resonant modes (Fuller et al. 2016), but is likely significantly smaller than the expected tidal heating if Enceladus is locked to one of the parallelly-evolving modes (as illustrated in Fig. 3 of Lainey et al. 2020).

3. The excitation of Tethys's 1° inclination is an important constraint on the evolution of the system before the current resonances. It is possible that Tethys's inclination is a consequence of Dione-Rhea 5:3 resonance crossing followed by Tethys-Dione secular resonance (Ćuk et al. 2016), but that scenario would require the current fast evolution of Rhea to be a very recent phenomenon. If Rhea's current rate of evolution is a long-term value, then Dione and Rhea never crossed their 5:3 resonance.

4. Another “fossil” feature of the inner system is the excited eccentricity of Mimas, which could have been excited by the past crossing of the Mimas-Dione 3:1 resonance (Ćuk and El Moutamid 2022).

5. In every modern theory of its tidal evolution, Rhea should have crossed the semi-secular ejection resonance with the Sun which should have excited Rhea's eccentricity and inclination well beyond the observed values. The damping of Rhea's inclination would be difficult to achieve, possibly indicating that Rhea may not be primordial but that it could have reaccreted close to its present orbit.

6. The current Titan-Hyperion 4:3 resonance has been evolved by tides and, depending on which estimate of Titan's tidal evolution rate we use, may be either 500 Myr old (using Lainey et al. 2020) or 4 Gyr old (using Jacobson 2022). In the former case Hyperion is not primordial and must be a product of a cataclysm of some kind within the last 1-1.5 Gyr.

7. Titan and Iapetus must have crossed their mutual 5:1 resonance at some point in the past. Once again, different estimates of Titan's migration rate place this crossing at either 50 Myr ago (Lainey et al. 2020) or 500 Myr ago (Jacobson 2022). The latter case of slower
evolution through the resonance leads to more dynamical excitation, so in this case Titan must have acquired its eccentricity after the 5:1 resonance with Iapetus, again indicating some kind of possibly catastrophic, recent event.

8. Titan's relatively high eccentricity and mechanically weak interior have long been thought to be problematic, indicating that the eccentricity had to be generated relative recently or it would have been damped. Mechanisms for this excitation usually involve some sort of catastrophe in which at least some moons are disrupted (Asphaug and Reufer 2013; Hamilton 2013; Ćuk et al. 2016; Wisdom et al. 2022)

9. A resonance between the precession of Saturn's spin axis and the orbit of Neptune almost certainly generated Saturn's obliquity (Ward and Hamilton 2004; Hamilton and Ward 2004; Saillenfest et al. 2021a,b). The substantial orbital expansion of Titan offers a natural origin for this resonance capture. Wisdom et al. (2022) find that Saturn is currently not in this spin-orbit resonance, but this claim depends on the value of Saturn's moment of inertia derived from interior models, which are not all in agreement. However, if the resonance is broken, this would strongly indicate a recent cataclysm in the outer part of the satellite system (beyond Titan), although the viability of the specific scenario of Wisdom et al. (2022) and the connection to the formation of the rings may need further confirmation.

The above points demonstrate that there are no easy answers to the questions posed by the Saturnian system. It appears that there is no simple unified explanation for the orbital evolution of all moons (past and observed), and while some moons may not be primordial, the age and extent of any disruption events is highly uncertain.

Fortunately, there are two observable quantities that may be constrained enough in the near future to help us narrow down the range of the system's possible histories. One is the precession rate of Saturn's spin axis. Current measured values lean toward Saturn being in the spin-orbit resonance with Neptune (Jacobson 2022), though non-resonant values are also permitted. If Saturn is determined to be in the spin-orbit resonance, that would remove the need for (but not fully rule out) the proposed cataclysm in the outer part of the system (which is more important for Saturn's precession than the inner satellites). The second observable quantity is the current rate of tidal evolution of Titan. At the time of writing there are two contradictory results in the literature (Lainey et al. 2020; Jacobson 2022) that are largely based on the same data. If a consensus could be reached on this issue, we would have better understanding of the physics of tidal evolution in the system. Titan is the largest moon and is relatively far from Saturn. As such it is most likely to be ancient and most likely to be locked to a resonant mode within Saturn. Confirmation of a fast-evolving Titan would demonstrate the reality of resonant locking; a slower migrating Titan could indicate that fast migration driven by resonant modes may not be permanent.

Apart from the two above mentioned observable quantities, progress may come from sources of data other than the satellite dynamics. Better understanding of Saturn's interior (Mankovich and Fuller 2021), the rings (Iess et al. 2019) and satellite surfaces (Ferguson et al. 2022, also see Rhoden et al. 2024, this collection) could all help us understand the history of the system. Furthermore, more theoretical work on the moons' past orbital evolution, while not assumption-free, could still produce valuable new insights on the possible histories of the system, especially if combined with constraints from other lines of evidence.

**Acknowledgements** The authors thank the International Space Science Institute and Dr. Valery Lainey for organizing the workshop "New Vision of the Saturnian System in the Context of a Highly Dissipative Saturn" that inspired this collection.

**Author Contribution** All authors have contributed to this review article.
**Funding** Open access funding provided by Università di Pisa within the CRUI-CARE Agreement. MĆ is supported by NASA Solar System Workings award 80NSSC22K0979. ME and MĆ are supported by NASA-SSW award 80NSSC23K0681. FN is supported by NASA-SSW award 80NSSC21K1802. MN is supported by NASA-CDAP award 80NSSC21K0531 and by the CRESST II agreement between NASA GSFC and Univ. Maryland, College Park (award number 80GSFC17M0002). AR is supported by NASA Cassini Data Analysis Program award 80NSSC21K0535.

**Data Availability** Not applicable.

**Code Availability** Not applicable.

**Declarations**

**Competing Interests** The authors have no relevant financial or non-financial interests to disclose.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**References**

Allan RR (1969) Evolution of Mimas-Tethys commensurability. Astron J 74:497. https://doi.org/10.1086/110827

Asphaug E, Reufer A (2013) Late origin of the Saturn system. Icarus 223(1):544–565. https://doi.org/10.1016/j.icarus.2012.12.009

Beuthe M, Rivoldini A, Trinh A (2016) Enceladus’s and dione’s floating ice shells supported by minimum stress isostasy. Geophys Res Lett 43(19):10–088

Boué G, Laskar J (2006) Precession of a planet with a satellite. Icarus 185(2):312–330

Boué G, Laskar J, Kuchynka P (2009) Speed limit on Neptune migration imposed by Saturn tilting. Astrophys J 702(1):L19–L22

Brasser R, Lee MH (2015) Tilting Saturn without tilting Jupiter: constraints on giant planet migration. Astron J 150(5):157

Brouwer D, Clemence GM (1961) Methods of celestial mechanics. Academic Press, San Diego

Buratti BJ, Mosher JA (1995) The dark side of Iapetus: additional evidence for an exogenous origin. Icarus 115(2):219–227. https://doi.org/10.1006/icar.1995.1093

Canup RM (2010) Origin of Saturn’s rings and inner moons by mass removal from a lost Titan-sized satellite. Nature 468(7326):943–946

Castillo-Rogez JC, Efroimsky M, Lainey V (2011) The tidal history of Iapetus: spin dynamics in the light of a refined dissipation model. J Geophys Res, Planets 116(E9):E09008. https://doi.org/10.1029/2010JE003664

Champenois S, Vienne A (1999a) Chaos and secondary resonances in the Mimas-Tethys system. Celest Mech Dyn Astron 74(2):111–146

Champenois S, Vienne A (1999b) The role of secondary resonances in the evolution of the Mimas-Tethys system. Icarus 140(1):106–121. https://doi.org/10.1006/icar.1999.6115

Charnoz S, Salmon J, Crida A (2010) The recent formation of Saturn’s moonlets from viscous spreading of the main rings. Nature 465(7299):752–754

Charnoz S, Crida A, Castillo-Rogez JC, Lainey V, Dones L, Karatekin Ö, Tobie G, Mathis S, Le Poncin-Lafitte C, Salmon J (2011) Accretion of Saturn’s mid-sized moons during the viscous spreading of young massive rings: solving the paradox of silicate-poor rings versus silicate-rich moons. Icarus 216(2):535–550

Chen EMA, Nimmo F (2008) Implications from Ithaca Chasma for the thermal and orbital history of Tethys. Geophys Res Lett 35(19):L19203. https://doi.org/10.1029/2008GL035402
Chen EMA, Nimmo F (2016) Tidal dissipation in the lunar magma ocean and its effect on the early evolution of the Earth-Moon system. Icarus 275:132-142. https://doi.org/10.1016/j.icarus.2016.04.012

Chen EMA, Nimmo F, Glatzmaier GA (2014) Tidal heating in icy satellite oceans. Icarus 229:11-30. https://doi.org/10.1016/j.icarus.2013.10.024

Chyba CF, Jankowski DG, Nicholson PD (1989) Tidal evolution in the Neptune-Triton system. Astron Astrophys 219(1-2):L23-L26

Colombo G (1966) Cassini's second and third laws. Astron J 71:891

Crida A, Charnoz S, Hsu H-W, Dones L (2019) Are Saturn's rings actually young? Nat Astron 3:967-970. https://doi.org/10.1038/s41550-019-0876-y

Ćuk M, El Moutamid M (2022) Three-body resonances in the Saturnian system. Astrophys J Lett 926(2):L18. https://doi.org/10.3847/2041-8213/ac501c

Ćuk M, Stewart ST (2012) Making the Moon from a fast-spinning Earth: a giant impact followed by resonant despinning. Science 338(6110):1047. https://doi.org/10.1126/science.1225542

Ćuk M, Dones L, Nesvorný D (2013) Titan-hyperion resonance and the tidal Q of Saturn. arXiv e-prints. arXiv:1311.6780

Ćuk M, Dones L, Nesvorný D (2016) Dynamical evidence for a late formation of Saturn's Moons. Astrophys J 820(2):97. https://doi.org/10.3847/0004-637X/820/2/97

Ćuk M, Dones L, Nesvorný D, Walsh KJ (2018) Secular resonance between iapetus and the giant planets. Mon Not R Astron Soc 481(4):5411-5421

Ćuk M, El Moutamid M, Tiscareno MS (2020) Dynamical history of the Uranian system. Planet Sci J 1(1):22. https://doi.org/10.3847/PSJ/ab9748

Czechowski L, Witek P (2015) Comparison of early evolutions of Mimas and Enceladus. Acta Geophys 63(3):900-921

Dermott SF, Malhotra R, Murray CD (1988) Dynamics of the Uranian and Saturnian satellite systems: a chaotic route to melting Miranda? Icarus 76(2):295-334

Dombard AJ, Cheng AF, McKinnon WB, Kay JP (2012) Delayed formation of the equatorial ridge on Iapetus from a subsatellite created in a giant impact. J Geophys Res, Planets 117(E3):E03002. https://doi.org/10.1029/2011JE004010

El Moutamid M (2022) How Saturn got its tilt and its rings. Science 377(6612):1264-1265

El Moutamid M, Ćuk M, Tiscareno M (2019) The orbital connection between Mimas and Enceladus. In: AAS/Division of Dynamical Astronomy meeting. AAS/Division of Dynamical Astronomy meeting, vol 51, p 102.04

Ferguson SN, Rhoden AR, Kirchoff MR (2020) Small impact crater populations on Saturn's moon Tethys and implications for source impactors in the system. J Geophys Res, Planets 125(9):e2020JE006400

Ferguson SN, Rhoden AR, Kirchoff MR (2022) Regional impact crater mapping and analysis on Saturn's moon dione and the relation to source impactors. J Geophys Res, Planets 127(6):e2022JE007204

French RG, Nicholson PD, Cooke ML, Elliot JL, Matthews K, Perkovic O, Tollestrup E, Harvey P, Chanover NJ, Clark MA, Dunham EW, Forrest W, Harrington J, Pipher J, Brahic A, Grenier I, Roques F, Arndt M (1993) Geometry of the Saturn system from the 3 July 1989 occultation of 28 Sgr and Voyager observations. Icarus 103(2):163-214.

Fuller J, Luan J, Quataert E (2016) Resonance locking as the source of rapid tidal migration in the Jupiter and Saturn moon systems. Mon Not R Astron Soc 458:3867-3879

Fuller J, Guillot T, Mathis S, Murray C (2024) Tidal dissipation in giant planets. Space Sci Rev 220

Goldreich P (1965) Inclination of satellite orbits about an oblate precessing planet. Astron J 70:5

Goldreich P (1966) History of the lunar orbit. Rev Geophys Space Phys 4:411-439

Goldreich P, Gold T (1963) On the eccentricity of satellite orbits in the solar system. Mon Not R Astron Soc 126(3):257-268

Greenberg R (1973) Evolution of satellite resonances by tidal dissipation. Astron J 78:338. https://doi.org/10.1086/111423

Hamilton DP (2013) Late major merger at Saturn: consequences for Titan and iapetus. In: AAS/Division for Planetary Sciences meeting abstracts. AAS/Division for Planetary Sciences meeting abstracts, vol 45, p 302.01

Hamilton DP, Ward WR (2004) Tilting Saturn. II. Numer Model Astron J 128(5):2510-2517

Hammond NP, Phillips CB, Nimmo F, Kattenhorn SA (2013) Flexure on dione: investigating subsurface structure and thermal history. Icarus 223(1):418-422

Hemingway DJ, Mittal T (2019) Enceladus's ice shell structure as a window on internal heat production. Icarus 332:111-131

Henrard J, Murigande C (1987) Colombo's top. Celest Mech 40:345-366

Howett CJA, Spencer JR, Pearl J, Segura M (2011) High heat flow from Enceladus' south polar region measured using 10-600 cm⁻¹ Cassini/CIRS data. J Geophys Res, Planets 116(E3):E03003. https://doi.org/10.1029/2010JE003718
Howett C, Miles G, Quick LC (2024) Constraining the mass and energy of dissipation systems. Space Sci Rev 220

Hubbard WB, Marley MS (1989) Optimized Jupiter, Saturn, and Uranus interior models. Icarus 78(1):102–118

Hussmann H, Spohn T (2004) Thermal-orbital evolution of Io and Europa. Icarus 171(2):391–410

Iess L, Rappaport NJ, Jacobson RA, Racioppa P, Stevenson DJ, Tortora P, Armstrong JW, Asmar SW (2010) Gravity field, shape, and moment of inertia of Titan. Science 327(5971):1367–1369

Iess L, Militzer B, Kaspi Y, Nicholson P, Durante D, Racioppa P, Anabtawi A, Galanti E, Hubbard W, Mariani MJ, Tortora P, Wahl S, Zannoni M (2019) Measurement and implications of Saturn's gravity field and ring mass. Science 364(6445):aat2965. https://doi.org/10.1126/science.aat2965

Jacobson RobertA (2022) The orbits of the main Saturnian satellites, the Saturnian system gravity field, and the orientation of Saturn's pole. Astron J 164(5):199

Lainey V, Karatekin Ö, Desmars J, Charnoz S, Arlot J-E, Emelyanov N, Le Poncin-Lafitte C, Mathis S, Remus F, Tobie G, Zahn J-P (2012) Strong tidal dissipation in Saturn and constraints on Enceladus' thermal state from astrometry. Astrophys J 752(1):14. https://doi.org/10.1088/0004-637X/752/1/14

Lainey V, Jacobson RA, Tajeddine R, Cooper NJ, Murray C, Robert V, Tobie G, Guillot T, Mathis S, Remus F, Desmars J, Arlot J-E, De Cuyper J-P, Dehant V, Pascu D, Thuillot W, Le Poncin-Lafitte C, Zahn J-P (2017) New constraints on Saturn's interior from Cassini astrometric data. Icarus 281:286–296. https://doi.org/10.1016/j.icarus.2016.07.014

Lainey V, Gomez Casajus L, Fuller J, Zannoni M, Tortora P, Cooper N, Murray C, Modenini D, Park R, Robert V, Zhang Q (2020) Resonance locking in giant planets indicated by the rapid orbital expansion of Titan. Nat Astron 4:1053–1058

Lainey V, Zannoni M, Iess L, Robert V, Fuller J, Mathis S, Guillot T (2024) Quantification of tides in giant planets from observations. Submitted to Space Sci Rev

Lari G, Saillenfest M, Fenucci M (2020) Long-term evolution of the Galilean satellites: the capture of Callisto into resonance. Astron Astrophys 639:A40

Laskar J (1990) The chaotic motion of the solar system: a numerical estimate of the size of the chaotic zones. Icarus 88(2):266–291

Levison HF, Walsh KJ, Barr AC, Dones L (2011) Ridge formation and de-spinning of Iapetus via an impact-generated satellite. Icarus 214(2):773–778. https://doi.org/10.1016/j.icarus.2011.05.031

Luan J (2014) Special mean motion resonance pairs: Mimas-Tethys and Titan-Hyperion. arXiv e-prints. https://doi.org/10.48550/arXiv.1410.2648. arXiv:1410.2648

Malamud U, Prialnik D (2013) Modeling serpentinization: applied to the early evolution of Enceladus and Mimas. Icarus 225(1):763–774

Malhotra R (1993) The origin of Pluto's peculiar orbit. Nature 365(6449):819–821

Mankovich CR, Fuller J (2021) A diffuse core in Saturn revealed by ring seismology. Nat Astron 5:1103–1109. https://doi.org/10.1038/s41550-021-01448-3

Mankovich CR, Dewberry JW, Fuller J (2023) Saturn's seismic rotation revisited. Planet Sci J 4(4):59

Meyer J, Wisdom J (2007) Tidal heating in Enceladus. Icarus 188(2):535–539. https://doi.org/10.1016/j.icarus.2007.03.001

Meyer J, Wisdom J (2008) Tidal evolution of Mimas, Enceladus, and Dione. Icarus 193(1):213–223

Moore WB, Schubert G (2000) The tidal response of Europa. Icarus 147(1):317–319

Mosqueira I, Estrada P, Turrini D (2010) Planetesimals and satellitesimals: formation of the satellite systems. Space Sci Rev 153(1):431–446. https://doi.org/10.1007/s11214-009-9614-6

Movshovitz N, Nimmo F, Korycansky DG, Asphaug E, Owen JM (2015) Disruption and reaccretion of mid-sized moons during an outer solar system late heavy bombardment. Geophys Res Lett 42(2):256–263

Movshovitz N, Fortney JJ, Mankovich C, Thorngren D, Helled R (2020) Saturn's probable interior: an exploration of Saturn's potential interior density structures. Astrophys J 891(2):109

Murray CD, Dermott SF (1999) Solar system dynamics. Cambridge University Press, Cambridge

Nakajima A, Ida S, Kimura J, Brasser R (2019) Orbital evolution of Saturn's mid-sized moons and the tidal heating of Enceladus. Icarus 317:570–582

Nesvorný D, Morbidelli A (2012) Statistical study of the early solar system's instability with four, five, and six giant planets. Astron J 144:117

Nesvorný D, Vokrouhlický D, Deienno R, Walsh KJ (2014) Excitation of the orbital inclination of Iapetus during planetary encounters. Astron J 148(3):52. https://doi.org/10.1088/0004-6256/148/3/52

Nettelmann N, Püstow R, Redmer R (2013) Saturn layered structure and homogeneous evolution models with different EOSs. Icarus 225(1):548–557

Neveu M, Rhoden AR (2019) Evolution of Saturn's mid-sized moons. Nat Astron 3(6):543–552

Nimmo F, Neveu M, Howett C (2023) Origin and evolution of Enceladus's tidal dissipation. Space Sci Rev 219(7):57. https://doi.org/10.1007/s11214-023-01007-4
Nixon CA, Lorenz RD, Achterberg RK, Buch A, Coll P, Clark RN, Courtin R, Hayes A, Iess L, Johnson RE et al (2018) Titan's cold case files-outstanding questions after Cassini-Huygens. *Planet Space Sci* 155:50–72

Ojakangas GW, Stevenson DJ (1986) Episodic volcanism of tidally heated satellites with application to Io. *Icarus* 66(2):341–358

O'Neill C, Nimmo F (2010) The role of episodic overturn in generating the surface geology and heat flow on Enceladus. *Nat Geosci* 3(2):88–91. https://doi.org/10.1038/ngeo731

Peale SJ (1969) Generalized Cassini's laws. *Astron J* 74:483

Peale SJ (1999) Origin and evolution of the natural satellites. *Annu Rev Astron Astrophys* 37:533–602. https://doi.org/10.1146/annurev.astro.37.1.533

Polycarpe W, Saillenfest M, Lainey V, Vienne A, Noyelles B, Rambaux N (2018) Strong tidal energy dissipation in Saturn at Titan's frequency as an explanation for iapetus orbit. *Astron Astrophys* 619:A133

Porco CC, Baker E, Barbara J, Beurle K, Brahic A, Burns JA, Charnoz S, Cooper N, Dawson DD, Del Genio AD, Denk T, Dones L, Dyudina U, Evans MW, Giese B, Grazier K, Helfenstein P, Ingersoll AP, Jacobson RA, Johnson TV, McEwen A, Murray CD, Neukum G, Owen WM, Perry J, Roatsch T, Spitale J, Squyres S, Thomas PC, Tiscareno M, Turtle E, Vasavada AR, Veverka J, Wagner R, West R (2005) Cassini imaging science: initial results on Phoebe and Iapetus. *Science* 307(5713):1237–1242. https://doi.org/10.1126/science.1107981

Porco CC, Helfenstein P, Thomas PC, Ingersoll AP, Wisdom J, West R, Neukum G, Denk T, Wagner R, Roatsch T, Kieffer S, Turtle E, McEwen A, Johnson TV, Rathbun J, Veverka J, Wilson D, Perry J, Spitale J, Brahic A, Burns JA, Del Genio AD, Dones L, Murray CD, Squyres S (2006) Cassini observes the active south pole of Enceladus. *Science* 311(5766):1393–1401. https://doi.org/10.1126/science.1123013

Renaud JP, Henning WG (2018) Increased tidal dissipation using advanced rheological models: implications for io and tidally active exoplanets. *Astrophys J* 857(2):98

Renaud JP, Henning WG, Saxena P, Neveu M, Bagheri A, Mandell A, Hurford T (2021) Tidal dissipation in dual-body, highly eccentric, and nonsynchronously rotating systems: applications to Pluto–Charon and the exoplanet TRAPPIST-1e. *Planet Sci J* 2(1):4

Rhoden AR, Walker ME (2022) The case for an ocean-bearing Mimas from tidal heating analysis. *Icarus* 376:114872

Rhoden AR, Ferguson SN, Castillo-Rogez JC, Martin E, Bland M, Zannoni M, Rambaux N, Salmon J (2024) Geologic constraints on the formation and evolution of Saturn's mid-sized moons. Submitted to *Space Sci Rev*

Rovira-Navarro M, Katz RF, Liao Y, van der Wal W, Nimmo F (2022) The tides of Enceladus' porous core. *J Geophys Res, Planets* 127(5):e2021JE007117

Rufu R, Canup RM (2020) Tidal evolution of the ejection resonance/quasi-resonance and the angular momentum of the Earth-Moon system. *J Geophys Res, Planets* 125(8):e06312. https://doi.org/10.1029/2019JE006312

Saillenfest M, Lari G (2021) Future destabilisation of Titan as a result of Saturn's tilting. *Astron Astrophys* 654:A83

Saillenfest M, Laskar J, Boué G (2019) Secular spin-axis dynamics of exoplanets. *Astron Astrophys* 623:A4

Saillenfest M, Lari G, Courtot A (2020) The future large obliquity of Jupiter. *Astron Astrophys* 640:A11

Saillenfest M, Lari G, Boué G (2021a) The large obliquity of Saturn explained by the fast migration of Titan. *Nat Astron* 5:345–349

Saillenfest M, Lari G, Boué G, Courtot A (2021b) The past and future obliquity of Saturn as Titan migrates. *Astron Astrophys* 647:A92

Sharma I, Burns JA, Hui CY (2005) Nutational damping times in solids of revolution. *Mon Not R Astron Soc* 359(1):79–92. https://doi.org/10.1111/j.1365-2966.2005.08864.x

Sinclair AT (1972) On the origin of the commensurabilities amongst the satellites of Saturn. *Mon Not R Astron Soc* 160:169. https://doi.org/10.1093/mnras/160.2.169

Squyres SW, Reynolds RT, Cassen PM, Peale SJ (1983) The evolution of Enceladus. *Icarus* 53(2):319–331. https://doi.org/10.1016/0019-1035(83)90152-5

Stickle AM, Roberts JH (2018) Modeling an exogenic origin for the equatorial ridge on Iapetus. *Icarus* 307:197–206. https://doi.org/10.1016/j.icarus.2018.01.017

Su Y, Lai D (2022) Dynamics of Colombo's top: tidal dissipation and resonance capture, with applications to oblique super-earths, ultra-short-period planets and inspiraling hot Jupiters. *Mon Not R Astron Soc* 509(3):3301–3320

Tajeddine R, Rambaux N, Lainey V, Charnoz S, Richard A, Rivoldini A, Noyelles B (2014) Constraints on Mimas' interior from Cassini ISS libration measurements. *Science* 346(6207):322–324

Thomas PC, Burns JA, Helfenstein P, Squyres S, Veverka J, Porco C, Turtle EP, McEwen A, Denk T, Giese B, Roatsch T, Johnson TV, Jacobson RA (2007) Shapes of the Saturnian icy satellites and their significance. *Icarus* 190(2):573–584. https://doi.org/10.1016/j.icarus.2007.03.012
Tian Z, Nimmo F (2020) Implications of second-order resonance for the thermal and orbital evolution of Mimas. Mon Not R Astron Soc 492(1):369–376. https://doi.org/10.1093/mnras/stz3427

Tian Z, Wisdom J, Elkins-Tanton L (2017) Coupled orbital-thermal evolution of the early Earth-Moon system with a fast-spinning Earth. Icarus 281:90–102. https://doi.org/10.1016/j.icarus.2016.08.030

Tortora P, Zannoni M, Hemingway D, Nimmo F, Jacobson RA, Iess L, Parisi M (2016) Rhea gravity field and interior modeling from Cassini data analysis. Icarus 264:264–273

Touma J, Wisdom J (1998) Resonances in the early evolution of the Earth-Moon system. Astron J 115(4):1653–1663. https://doi.org/10.1086/300312

Tremaine S (1991) On the origin of the obliquities of the outer planets. Icarus 89(1):85–92

Tremaine S, Touma J, Namouni F (2009) Satellite dynamics on the Laplace surface. Astron J 137(3):3706–3717

Tsiganis K, Gomes R, Morbidelli A, Levison HF (2005) Origin of the orbital architecture of the giant planets of the solar system. Nature 435(7041):459–461

Tyler RH (2008) Strong ocean tidal flow and heating on moons of the outer planets. Nature 456(7223):770–772. https://doi.org/10.1038/nature07571

Vokrouhlický D, Nesvorný D (2015) Tilting Jupiter (a bit) and Saturn (a lot) during planetary migration. Astrophys J 806(1):143

Ward WR (1975) Tidal friction and generalized Cassini's laws in the solar system. Astron J 80:64–70

Ward WR (1981) Orbital inclination of Iapetus and the rotation of the Laplacian plane. Icarus 46(1):97–107. https://doi.org/10.1016/0019-1035(81)90079-8

Ward WR, Hamilton DP (2004) Tilting Saturn. Astron J 128(5):2501–2509

White OL, Schenk PM, Bellagamba AW, Grimm AM, Dombard AJ, Bray VJ (2017) Impact crater relaxation on dione and Tethys and relation to past heat flow. Icarus 288:37–52

Wisdom J (1980) The resonance overlap criterion and the onset of stochastic behavior in the restricted three-body problem. Astron J 85:1122–1133. https://doi.org/10.1086/112778

Wisdom J (2008) Tidal dissipation at arbitrary eccentricity and obliquity. Icarus 193(2):637–640

Wisdom J, Peale SJ, Mignard F (1984) The chaotic rotation of Hyperion. Icarus 58(2):137–152. https://doi.org/10.1016/0019-1035(84)90032-0

Wisdom J, Dbouk R, Militzer B, Hubbard WB, Nimmo F, Downey BG, French RG (2022) Loss of a satellite could explain Saturn's obliquity and young rings. Science 377(6612):1285–1289

Zahnle K, Schenk P, Levison H, Dones L (2003) Cratering rates in the outer solar system. Icarus 163(2):263–289. https://doi.org/10.1016/S0019-1035(03)00048-4

Zhang K, Nimmo F (2009) Recent orbital evolution and the internal structures of Enceladus and Dione. Icarus 204(2):597–609. https://doi.org/10.1016/j.icarus.2009.07.007

Zhang K, Nimmo F (2012) Late-stage impacts and the orbital and thermal evolution of Tethys. Icarus 218(1):348–355. https://doi.org/10.1016/j.icarus.2011.12.013

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Authors and Affiliations

Matija Ćuk¹ ID • Maryame El Moutamid² ID • Giacomo Lari³ ID • Marc Neveu⁴,⁵ ID •

Francis Nimmo⁶ ID • Benoît Noyelles⁷ ID • Alyssa Rhoden⁸ ID • Melaine Saillenfest⁹ ID

<figure><img src="image_9.png" /></figure>

G. Lari
giacomo.lari@unipi.it

M. Ćuk
mcuk@seti.org

M. El Moutamid
maryame@astro.cornell.edu

M. Neveu
marc.f.neveu@nasa.gov

F. Nimmo
fnimmo@ucsc.edu

<figure><img src="image_10.png" /></figure>
B. Noyelles
benoit.noyelles@univ-fcomte.fr

A. Rhoden
alyssa@boulder.swri.edu

M. Saillenfest
melaine.saillenfest@obspm.fr

1. Carl Sagan Center, SETI Institute, 339 N Bernardo Ave, Mountain View, 94043, CA, USA

2. Cornell Center of Astronomy and Planetary Sciences, Cornell University, Space Sciences Building, Ithaca, 14850, NY, USA

3. Dipartimento di Matematica, Università di Pisa, Largo Bruno Pontecorvo 5, 56127, Pisa, Italy

4. Department of Astronomy, University of Maryland, 4296 Stadium Dr., College Park, 20742, MD, USA

5. Planetary Environments Laboratory, NASA Goddard Space Flight Center, 8800 Greenbelt Rd., Greenbelt, 20771, MD, USA

6. Department of Earth and Planetary Sciences, University of California Santa Cruz, 1156 High St, Santa Cruz, CA 95064, USA

7. Institut UTINAM UMR 6213 / CNRS, Univ. of Franche-Comté, OSU THETA, BP 1615, 25010, Besançon Cedex, France

8. Department of Space Studies, Southwest Research Institute – Boulder, 1050 Walnut St., Boulder, 80302, CO, USA

9. IMCCE, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Université, Université de Lille, 77 av. Denfert-Rochereau, 75014, Paris, France