DYNAMICAL EVIDENCE FOR A LATE FORMATION OF SATURN'S MOONS

One of the most important results of the Cassini mission to Saturn was the discovery of a high heat flow from the south polar region of Enceladus (Porco et al. 2006). This heat flow is estimated at 10–15 GW (Howett et al. 2011, 2014), which is by an order of magnitude in excess of theoretical predictions assuming a dynamically steady state (Meyer & Wisdom 2007). A dynamical equilibrium means that Enceladus maintains a constant orbital eccentricity, which corresponds to an equilibrium between excitation by the 2:1 orbital resonance with Dione and damping by tidal dissipation within Enceladus. When calculating this dynamical equilibrium, Meyer & Wisdom (2007) assumed that the strength of tidal dissipation within Saturn (which pushes Enceladus into the orbital resonance) can be quantified by a tidal quality factor Q = 18,000, as suggested by Sinclair (1983). Tidal Q = 18,000 corresponds to the strongest tidal response by Saturn that is still consistent with Mimas not evolving past its current location over 4.5 Gyr (Sinclair 1983; Murray & Dermott 1999). Any larger tidal evolution rates (corresponding to smaller tidal Q values) would place Mimas within Saturn's rings less than 4.5 Gyr ago.

There are other possible approaches to constraining the tidal response of Saturn. One is direct observation of the moons' tidal acceleration. Secular expansion of satellite orbits with orbital periods longer than the planet's spin period is predicted by tidal theories (Murray & Dermott 1999) but was previously thought to be too slow for detection. Recently, Lainey et al. (2012) reanalyzed a large number of historical observations of Saturn's moons and found tidal recession corresponding to Q = 1700 ± 500. More recent preliminary results using Cassini-based and ground-based observations of the satellites apparently support the idea of such fast tidal evolution of Saturn's moons (Lainey et al. 2015). Lainey et al. (2012) also point out that Q ≃ 1700 matches what is required for Enceladus to be in tidal equilibrium and still produce the observed heat flux. However, such fast tidal rates would require the inner moons like Mimas to have formed less than about 0.5 Gyr ago (as they would have otherwise migrated past their present positions, assuming that the Q of Saturn did not change over time). Theoretical models of Saturn's interior dissipation are also available, but they are not conclusive to this debate, as there is a basic disagreement about whether the dissipation happens in the core or the envelope (Ogilvie 2014; Remus et al. 2015; Storch & Lai 2015).

We can also constrain Saturn's tidal response by theoretical modeling of present and past orbital resonances between different moons. Low Q (i.e., a high rate of dissipation) of Saturn on the order of Q = 1700 would make the apparent highly evolved state of the Titan–Hyperion resonance a natural product of tidal evolution over the age of the solar system (Greenberg 1973a; Ćuk et al. 2013), whereas in the classical picture the resonance had to be primordial, as it would not have not changed significantly over 4.5 Gyr (Peale 1999). Here we address the two most important orbital resonances that have likely happened in the past: the 3:2 mean-motion resonance (MMR) between Tethys and Dione, and the 5:3 MMR between Dione and Rhea. These resonances should have happened close in time about 1.1 Gyr ago if Q = 18,000 (Murray & Dermott 1999; Meyer & Wisdom 2007), or about 100 Myr ago if Q = 1700. The 3:2 Tethys–Dione resonance has previously been studied using semianalytical methods, which found substantial excitation of the two moons' orbital eccentricities but ignored orbital inclinations (Zhang & Nimmo 2012). This likely eccentricity excitation was also suggested to be related to geological features on Tethys that indicate past tidal heating (Chen & Nimmo 2008). Unlike eccentricities, which damp over time by tides, the moons' orbital inclinations cannot be easily erased once they are excited (Chen et al. 2014). Therefore, it is essential to probe the effects of possible past resonances on the moons' inclinations, using methods that account for the full dynamics of the system.

In this work we use two different numerical integrators, simpl and complex. Simpl is relatively efficient but cannot model close approaches between massive bodies, and it is applicable only to the situations where the moons' orbits are well separated. Complex is capable of integrating close encounters between moons but has a downside of being much slower than simpl.

With a goal of studying the long-term tidal evolution of Saturn's moons, we wrote a new symplectic integrator, simpl (Symplectic Integrator for Moons and PLanets). This is a combination of two Levison–Duncan-type integrators (Levison & Duncan 1994), one for the planets and one for the moons of one planet. This type of integrator is a "democratic heliocentric" (i.e., uses heliocentric coordinates and does not require strict ordering of planets) version of the mixed-variable symplectic integrator of Wisdom & Holman (1991), which assumes that bodies follow Keplerian orbits and implements perturbations as discrete "kicks" executed in Cartesian coordinates. We adapted the specific algorithm from Chambers et al. (2002), who designed it to model the dynamics of planets around one component of a binary star. simpl directly includes mutual moon–moon interactions, as well as solar and planetary perturbations (direct planetary perturbations can be switched off to speed up computations; back-reaction of satellites on planets is ignored). We also include the planet's oblateness, tides raised on the planet by the satellite (as a simple tangential acceleration on the moon in the plane of its orbit, not changing the planet's rotation) and satellite tides (implemented as radial-only kicks), and artificial planetary or satellite migration (to approximate interaction with protoplanetary or protosatellite disks). When calculating tidal accelerations, we assumed that the tidal quality factor Q of Saturn was independent of the satellite's orbital frequency (i.e., all the moons experienced the same Q for Saturn). We extensively validated simpl against published results and analytical estimates and found that all precessional motions and secular tidal effects are reproduced correctly.

To integrate the orbits of moons like Tethys and Dione with orbital periods of 2–3 days, we used a 0.1-day (2.5 × 10⁻⁴ yr) integration time step. This is the time step we used in most simulations shown in this paper, unless otherwise noted. We always included Titan in our simulations, and we usually included only Jupiter and Saturn in the heliocentric part of the integration.

complex (Crossing Orbit Moon and PLanet EXtrapolator) is a more basic type of symplectic integrator, in which the Hamiltonian is not separated into Keplerian and perturbation parts (as in simpl) but into kinetic and potential energy. The exact implementation has been taken from Forest & Ruth (1990) and uses a fourth-order algorithm. complex uses many of the same subroutines as simpl and the same input and output files, except that simpl's Keplerian step has been replaced by a sequence of symplectic mappings in Cartesian coordinates. The subroutines in simpl that implement perturbations from other bodies, the planet's figure, and tidal deformations have been incorporated into the symplectic mapping, in the substep where the momenta are being modified as a result of potential energy.

complex has no preference for Keplerian orbits and can therefore resolve close approaches between the moons naturally. The downside of complex is its much slower speed, requiring about 600 time steps for the shortest orbital period. For the sake of simplicity, we decided to use the same routines to integrate the satellites and the planets in complex, although the planetary orbits could in principle be done more efficiently. Overall, complex is about two orders of magnitude slower than simpl, and we used it only when the satellite system went through an instability. Short collision times between moons (centuries to millennia) keep the relative inefficiency of complex from becoming a major issue. The switching between simpl and complex was done manually, meaning that we go back to the state of a simpl simulation just before the instability and continue the simulation using complex.

While the time step of complex was sufficiently small to resolve close approaches, it was still too large to reliably resolve collisions. Therefore, we shrunk the time step by an additional factor of 20 when two moons were between 10 combined radii from each other. Once again, relatively large collisional cross sections of the moons compared to the sizes of their orbits made this relatively inefficient method workable. Note that this approach would not work for simulations of planets interacting with small bodies, where special treatment for the two bodies undergoing the close encounter (Levison & Duncan 1994) makes more sense than slowing down the whole integration.

The period ratio of Dione and Tethys is currently about 1.45, so, as Tethys evolves outward much faster than Dione, these two moons should have crossed their mutual 3:2 resonance at some point in the past. In this section we use the symplectic numerical integrator simpl to explore the Tethys–Dione 3:2 resonance crossing. The benefit of the numerical approach is that all subresonances, including the inclination-related ones, are fully accounted for. Figure 1 shows a typical evolution of eccentricities and inclinations of Tethys and Dione during the resonance crossing, assuming that both moons' orbits were close to circular and equatorial before the resonant encounter. The moons are originally captured in a subresonance that leads to a monotonic growth of the inclination of Tethys. Once Tethys's inclination grows to about 4°, the resonance becomes chaotic and eccentricities and inclinations of both moons vary rapidly. Eventually, the resonant lock breaks, leaving both orbits significantly eccentric and inclined. While the eccentricities can damp owing to eccentricity tides, the inclinations of Saturn's midsized moons are not expected to evolve appreciably over the age of the solar system (Chen et al. 2014). The final inclination of Dione is always substantial, typically comparable to a degree, while the observed inclination of Dione is only 0028 (Table 1).

Zoom In Zoom Out Reset image size

We used Q/k₂ = 4000 for Saturn in almost all the simulations is this paper, which is equivalent to Q = 1360, as Saturn's Love number is estimated to be k₂ = 0.341 (Gavrilov & Zharkov 1977).³ For the satellites, we used Q = 100 and k₂ = 0.05 in the simulation shown in Figure 1, which implies at least some melting, as expected when eccentricities reach a few percent (Chen & Nimmo 2008). A more detailed model would have tidal parameters varying with tidal heating, but our integrator did not include that capability. Instead, we explored a range of constant tidal parameters for the satellites. Given that Enceladus seems to currently have very low Q/k₂ (if Enceladus's orbit is in equilibrium as proposed by Lainey et al. 2012), we also integrated the resonance crossing using Q = 100, k₂ = 0.5 (Figure 2). Despite strong tidal dissipation, the inclination of Dione is still excited during the crossing. One can in principle imagine an even stronger tidal response from Tethys or Dione (with Q/k₂ ≃ 10), but these would cause Tethys to move inward as a result of satellite tides being faster than the outward migration due to Saturn tides (for eccentricities of the order e = 0.01). Such a reversal of migration would effectively prohibit Tethys from crossing the resonance with Dione and reaching the present configuration (see Equation (1)).

Zoom In Zoom Out Reset image size

We found it impractical to directly model resonance crossing in the "classical" case when Q = 18,000 for Saturn, as the integrations would need to be extended to many hundreds of Myr. However, chaotic excitation of both the eccentricity and inclination is expected to be even stronger when the resonance is crossed more slowly. If we assume a larger tidal Q for Saturn, a wider range of tidal responses from Tethys would overwhelm Saturn's tides and lead to divergent migration, making the resonance an impassable barrier for Tethys. The ratio between the migration rate due to tidal dissipation within the satellite and the migration rate due to dissipation within the planet is

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{a}}_{\mathrm{sat}}}{{\dot{a}}_{\mathrm{pl}}}=-7{e}^{2}\displaystyle \frac{{Q}_{\mathrm{pl}}^{*}}{{Q}_{\mathrm{sat}}^{*}}{\left(\displaystyle \frac{{M}_{\mathrm{pl}}}{{M}_{\mathrm{sat}}}\right)}^{2}{\left(\displaystyle \frac{{R}_{\mathrm{sat}}}{{R}_{\mathrm{pl}}}\right)}^{2},\end{eqnarray} \tag{ 1 }$$

where subscripts "sat" and "pl" refer to the satellite and the planet, M, R, and e are mass, radius, and orbital eccentricity, respectively, while Q* = Q/k₂ is the measure of tidal dissipation. For e = 0.01 (typical of eccentricities excited by the Tethys–Dione resonance) and the size and mass of Tethys, we get ${Q}_{\mathrm{pl}}^{*}/{Q}_{\mathrm{sat}}^{*}$ = 0.03 for outward migration to stop. If ${Q}_{\mathrm{pl}}^{*}\;=\;4000$, then ${Q}_{\mathrm{sat}}^{*}\gt 120$ is necessary for Tethys to continue outward migration. This limit is comparable to the current tidal response of Enceladus inferred by Lainey et al. (2012) by assuming an equilibrium. However, if ${Q}_{\mathrm{pl}}^{*}\;$ = 50,000, then ${Q}_{\mathrm{sat}}^{*}\gt 1500$ for the migration to stay convergent, which likely excludes any major melting within Tethys, despite significant eccentricity (if Q = 100, then k₂ would have to be less than 0.07, indicating a relatively rigid body). Without going into modeling the actual physical response of the satellite, we can just use the above considerations to exclude the parameter space with ${Q}_{\mathrm{pl}}^{*}/{Q}_{\mathrm{sat}}^{*}\gt 0.03$ from our consideration, as those parameters would forever trap Tethys on the inner side of the 3:2 resonance with Dione. The only exception is ${Q}_{\mathrm{pl}}^{*}\gt 2\times {10}^{5}$, in which case the resonance was never crossed, regardless of satellite tides. However, this extremely slow evolution rate is about 40 times too slow to account for the heating of Enceladus, and we prefer the interpretation that the resonance was never crossed because the moons are young.

It is possible that the inclination of Tethys was already excited before it encountered the resonance with Dione. Analytical considerations open a possibility that a pre-inclined Tethys could have conceivably avoided capture into the inclination-type resonance evident in Figure 1 and the associated inclination increase. We tested this possibility using initial inclinations of 1° and 03 for Tethys in multiple simulations. Figures 3 and 4 each plot the results of one simulation with initial i = 1° and i = 03 for Tethys, respectively. For initial i = 03 resonance capture still happens, although it usually lasts a much shorter time than in Figure 1, while for initial i = 1° the moons enter the chaotic behavior at the very beginning of the resonance crossing. In all cases the resulting final inclination of Dione is much higher than the observed one (Figure 5). We also conducted simulations that start with one moon's eccentricity being 0.01. While simulations with initially eccentric Dione did lead to the onset of chaos earlier in the simulation, those with initially eccentric Tethys show no difference from the circular case, chiefly because most of this initial eccentricity damps by the time chaos is triggered (despite Q/k₂ = 10⁴ being used in this run). In all cases the inclinations of both Tethys and Dione are excited well in excess of the observed values (Figure 5).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The numerical experiments shown in this section suggest that the inclinations of both Tethys and Dione are excessively excited by the crossing of their 3:2 resonance, regardless of the assumed initial conditions and tidal parameters. In the next section we will discuss possible resonances that could have affected Dione and Tethys later on, but they cannot make a past 3:2 resonance crossing compatible with the present orbits. The solid line in Figure 5 plots the combinations of inclinations that have the observed out-of-plane angular momentum deficit⁴ (AMD), and they are all well short of the inclination AMD seen among the simulation outcomes. Capture into additional MMRs could only increase the AMD of Tethys and Dione, while the secular resonance involving these two moons (see next section) would conserve AMD associated with inclinations. Likewise, chaotic interaction within an MMR is extremely unlikely to decrease the combined inclinations but tends to redistribute inclinations on short timescales and increase AMD over time. This applies to subsequent resonances with Enceladus and Mimas; Mimas does have a large inclination but small AMD compared with Tethys owing to its small mass, and it could not have absorbed out-of-plane AMD from Tethys and Dione. The fact that both Tethys and Dione are always excessively inclined after their 3:2 MMR crossing is a strong constraint on their relative locations at formation, and we will discuss the implications in the following sections.

In the relative chronology of resonance crossings in the Saturnian system, the Dione–Rhea 5:3 MMR would happen shortly after the 3:2 Tethys–Dione resonance. We have concluded in the previous section that the passage through 3:2 Tethys–Dione resonance did not happen, and here we explore the possibility that Tethys formed outside the 3:2 resonance with Dione and Dione formed interior to the 5:3 resonance with Rhea, and that Dione and Rhea have since crossed their 5:3 resonance.

Figure 6 shows one simulation of the 5:3 Dione–Rhea MMR crossing. In the beginning Dione is captured in the inclination-type resonance, which eventually breaks and a chaotic phase follows, during which the eccentricities and inclinations of both bodies are excited. This part of the evolution is expected, and while it could explain the excitation of Rhea's inclination (i = 035), it is seemingly at odds with the very low inclination of Dione (i = 003). While the eccentricities are also excited, they can be expected to damp later through satellite tides. However, at 24 Myr we see a striking exchange of eccentricities and inclinations between Tethys and Dione. Both the eccentricity and inclination of Tethys get excited this way, while those of Dione are reduced. We find that this unexpected dynamical event is caused by a previously unidentified Tethys–Dione secular resonance.

Zoom In Zoom Out Reset image size

The secular resonance seen in these simulations has the resonant argument Ψ = ϖ₂ − ϖ₁ + Ω₂ − Ω₁, where ϖ and Ω are the longitudes of the pericenter and the node, and the subscripts 1 and 2 refer to Tethys and Dione, respectively. The pattern of signs in the resonant argument indicates that in this resonance the eccentricity and inclination of each moon should be correlated, and anticorrelated with those of the other moon, consistent with numerical integrations showing Tethys getting more inclined and eccentric, with Dione's orbit becoming more circular and planar. The resonance appears to break when either the eccentricity or inclination of Dione reaches zero, with the latter case being able to explaining the very low present inclination of Dione. Note that in Figure 6 it is eccentricity, rather than inclination, that reaches zero first and breaks the secular resonance.

The type of secular resonance seen here affecting Tethys and Dione has not, to our knowledge, been described in the literature. Therefore, here we briefly discuss why it arises soon after the Dione–Rhea 5:3 MMR. The resonant argument Ψ is a combination of the arguments ϖ + Ω for each moon, which are typically very slowly varying for a regular satellite, as $\dot{\varpi }\simeq -\dot{{\rm{\Omega }}}$ for low-e, low-i orbits chiefly perturbed by planetary oblateness. Usually, the leading term in $\dot{\varpi }+\dot{{\rm{\Omega }}}$ is $9/4{J}_{2}^{2}{(a/R)}^{4}n$ (Murray & Dermott 1999), where J₂ is the standard oblateness moment and R the radius of Saturn, and a and n are the moon's semimajor axis and mean motion. This result does not depend on the definitions of orbital elements used (geometric versus osculating), which affect individual expressions for $\dot{{\rm{\Omega }}}$ and $\dot{\varpi }$ but cancel out in their sum (Greenberg 1981). This term currently leads to periods of about 5000 and 20,000 yr for ϖ + Ω of Tethys and Dione, respectively, making $\dot{{\rm{\Psi }}}\lt 0$. Clearly, another source of precession is necessary to bring these two moons into a secular resonance.

Precession due to oblateness is also affected by eccentricity and inclination. From Danby (1992) we have

$$\begin{eqnarray}\dot{\varpi }+\dot{{\rm{\Omega }}}=\dot{\omega }+2\dot{{\rm{\Omega }}} & = & \displaystyle \frac{3{{nJ}}_{2}}{{a}^{2}{(1-{e}^{2})}^{2}}\left(1-\displaystyle \frac{5}{4}{\mathrm{sin}}^{2}i-\mathrm{cos}i\right)\\ & \simeq & \displaystyle \frac{3{{nJ}}_{2}}{{a}^{2}}\left(-\displaystyle \frac{3}{4}{\mathrm{sin}}^{2}i+O({i}^{4})\right).\end{eqnarray} \tag{ 2 }$$

Here ω = ϖ − Ω is the argument of pericenter, and e and i are the moon's orbital eccentricity and inclination (to the planet's equator), respectively. For i = 1°, we find that the quantity in parentheses dependent on inclination is only −2 × 10⁻⁴, which for Tethys would produce retrograde precession of $\dot{{\rm{\Omega }}}+\dot{\varpi }$ with a period of 11,000 yr. This is a potentially significant term, but it would add a negative term to $\dot{{\rm{\Psi }}}$ before the resonance is encountered, as it is Dione, rather than Tethys, that must be inclined before the resonance capture occurs. Therefore, this term cannot be the cause of the resonance being encountered in the first place. However, inclination-dependent terms may be important for maintaining the resonance, as they will slow down the excess apsidal precession of Tethys as it acquires inclination.

Since the secular resonance is encountered near Dione's 5:3 resonance with Rhea, it appears that near-resonant terms are important for establishing the secular commensurability. With this being the second-order resonance, six different resonant terms can affect the precession of Dione. After some numerical testing, we have established that this additional precession is sensitively dependent on Dione's distance from the outermost subresonance (associated with the ${e}_{2}^{2}$ term) and is also strongly dependent on Dione's eccentricity. These terms appear to accelerate the precession of Dione's longitude of pericenter but do not affect the node as much, producing a significant change in the secular resonant argument Ψ.

Figure 7 shows a simulation where Dione and Rhea are placed just outside their 5:3 resonance (i.e., after an assumed 5:3 resonance passage), with Tethys having a near-circular and planar orbit and Dione having e = 0.02 and i = 08 (Rhea, not plotted, had initial e = 0.006 and i = 03). About 1.4 Myr into the simulation the secular resonance is established, and it lasts until about 8.5 Myr, when the average inclination of Dione drops below 01. While this is a somewhat artificial setup (we assume the outcome of the 5:3 Dione–Rhea resonance, without simulating it), it clearly shows the potential of the secular resonance to produce the current orbits of Tethys and Dione, with the right initial conditions.⁵ We find that resonance capture is almost certain (assuming Q/k₂ = 4000) if the eccentricity of Dione is larger than 0.01. The current low inclination of Dione indicates that, if it was reduced by this secular resonance, the resonance had to break when the inclination of Dione reached zero (as opposed to the situation shown in Figure 6, where eccentricity reaches zero first). As ${\dot{e}}_{2}\;=\;{\dot{i}}_{2}$ in this secular resonance, and the angular momentum is conserved, this means that the initial ${e}_{2}\gt {({a}_{1}/{a}_{2})}^{1/4}{({m}_{1}/{m}_{2})}^{1/2}\mathrm{sin}{i}_{1}\;=\;0.012$. This gives us a requirement for post-5:3 resonance orbits of Dione and Rhea: e₂ > 0.012, e₃ unconstrained, i₂ = 07, and i₃ = 033 (current value). We assume low eccentricity and inclination for Tethys at this time, as higher eccentricities of Tethys often lead to greater chaos and sometimes overlap between the Tethys–Dione 3:2 and Dione–Rhea 5:3 resonances, which is long-lasting, is hard to break, and excites inclinations beyond the observed values.

Zoom In Zoom Out Reset image size

While the required initial conditions for the secular resonance fall within the range of stochastic outcomes of the 5:3 Dione–Rhea resonance, we could not find a combination of initial conditions and tidal parameters that lead to such a state with high probability. We find that the initial capture into the ${i}_{2}^{2}$ subresonance of the 5:3 MMR (involving only the inclination of Dione) is certain for very low initial inclinations (i₂ ≃ 002). This subresonance, however, tends to break only at i₂ = 14, which is well above the value required for the secular resonance. An initial inclination i₂ = 01 leads to a jump through the ${i}_{2}^{2}$ subresonance and capture into the i₂i₃ subresonance, in which the inclinations of both Dione and Rhea grow. This resonance usually breaks in the i₂ = 02–07 range, approaching the target inclination needed for the secular resonance, while the inclination of Rhea reaches i₃ = 01–04, approaching the current values. Figure 8 shows one simulation using initial i₂ = 01 that results in a very low inclination of Dione and a substantial inclination of Tethys after the secular resonance.

Zoom In Zoom Out Reset image size

Figure 9 shows the outcome of several sets of simulations starting with i₂ = 01, with all other eccentricities and inclinations (except those of Titan) being close to zero. It is clear from the figure that the amount of excitation needed to explain the present system roughly matches the outcomes of the Dione–Rhea MMR. While no single simulation exactly reproduces the current orbits, very low inclinations of Dione and significant inclinations of Tethys are often produced through the secular resonance.

Zoom In Zoom Out Reset image size

While the mechanism involving the 5:3 Dione–Rhea resonance followed by the Tethys–Dione secular resonance is a very promising avenue for producing the observed inclinations, we were not able to achieve this outcome with any reasonable probability using simpl. One possible explanation is that the process is so stochastic that we cannot expect to exactly reproduce the real history of the system. We think it is still too early to reach such a conclusion, and we hope to explore this question in future work. The major limitations of our present approach include the following: (1) we assume that the tidal properties of the satellites are constant over time, despite varying eccentricities and tidal heating; (2) we assume that Saturn's tidal Q is frequency independent, which has implications for the exact timing of recent resonance crossings; and (3) while our integrator allows it, we did not explore different initial tidal parameters for the three moons (i.e., we assumed that they were all equally rigid or dissipative, which is likely to be wrong). A clear step forward would be incorporation of a geophysical model into the integrator that would adjust tidal parameters depending on the tidal heating and other processes the moon is subject to.

The chief conclusion of this section is that the past 5:3 MMR between Dione and Rhea likely did happen, unlike the 3:2 Tethys–Dione resonance discussed in the previous section. This gives us a very well defined relative age of Tethys, Dione, and Rhea, which is related to their absolute age through the controversial tidal Q of Saturn. As for the absolute ages, the classical Q = 18,000 (Meyer & Wisdom 2007) and the "fast" Q = 1700 (Lainey et al. 2012) would place the formation of Tethys, Dione, and Rhea at about 1.1 Gyr and 100 Myr ago, respectively, with both dates being much younger than the age of the solar system. In the following sections, we will briefly explore the possible dynamical mechanisms that could have led to the late formation of the present set of midsized icy moons of Saturn.

As we are proposing that the moons of Saturn (out to Rhea) formed much later than the planet, we need to address the plausibility and basic physics of such an event. First of all, the flatness of Saturn's regular satellite system indicates formation from a thin disk. Since we think that Tethys, Dione, and Rhea all formed close to their present orbits, this disk would need to extend out to about 10 Saturn radii and compose the combined mass of the midsized moons. Since the Titan–Hyperion resonance appears very old (Greenberg 1973a; Ćuk et al. 2013), we assume that Titan (and other more distant moons) predated this protosatelite disk from which Tethys, Dione, and Rhea formed. This does not necessarily mean that Titan must be primordial, as it may have originated during another, more ancient satellite instability event (Hamilton 2013), which is outside the focus of this paper. As Saturn has youthful and dynamic rings surrounded by numerous small moons (Cuzzi et al. 2010), we propose that the rings also formed in the same event as Tethys, Dione, and Rhea. This means that the protosatellite disk discussed here extended all the way inward to the Roche limit (possibly owing to viscous spreading). Unlike a primordial circumplanetary disk, this disk would be composed exclusively of solids, save for short-lived vapor created in collisions. In this section we discuss the likely lifetime and fate of this disk, while in Sections 7 and 8 we will address the origin of the disk.

Burns & Gladman (1998) used numerical simulations to study the fate of ejecta in the Saturnian system in order to constrain possible hazards for the Cassini spacecraft. They found that the moons accrete any material initially found on crossing orbits in mere centuries (also see Alvarellos et al. 2005), and theoretical predictions of their growth rate from debris are also on the order of millennia (see Goldreich et al. 2004; Chiang & Laughlin 2013). We can use the expression in Chiang & Laughlin (2013) (also see Goldreich et al. 2004) to calculate the coagulation time:

$$\begin{eqnarray}&&t\simeq \displaystyle \frac{1}{{F}_{\mathrm{grav}}}\displaystyle \frac{{\rho }_{\mathrm{core}}{R}_{\mathrm{core}}}{{\sigma }_{\mathrm{solid}}{\rm{\Omega }}}\end{eqnarray} \tag{ 3 }$$

where ρ and σ are the volume and surface densities of a solid core and the disk, respectively, and F_grav is the gravitational focusing factor (we will assume F_grav = 1). If we estimate σ_solid = 4 × 10²¹ kg π⁻¹ (6 × 10⁸ m)⁻² ≃ 4000 kg m⁻², ρ_core = 1200 kg m⁻³, and Ω = 2 × 10⁻⁵ s⁻¹ (corresponding to a 4-day orbital period), we get

$$\begin{eqnarray}&&t\simeq 0.5\;\mathrm{yr}\;{\mathrm{km}}^{-1}{R}_{\mathrm{core}}.\end{eqnarray} \tag{ 4 }$$

Therefore, it would take less than 1000 yr to grow 1000 km satellites, in agreement with the numerical integrations of debris removal by Burns & Gladman (1998). However, there are serious problems with this approach. The first is that the present satellites have well-separated Hill spheres and would therefore not sweep up all of the material between them. Assuming a uniform disk, we can estimate the number (n) and size (r) of the satellites assuming that they are separated by 5 Hill radii:

$$\begin{eqnarray}&&r=5{{nR}}_{\mathrm{hill}}=5{nr}{\left(\displaystyle \frac{M}{3{M}_{\mathrm{pl}}}\right)}^{1/3}=5{nr}{\left(\displaystyle \frac{{M}_{\mathrm{disk}}}{3{{nM}}_{\mathrm{pl}}}\right)}^{1/3}\end{eqnarray} \tag{ 5 }$$

so therefore

$$\begin{eqnarray}&&n={5}^{-3/2}{\left(\displaystyle \frac{3{M}_{\mathrm{pl}}}{{M}_{\mathrm{disk}}}\right)}^{1/2}\simeq 60,\end{eqnarray} \tag{ 6 }$$

where we have assumed M_disk ≃ M_Rhea + M_Dione + M_Tethys. This makes the mass of a typical protosatellite M = M_disk/n = 0.7 × 10²⁰ kg, comparable to that of Enceladus. These protosatellites would take significantly longer to merge into the larger moons we observe today. The mechanism for their destabilization and mergings would likely be chaotic interaction between the moons, many of which would likely be in overlapping MMRs with several others. Additionally, there would be some fragmentation during their accretion and mergings, and rings of debris would likely occupy quasi-stable areas between neighboring moons' cleared zones. Interaction with this debris would also affect the protosatellites' orbits, leading to orbit crossings and mergers. Direct numerical modeling is clearly needed to explore the evolution of such a system.

Another caveat with accretion is outside perturbations, mainly from Titan. As Burns & Gladman (1998) show, the orbits close to Titan's 3:1 and 2:1 resonances are significantly more perturbed than those in between midsized moons. This excitation would inhibit accretion, at least in some locations, potentially leading to a relatively long-lived debris ring outside the orbit of Rhea. Eventually, this disk would spread and disperse owing to interparticle scatterings and collisions. Crida & Charnoz (2012) find that the timescale for the lifetime of planetary rings, based on numerical experiments, can be expressed as

$$\begin{eqnarray}&&{T}_{\mathrm{disk}}=\left(\displaystyle \frac{\pi }{52{({R}_{{\rm{H}}}/{R}_{\mathrm{part}})}^{5}}\right)\displaystyle \frac{{M}_{\mathrm{pl}}^{2}}{{\sigma }_{\mathrm{disk}}^{2}{r}_{\mathrm{disk}}^{4}}{T}_{\mathrm{orbit}}\end{eqnarray} \tag{ 7 }$$

where R_part and R_H are a particle's physical and Hill radius, respectively. Assuming the disk radius of r = 6 × 10⁸ m (10 Saturn radii) and a density of 500 kg m⁻³ for disk particles, (R_H/R_part) ≃ 6, and T_orbit = 5 days,

$$\begin{eqnarray}&&{T}_{\mathrm{disk}}=3\times {10}^{11}\;{\rm{yr}}\displaystyle \frac{1}{{\sigma }_{\mathrm{disk}}^{2}{[\mathrm{kg}{{\rm{m}}}^{-2}]}^{2}}.\end{eqnarray} \tag{ 8 }$$

If we assume a surface density σ_disk = 4000 kg m⁻², the spreading time according to the above formula is about 2000 yr. After the surface density decreases by a factor of a few, the disk can last tens of kiloyears. At that point it would still contain a mass on the order of M = 10²¹ kg, making it possible to force protosatellites (or satellite fragments) of comparable mass to migrate as the disk is spreading. While the material on the inner edge would likely be accreted by Rhea, the large gap between the inner system and Titan would leave space for any protosatellites (or fragments of disrupted satellites) to migrate and potentially interact with Titan through MMRs.

The outer part of the disk, outside the 3:1 resonance with Titan (and particularly around the 2:1 resonance at 12.7 Saturn radii), would experience strong forced eccentricity due to Titan and may endure as a debris disk for a much longer time (Burns & Gladman 1998). A satellitesimal (or a fragment) that is ejected outward may be able to migrate by interacting with the spreading disk and become trapped in an MMR with Titan. The migration of such an object due to disk spreading could be much faster than that due to tidal migration and could excite Titan's eccentricity to the observed value in a relatively short amount of time. In the simulation shown in Figure 10, we introduced an artificial fast migration of a 4 × 10²⁰ kg moon (0.3% of Titan's mass) starting around 13 Saturn radii, just outside Titan's 2:1 resonance. Our hypothesis is that there would be a spreading ring of excited, nonaccreting debris around the location of the resonance. The initial eccentricity of the moon was set to e = 0.1, as we would expect this moon's (or fragment's) orbit to be excited. The moon in question migrates out until it is caught in the 5:3 resonance with Titan. In the first 1000 yr after resonance capture, it is the eccentricity of the inner moon that grows, indicating a Lindblad resonance (involving the smaller body's eccentricity). However, very soon the system settles into a "co-rotation"-type resonance (involving the larger body's eccentricity) and stays in it for the rest of the simulation. This behavior is unexpected and is likely related to the relatively rapid rate of evolution, which stabilizes formally overlapping subresonances.

Zoom In Zoom Out Reset image size

In 50,000 yr, the eccentricity of Titan grows to about 0.05, well above its observed value. Eventually, the resonance would break as the eccentricities of both bodies become high (the inner moon is indirectly affected by Titan's eccentricity), or the migration slows down and long-term chaos takes over. Figure 11 shows a similar simulation involving a twice larger inner moon (M = 8 × 10²⁰ kg) and a similar capture into the 8:5 resonance with Titan. After the resonance breaks, the inner moon collides with either Titan or the disk.

Zoom In Zoom Out Reset image size

Reabsorption of the protosatellite into the disk would have no consequences for Titan or Hyperion. However, close encounters followed by a collision with Titan would perturb the orbits of both Titan and Hyperion. If the Titan–Hyperion resonance is several gigayears old (Greenberg 1973a; Ćuk et al. 2013), its continued existence might be in conflict with any destructive events in recent history. We therefore used complex to run a set of 100 simulations of scatterings between Titan, Hyperion, and a 9 × 10²⁰ kg satellite. The "rogue" satellite was initially placed on an eccentric orbit with a = 13.2 R_S and e = 0.53, with each run having a different longitude of pericenter or a different initial longitude for the inner satellite. Titan and Hyperion were initially placed deep in the resonance, with Titan on a circular orbit and Hyperion having e = 0.1. In 85 (out of 100) 1000 yr simulations using complex, Titan and the inner moon collided, and in 76 of these Titan and Hyperion maintained a stable 4:3 resonance. The stability of the resonance was determined by integrating post-merger Titan and Hyperion for 10,000 yr using simpl.

The average orbits of Titan and Hyperion in the post-merger simpl simulations are shown in Figure 12. The eccentricity of Titan is excited to about 0.01–0.02 on average (about half of its present value), requiring a different excitation mechanism (such as the MMR described in this section) to explain its present orbit. The main (resonant) component of Hyperion's eccentricity stays close to 0.1, with the other ("slow") component being proportional to Titan's eccentricity. The resonant argument is excited, and the present libration semiamplitude of 59° is firmly within the ensemble of outcomes.

Zoom In Zoom Out Reset image size

We conclude that a resonant interaction with a moon as small as 4 × 10²⁰ kg, migrating owing to disk spreading, can excite Titan's eccentricity to the current value. Even if this moon subsequently impacts Titan, the Titan–Hyperion MMR is likely to survive this event. We numerically confirm resonance survival for a "rogue moon" mass of 9 × 10²⁰ kg, and it is certain that collisions with smaller moons would have an even less disruptive effect on the Titan–Hyperion resonance. An eccentricity excitation only about 100 Myr ago would explain the surprisingly large free eccentricity of Titan despite its likely dissipative surface and interior (Sagan & Dermott 1982; Sohl et al. 1995).

In the previous sections we find that Saturn's present-day moons interior to Titan's orbit likely accreted much after the planet's formation, presumably from a disk, which covered the locations of the moons from Mimas to Rhea. As we discussed in Section 5, this disk would last only millennia, possibly somewhat longer in locations where Titan's perturbations inhibit accretion. Therefore, we think that the disk was a temporary feature, and we propose that it originated in the destruction of a previous satellite system. While we argue that Saturn's present midsized satellite system is relatively young, its expected lifetime is measured in Gyr,⁶ even in the case of "fast tides" proposed by Lainey et al. (2012). Therefore, we argue that the most long-term arrangement of the mass currently in the midsized moons of Saturn is a system of satellites. As this previous system of satellites would have had the extent and total mass of the present midsized icy moon system, we will assume that it also had most of its mass in the outermost two satellites, like Dione and Rhea today (this architecture is also most stable against tidal evolution).

Is there a way of destabilizing a midsized satellite system? We expect that any perturbations that can lead to large orbital eccentricities and moon–moon collisions would have to be resonant in nature. In this and the following sections, we will discuss several candidate resonances for destabilization of the putative previous-generation satellite system of Saturn. Assuming Q ≃ 1700 for Saturn, a Rhea-sized moon would evolve to about the present location of Rhea over 4.5 Gyr, while a somewhat smaller moon (like Tethys or Dione) would converge on positions very close to, but interior to, Rhea's present orbit. This arrangement, with a Rhea-sized moon close to Rhea's present orbit and a "proto-Dione" on a converging interior orbit, is the starting point for our studies of past dynamical instability.

The simplest case would be an MMR like that of Tethys and Dione described in Section 3. However, the relatively small masses of the inner moons make it impossible to excite their eccentricities to the high values required for orbit crossing (at least e = 0.1–0.2). Tidal evolution while the moons are captured in an MMR can in principle lead to monotonic eccentricity growth, but tidal dissipation within the moons can easily limit eccentricity growth for realistic physical parameters. A chaotic resonance between a large and a much smaller moon could excite the smaller moon's eccentricity enough to cause a collision, but such a mismatched collision is unlikely to lead to disruption of a Dione- or Rhea-sized moon (Leinhardt & Stewart 2012).

Resonances between a midsized satellite and Titan have the potential of being more powerful, but they are also unlikely to have been the cause of a past instability. We find that the 4:1 resonance with (noneccentric) Titan at about 8 Saturn radii leads to temporary capture, after which the inner moon acquires an eccentricity of e ≤ 0.01 and i = 03. Another major resonance close to the midsized moons is the 3:1 with Titan, at 10 Saturn radii. This resonance is strongly chaotic, making the inner moon significantly eccentric and plausibly leading to instability. However, this resonance overlaps with the 4:1 resonance with Hyperion (which is in the 3:4 resonance with Titan). As Hyperion is small and does not damp eccentricity, a Rhea-sized moon in the 3:1 resonance with Titan quickly destabilizes Hyperion. Owing to the survival of Hyperion, we can exclude the 3:1 resonance with Titan as the source of a past instability.

Another set of resonances in the region inhabited by the midsized satellites is the so-called evection-type resonances. These are commensurabilities between a satellite's apsidal precession and the Sun's (or another distant perturber's) apparent planetocentric orbital period. Given Saturn's oblateness and heliocentric orbit, the strongest of these resonances happen for moons at about 8.3 Saturn radii, which have an apsidal precession period equal to Saturn's 29.5 yr orbital period. Additional harmonics with evection happen at 6.9 Saturn radii (precession is 1/2 of Saturn's year) and 7.3 Saturn radii (precession equals 2/3 of Saturn's year).

Figure 13 shows a moon encountering the main evection resonance assuming that its tidal response stays limited (M = 16 × 10²⁰ kg, Q = 100, k₂ = 0.05 for the satellite). The eccentricity of the moon is strongly excited, and so is the inclination. However, since this resonance affects only the eccentricity but not the mean motion of the satellite, and the satellite tides are too weak to reverse the evolution (as we assumed k₂ = 0.05), the satellite eventually passes the resonance without being destabilized (as long as it has no very close neighbors). After crossing the resonance, both eccentricity and inclination are high. While eccentricity can damp, inclination will not, implying that any moon that passed through the evection resonance should have an inclination of several degrees. This is well in excess of the current inclination of Rhea (the only moon beyond the evection resonance), implying that an intact Rhea with these tidal properties could not have formed close to the rings and migrated out to its present location, as suggested by Charnoz et al. (2011).

Zoom In Zoom Out Reset image size

If we assume large-scale melting within the moon whose orbit is excited by the evection resonance, its orbital behavior changes. As shown in Figure 14 (where we assumed a much larger k₂ = 0.5 for the moon), the evection resonance excites a large eccentricity, which then damps and in the process makes the moon migrate inward. After the eccentricity is damped, the outward migration continues and the cycle repeats. Each next "bounce" from the resonance will not be identical, as the inclination also evolves. We see no reason why this process cannot continue for 1 Gyr or more, trapping the satellite close to the original location of the evection resonance.

Zoom In Zoom Out Reset image size

The difference between the scenarios shown in Figures 13 and 14 depends only on the internal properties of the satellite. Regardless of the outcome (resonance crossing or long-term stalling of outward tidal evolution), the evection resonance has a potential of exciting satellite orbits to eccentricities of e = 0.1 or more in very short amounts of time. Before considering the effect of evection on two interacting satellites, we will briefly address some higher-order harmonics of evection.

The evection resonance in the Saturnian system is notably chaotic and tends to affect both eccentricity and inclination. This is in contrast to the evection resonance in the early Earth–moon system, where it affects only eccentricity, and can lead to stable capture (Touma & Wisdom 1998; Ćuk & Stewart 2012) or a stable near-resonant equilibrium (Wisdom & Tian 2015). We find that this complex behavior is caused by the evection resonance's overlap with its nearby harmonics, which affect inclination. The principal inclination harmonic has the resonant argument λ_Sun − ϖ_Sun + Ω₂ − Ω_Sun, where Ω₂ is the moon's longitude of ascending node, while Ω_Sun is the node of the Sun's apparent orbit measured relative to the moons' Laplace plane (in practice, this angle is equivalent to the equinox of Saturn). This harmonic can be seen as a combination of the annual equation and the principal secular perturbation in inclination (Brouwer & Clemence 1961). This argument's strength varies with Saturn's eccentricity, and therefore it cannot lead to stable long-term capture, but it can induce significant variation of inclination, as seen in Figures 13 and 14. In Figure 14, the inclination is stable after about 100 Myr, indicating that the moon may be in a state where it is affected only by the principal evection resonance.

The strongest among the 2:1 evection harmonics at 6.9 Saturn radii involves not only the perturbed satellite and the Sun but also Titan. This is due to its resonant argument 2λ_Sun − ϖ − ϖ_T, where ϖ designates longitude of pericenter and λ the mean longitude. This argument combines the basic evection term 2λ_Sun − 2ϖ with the moon–Titan secular term ϖ − ϖ_T (for a summary of different perturbation terms see Brouwer & Clemence 1961). We find that midsized icy moons can get captured in this resonance, and that their further evolution through the resonance increases the eccentricity of Titan. This resonance may seem like an attractive way of generating Titan's eccentricity (Ćuk et al. 2013), but the resonance always breaks as a result of satellite tides when the smaller moon's eccentricity reaches a few percent. The increase in Titan's eccentricity is always a factor of a few too small to explain the observed value, and this resonance never seems to lead to instability.

The strongest 3:2 evection harmonic at 7.3 Saturn radii has the resonant argument 3λ_Sun − 2ϖ − ϖ_Sun. This argument can be viewed as a combination of the evection and the annual equation λ_Sun − ϖ_Sun (Brouwer & Clemence 1961). As its argument includes ϖ_Sun, the 3:2 evection harmonic term's strength is proportional to Saturn's heliocentric eccentricity, which varies by an order of magnitude (0.01–0.09) during a 50,000 yr secular cycle (Murray & Dermott 1999). This means that the resonance is weak during a low point in Saturn's eccentricity, and our simulations show that a capture into this resonance breaks relatively fast, with the moon moving past the resonance with an eccentricity of a few percent and an inclination of a degree or two. We find that this is an unlikely candidate for instability.

In order to produce a disk, we need a collision between two large moons, probably somewhat similar to Dione and Rhea. As the evection resonance is the most promising source of excitation, we envisage a scenario in which the outer moon is in or near the evection resonance distance, while the inner moon is somewhat closer to Saturn. This configuration can arise either by moons migrating together in an MMR or by the inner moon catching up to the outer one. In general, similar-sized moons would tend to converge on close orbits owing to tides, and the real Dione and Rhea (mass ratio 1:2.4) are evolving toward a configuration more compact than their 4:3 mutual resonance. Another possibility is that the outer moon is caught in a cycle of evection resonance "bounces" (Figure 14), while the inner moon is catching up.

In Figure 15 we show an integration (using simpl) of convergent migration of a "proto-Dione" (with a combined mass of Dione and Tethys) and "proto-Rhea" (with the mass of Rhea) into their mutual 4:3 resonance. Once the outer moon enters the evection resonance, the satellites become dynamically coupled to Saturn's heliocentric orbit and their eccentricities are ultimately excited to large values (e = 0.1–0.2), leading to orbit crossing and a collision.

Zoom In Zoom Out Reset image size

In this simulation, the 4:3 MMR capture happens first, and only then do the moons enter the evection resonance. While the moons are in the 4:3 MMR, they are also in the secular resonance ϖ₂ − ϖ₁ (i.e., apsidal anti-alignment; Figure 16). Owing to the stabilizing effects of the secular resonance, the moons are capable of simultaneously being in stable e₁ and e₂ subresonances of the 4:3 resonance, without the chaos usually accompanying resonance overlap (i.e., the arguments 4λ₂ − 3λ₁ − ϖ₁ and 4λ₂ − 3λ₁ − ϖ₂ are both librating). This configuration is apparently stable until the outer satellite reaches the evection resonance.

Zoom In Zoom Out Reset image size

The next stage of their evolution, starting at about 24 Myr, is chaotic, and during this time the moons are not in either 4:3 or (apsidal) secular resonance. Chaos is driven by both the evection resonance proper and its associated harmonics. While the main evection resonance 2λ_Sun − 2ϖ₂ and its inclination-related harmonic λ_Sun − ϖ_Sun + Ω₂ − Ω_Sun affect the outer satellite, the 2:1 evection harmonic λ_Sun − ϖ₂ − ϖ_T affects the inner satellite (Ω is the longitude of the node, and the subscript T refers to Titan). Figure 17 shows the 10 kyr interval 26 Myr into the simulation, during which both the main evection argument and its 2:1 harmonic reverse the direction of circulation, and this behavior is correlated with the eccentricities of the moons (which continually exchange angular momentum). The presence of the 2:1 evection harmonic at the location of the inner moon is a simple consequence that two moons around an oblate body have an apsidal precession period ratio of about two when they are close to the 4:3 resonance. Since $\dot{\varpi }\propto {a}^{-7/2}$,

$$\begin{eqnarray}&&\displaystyle \frac{\dot{{\varpi }_{1}}}{\dot{{\varpi }_{2}}}={\left(\displaystyle \frac{{a}_{2}}{{a}_{1}}\right)}^{3.5}={\left({\left[\displaystyle \frac{4}{3}\right]}^{2/3}\right)}^{7/2}=1.96;\end{eqnarray} \tag{ 9 }$$

therefore, if the outer moon's pericenter precesses at a rate close to Saturn's mean motion, the inner moon would naturally precess at about twice that rate (barring other complications), exposing itself to perturbations arising from 2:1 near-resonance with the heliocentric motion. The importance of the 2:1 harmonic depends on the eccentricity of Titan, which was 4 × 10⁻³ in this simulation (i.e., one order of magnitude lower than observed). If Titan efficiently damped any eccentricity it might have had before the last rearrangement of the system, this resonance may have been much less prominent. In general, we find that the main evection resonance is the most important one for the eventual fate of the system.

Zoom In Zoom Out Reset image size

The two moons, while having large and rapidly varying eccentricities and inclinations, remain stable for several Myr. Eventually the pair enters a state in which the 2:1 harmonic is less important and the evection resonance dominates. Just before instability, the outer moon gets a major eccentricity boost from the evection resonance (Figure 18). Soon after that, close encounters begin and the two orbits are irreversibly coupled. At about 26.6 Myr, as the orbits become crossing, we pass the simulation from simpl to complex. We find that the collision typically occurs in about 1000 yr. This collision, at about 3 km s⁻¹, is just energetic enough to be erosive and can lead to the formation of a disk (Leinhardt & Stewart 2012).

Zoom In Zoom Out Reset image size

We also achieved instability in a simulation where the moons are initially not in 4:3 resonance, but on orbits converging as a result of tidal evolution (this time the masses are M₁ = 11 × 10²⁰ kg, M₂ = 16 × 10²⁰ kg). In this simulation we placed the outer moon just interior to the evection resonance, so dramatic effects happen soon. However, given that a moon could be trapped at this distance by the evection resonance for a long time (Figure 14), it is fully plausible to have the inner moon converge on the outer moon while the latter is trapped by the evection resonance.

While the collision resulting from the simulation shown in Figure 15 is marginally energetic enough to disrupt the moons (Leinhardt & Stewart 2012), in reality the moons' inclinations may be higher by the time destabilization occurs. A previous residence in either the evection or 4:3 resonance, as well as crossing of any of the other evection harmonics or any other MMRs (with one another or with other satellites), would lead to inclinations of several degrees before the beginning of the last resonance, increasing the likely collision velocity.

We used direct numerical simulations to explore a past hypothetical 3:2 MMR crossing between Saturn's moons Tethys and Dione, as well as a 5:3 resonance crossing between Dione and Rhea. As the current Tethys/Dione and Dione/Rhea orbital period ratios are just above 2/3 and 3/5, respectively, and tidal recession from Saturn is expected to be fastest for Tethys and slowest for Rhea, these resonances have generally been thought to have been crossed at some point in the past. We find that the 3:2 resonance crossing always leads to excessive excitation of the orbital inclinations of both Tethys and Dione, which is inconsistent with their current orbits. Using the same approach, we conclude that the 5:3 resonance between Dione and Rhea, which closely follows the Tethys–Dione 3:2 resonance in the timeline of tidal evolution, most likely did happen. We find that the 5:3 Dione–Rhea resonance, immediately followed by a previously unknown Tethys–Dione secular resonance, is the most likely source of the current inclinations of Tethys and Rhea.

We can therefore state that Tethys and Dione evolved tidally by only a modest amount over their lifetimes, which is only about one-quarter of the tidal evolution envisaged in Murray & Dermott (1999). There are two possible interpretations: either tidal evolution of Saturn's moons has been very slow, or Saturn's midsized moons are significantly younger than the solar system. While both interpretations are consistent with the lack of the past Tethys–Dione resonance, we favor the idea that the moons are young, possibly as young as 100 Myr. Young moons would be consistent with the astrometric analysis of Lainey et al. (2012, 2015), equilibrium tidal heating of Enceladus (Porco et al. 2006; Meyer & Wisdom 2007; Howett et al. 2011), and the evolved state of the Titan–Hyperion resonance (Greenberg 1973a; Ćuk et al. 2013). Therefore, we propose that Tethys, Dione, and Rhea all formed in one large event about 100 Myr ago. The Trojan moons of Tethys and Dione that share their inclinations must have formed even more recently, after their passage through the secular resonance. Mimas, Enceladus, and the rings could have formed at the same epoch as Tethys, Dione, and Rhea, or could be younger yet. Tidal evolution of Mimas and Enceladus and their current resonances with Tethys and Dione are a complex subject (Meyer & Wisdom 2008, and references therein) and will be addressed in future work.

Our results require most craters on Rhea and the moons interior to its orbit to be produced by planetocentric impactors. Many authors have considered the cratering record of Saturn's midsized moons to result from both heliocentric and planetocentric impactors, with the latter being especially important on Mimas and Enceladus (Dones et al. 2009, p. 613; Kirchoff & Schenk 2010). Differences in crater size–frequency distributions (SFDs) between moons do not contradict our scenario, as each moon would be primarily cratered by debris on nearby orbits. Bierhaus et al. (2012, 2015) have shown that differences in surface gravity and impact velocity can produce a wide range of crater SFDs for the moons orbiting a planet. Impact speeds for planetocentric impactors onto the moons will be roughly an order of magnitude smaller than for comets from heliocentric orbit, so larger planetocentric impactors are needed to make a given crater. However, the required impactor sizes remain plausible, as they are still smaller than the crater sizes, even for basins in the Saturn system. For instance, using Keith Holsapple's tool at http://keith.aa.washington.edu/craterdata/scaling/index.htm, we estimate that a planetocentric impactor with a diameter of 40 km impacting Mimas at 3 km s⁻¹ could procuce the Herschel basin, which has a diameter of 130 km. Planetocentric impactors are also fully compatible with more than one cratering episode: the relative youth of Trojan companions and their apparent origin as debris from the larger moons imply at least one intense cratering episode following the resonance crossings described in this paper.

The absence of a past 3:2 resonance between Tethys and Dione also rules out the simplest version of the scenario proposed by Charnoz et al. (2011) and Crida & Charnoz (2012), in which the moons formed sequentially out of past massive rings and then migrated outward. That hypothesis implies that Tethys had to "catch up" with Dione, invariably crossing the 3:2 resonance. Our results do not exclude the possibility that the previous generation of moons did form this way but were destroyed and re-accreted close to their present positions. Other ultimate origins of the midsized satellite system, such as those proposed by Canup (2010) and Asphaug & Reufer (2013), or a late origin of Titan (Hamilton 2013), can certainly be reconciled with a recent episode of re-accretion. Disruptions, followed by formation of a debris disk and re-accretion, erase much of the diagnostic information about the moons' formation and put few constraints on the ultimate origin of the system.

Using order-of-magnitude arguments, we conclude that the disk that gave rise to the present midsized moons can only be short-lived and must be a product of disruption of a previous generation of midsized moons. We explored different configurations of Saturn's past satellite system in order to find a likely source of such an instability. The most likely trigger for a global instability is the evection (semisecular) resonance between the largest past moon (assumed to be Rhea sized) and the Sun, located just interior to the present orbit of Rhea. If most of the current icy moons' mass was concentrated in two moons on converging orbits (somewhat similar to Dione and Rhea), with the outer one evolving into the evection resonance, we find that a large-scale instability and mutual collisions are a common outcome. We also propose disk-driven migration of small moons (at the edge of the resulting disk) into orbital resonances with Titan as a plausible source of Titan's eccentricity.

There are three possible avenues toward testing the hypothesis of the young Saturnian midsized satellite system. While preliminary results appear to support fast migration (Lainey et al. 2015), further work on the satellite astrometry using Cassini data should be able to decisively confirm or falsify the findings of Lainey et al. (2012). Second, better understanding of the geophysics of Enceladus may be able to clarify whether its energy source is indeed tidal and whether the dynamical equilibrium hypothesis is appropriate. Third, this hypothesis requires that the majority of craters on Saturn's moons interior to Titan (many of which have heavily cratered surfaces) must be the product of planetocentric impacts, as there would not have been enough time for them to accumulate sufficient numbers of cometary impacts (Dones et al. 2009, p. 613). We hope that further work on crater and impactor population modeling will be able to put stronger constraints on impactor sources.

M.Ć. acknowledges support from NASA's Outer Planets Research Program, awards NNX11AM48G and NNX14AOO38G. We thank an anonymous reviewer for very insightful comments that greatly improved the paper.