A Past Episode of Rapid Tidal Evolution of Enceladus?

Janus and Epimetheus are inner moons of Saturn that are currently in a horseshoe orbital configuration which undergoes a reversal every four years (Murray & Dermott 1999). Gravitational interactions with Saturn's rings make the width of the "horseshoe" decrease over time, and the timescale for this evolution is on the order of 10⁷ yr (Lissauer 1985). For our purposes Janus and Epimetheus (and their presumed parent body) are interesting because they are expected to migrate very fast due to ring torques, potentially crossing resonances with other moons. Tajeddine et al. (2017) estimate that Janus migrates about 40 km Myr⁻¹ due to ring torques, and the two moons' putative parent body (if they originated in a breakup) would have migrated even faster, at 50 km Myr⁻¹. Janus and Epimetheus are about 1500 km exterior to the 2:1 MMR with Enceladus. Ignoring the tidal evolution of Enceladus for the moment, it appears that the precursor of Janus and Epimetheus should have crossed this resonance about 30 Myr ago, putting constraints on the age and evolution of the system.

We simulated this resonance crossing between coorbitals and Enceladus using the numerical integrator simpl, previously used by Ćuk et al. (2016). Enceladus was assumed to be in its present orbital resonance with Dione. At first we assumed that Janus and Epimetheus were in a horseshoe configuration before the resonance, and gave them their present eccentricities and inclinations. The left-hand panels in Figure 1 show a typical outcome of such simulations, in which the horseshoe configuration is disrupted and Janus and Epimetheus experience close encounters (our integrator does not test for collisions). This destabilization happens well before the core of the 2:1 resonance with Enceladus is reached, implying that the horseshoe is disrupted by near-resonant perturbations.

Zoom In Zoom Out Reset image size

Next we assumed that the two moons were in corbitals before the resonance, but in a Trojan or "tadpole" configuration (Murray & Dermott 1999); typical results are shown in the right-hand side panels of Figure 1. Interestingly, near-resonant perturbations always convert the tadpole configuration into a horseshoe, which is stable from there on. However, the moons are then captured into a 2:1 resonance with Enceladus, in which their eccentricities grew over time while the semimajor axis is locked to that of Enceladus. This eccentricity growth is largely unaffected by any tidal damping, as the small sizes and demonstrably rigid natures of the irregularly shaped Janus and Epimetheus do not allow for significant tidal deformation or dissipation over the 10⁷ yr timescales in which we are interested. We expect the resonance eventually to break and Janus and Epimetheus to collide. As we expect the moons to have large eccentricities at the time of the collision while their semimajor axes are initially within resonance with Enceladus, the resulting debris must reaccrete significantly interior to the resonance (due to angular momentum conservation).

We also ran a simulation assuming that Janus and Epimetheus were contained within one body during the resonance crossing with Enceladus (which may have subsequently broken up). Our simulations of the 2:1 resonance crossing with Enceladus show that the precursor is always captured in the resonance (Figure 2). Note that we assumed the current orbital eccentricity of Janus for the progenitor, and the present-day configuration of the Enceladus–Dione resonance. In most simulations, proto-Janus is captured into a stable corotation resonance in which the eccentricity of Enceladus is slightly higher than equilibrium eccentricity in the resonance with Dione (left-hand panels in Figure 2). In some cases this resonance is broken and proto-Janus is captured into the e_E (i.e., Lindblad) subresonance in which the eccentricity of proto-Janus grows over time, like we saw for the initial tadpole configuration. Given the large orbital precession rates of the innermost moons, all resonances are well separated and there is no obvious way of breaking this lock purely through dynamics. However at eccentricities in the 0.05–0.1 range proto-Janus would cross the orbits of Prometheus and Pandora (if they were present at that epoch) or the outer edge of the rings. Resulting collisions would presumably break the resonant lock between proto-Janus and Enceladus. Once again, angular momentum conservation dictates that any reaccretion of proto-Janus's debris happens interior to the resonance with Enceladus.

Zoom In Zoom Out Reset image size

If we assume that the coorbitals were produced in the breakup of a single progenitor, we can use their orbits to constrain the orbit of this progenitor. Using the current orbits of the coorbitals, and under the assumption of conservation of momentum at separation, we find that the prebreakup parent body of Janus and Enceladus had an eccentricity e ≤ 0.006, much less than expected from the past capture into the 2:1 resonance with Enceladus. Eccentricity can be erased through reaccretion, but as we discuss above, disruption and reaccretion move the resulting new moons interior to the resonance with Enceladus. Therefore, either the tidal evolution of Enceladus has to be comparable with the orbital evolution of the horseshoe pair, or the inner Saturnian system has to be rather young (a combination of these two factors is also possible).

If the resonant pair Enceladus and Dione are evolving through equilibrium tides with Q/k₂ = 5000 for Saturn (Lainey et al. 2012), their joint orbital evolution timescale would be $a/\dot{a}\approx 10\,\mathrm{Gyr}$. Enceladus being locked to an inertial mode would produce the same timescale (even if the amount of tidal heating would be very different, as all of Dione's orbital evolution would now be due to Enceladus). This rate of evolution would only push the resonance with the Janus–Epimetheus precursor back to 40 Myr. Locking of Enceladus to a Saturn-frame mode would lead to $a/\dot{a}\simeq 130\,\mathrm{Gyr}$, which would produce negligible migration compared to the ring-torque-driven evolution of the coorbitals' precursor, keeping the apparent age of the resonant encounter at 30 Myr ago. Enceladus's rate of evolution through equilibrium tides (but not resonance locking) depends on its MMR with Dione. Before the resonance Enceladus would have had an equilibrium tidal evolution timescale $a/\dot{a}=5\,\mathrm{Gyr}$ (assuming Q/k₂ = 5000 for Saturn), so depending on the age of the Enceladus–Dione resonance, the progenitor–Enceladus resonance would have happened between 40 and 75 Myr ago.

Lainey et al. (2020) propose that a model in which equilibrium tides are weak and the moons' orbits evolve solely through resonance locking may result in a long-term stable system. In this model, the moons' orbits are not expected to converge and enter MMRs with each other. In this way dynamical excitation and possible destabilization through MMRs are avoided. The simulations discussed here show that the inclusion of ring torques introduces relatively rapid convergence of the orbits of some of the inner moons. A young age (<50 Myr), at least for the Janus–Epimetheus parent body (JEPB) if not for other moons, appears inevitable. If we assume the presence of both equilibrium tides and resonant tidal response (Section 5.3) the actual maximum age for the JEPB may be somewhat older than 75 Myr but difficult to calculate precisely. Also, it appears that the JEPB formed close to its present location, and that the idea that Janus and Epimetheus evolved from the rings smoothly to their present positions (Crida & Charnoz 2012) may not be tenable given the apparent impossibility of the pair (or their parent body) crossing the 2:1 resonance with Enceladus.

Mimas and Tethys are currently in a 4:2 inclination-type resonance that involves the inclinations of both moons (the resonant argument is 4λ_Θ − 2λ_M − Ω_M − Ω_Θ), where the subscripts M and Θ refer to Mimas and Tethys, respectively. A large libration amplitude of the resonance has in the past been used to suggest that the resonance capture was a low-probability event (Allan 1969; Sinclair 1972). Our work so far (under the assumption of equilibrium tides) suggests that there are at least two distinct dynamical pathways to resonance capture, both of which require a later increase of the resonance argument libration width due to a third body.

If the initial inclination of Mimas was very low (i_M ≃ 0001), the lack of capture into the preceding i_M ² subresonance and capture into i_M i_Θ subresonance are highly likely events (assuming i_Θ ≈ 1° before the resonance). The right-hand panels in Figure 3 shows a typical capture into the Mimas–Tethys resonance assuming equilibrium tides with Q/k₂ = 4000 for Saturn. This rate of tidal evolution implies that the current Mimas–Tethys resonance is only 20–25 Myr old, consistent with the relatively young proposed age of the system (Ćuk et al. 2016). The physical reason behind the lack of capture into the i_M ² subresonance appears to be the nonadiabatic nature of the resonance crossing, as proposed by Sinclair (1972). While the tidal evolution is smooth and slow by most criteria, the extremely narrow width of the i_M ² subresonance at very low inclinations seems to allow for nonresonant crossing.

Zoom In Zoom Out Reset image size

Alternatively, if we assume a preresonant inclination of i_M = 01 for Mimas, there is a high probability (about 70%–80%) of missing capture into the i_M ² subresonance, and a similarly high probability of capture into the i_M i_Θ subresonance. In this case the resulting libration amplitude of the resonant argument is about 60°, still short of the observed one but much larger than in the low-inclination limit. In this regime the capture into i_M ² is unlikely because of the large initial inclination (as is common for resonant encounters; Murray & Dermott 1999), while the capture into i_M i_Θ subresonance is somewhat more likely due to large i_Θ (also, a negative kick to the i_M when crossing i_M ² also makes the next capture more likely).

Figure 3 shows two of the many simulations we did to study capture into the 4:2 Mimas–Tethys resonance in isolation, meaning that we did not simulate the system's evolution before and after this encounter. Initially we expected the two subresonances of the 4:2 MMR (i² _M and i_M i_Θ) to be the only relevant inclination-type terms of the 2:1 Mimas–Tethys commensurability. Later simulations of the comprehensive recent dynamical history of the system have revealed additional relevant inclination-type resonances. In particular, a third-order harmonic of the 2:1 Mimas–Tethys resonance with the argument 2λ_Θ − λ_M + ϖ_M − Ω_M − Ω_Θ is surprisingly strong, and can lead to capture for low initial inclinations of Mimas (see Figure 4). This subresonance (which we designate e_M i_M i_Θ) is encountered well before the i_M ² and i_M i_Θ subresonances shown in Figure 3. The reason for the unexpected strength of the e_M i_M i_Θ resonant term is the high inclination of Tethys and high eccentricity of Mimas, which we assumed both predate these moons' 4:2 MMR. These assumptions were based on a lack of available mechanisms of exciting e_M after the capture into the current 4:2 Mimas–Tethys resonance, and the necessity of large i_Θ for the capture to occur. We find that an initial i_M ≈ 01 is necessary in order to avoid capture into the e_M i_M i_Θ subresonance. Therefore, while both low- and high-inclination routes to capture into the observed i_M i_Θ resonance are possible in isolation, additional constraints from the e_M i_M i_Θ crossing make the high-inclination path the only likely capture mechanism.

Zoom In Zoom Out Reset image size

The captures into resonance shown in Figure 3 result in a libration amplitude of the resonant argument that is smaller than observed. In general, we find that that the libration amplitude of the Mimas–Tethys resonance is highly vulnerable to perturbations from other moons (including Enceladus, Section 4.2), so this resonance's libration amplitude is less indicative of the capture mechanism than previously thought. We note that the difference between the postcapture libration amplitude and the observed value is less pronounced in our preferred high-inclination route to Mimas–Tethys 4:2 MMR compared to the low-inclination route. Clearly more detailed study of the system's very recent evolution, accounting for the influence of all the moons (including Tethys's small Trojan companions), is necessary to understand better the evolution of the libration amplitude of the Mimas–Tethys resonance.

All hypotheses of the origin of the Mimas–Tethys resonance must include a convergent evolution of these two moons. This convergent evolution is expected in the case of equilibrium tides, but interestingly only in the "constant Q" rather than the "constant time-lag" model (as the synodic period of Mimas is relatively long). In the case of evolution by both moons being locked to inertial waves there should in principle be no relative orbital evolution, either convergent or divergent, while resonance locking to Saturn-frame modes would result in divergent evolution. If only Mimas is in resonant lock, then a resonance would form and should be about 20 Myr old (in the case of inertial waves) or about 300 Myr old (in case of locking to Saturn-frame modes). Locking of Mimas to a Saturn-frame mode would require a "background" Q ≥ 100,000, as faster equilibrium tidal migration would allow Mimas to "outrun" the mode.

Another common assumption in all scenarios on the origin of the Mimas–Tethys resonance requires Tethys to have a prior inclination of about i_Θ ≈ 1°. This relatively large inclination requires some kind of past dynamical interaction with other moons, most likely a resonance. Ćuk et al. (2016) have found this inclination to be a plausible consequence of Dione–Rhea 5:3 resonance crossing, followed by the (associated) Tethys–Dione secular resonance. We note that this scenario requires convergent orbital evolution of Dione and Rhea, at least in the past (see Section 3.3).

The orbit of Rhea is directly observed to evolve on a 6 Gyr timescale (Lainey et al. 2017, 2020; Jacobson 2022). The rate of orbital evolution of Titan may be similarly fast (Lainey et al. 2020), but this is still disputed (Jacobson 2022). Observational results for Rhea are clearly not consistent with the expectations based on equilibrium tides. Locking to resonant modes within Saturn has been proposed as an explanation for these observations (Fuller et al. 2016), with the apparent rate of the moon's orbital evolution being determined by the drift in frequency of the resonant peaks in the tidal response of the planet. After Rhea's fast evolution was first observed (Lainey et al. 2017), the assumption that the peaks in the tidal response have a constant absolute frequency drift in Saturn's rotating frame (Fuller et al. 2016) implied that Titan's own drift timescale should be $a/\dot{a}\lt 2\,\mathrm{Gyr}$. The observed tidal evolution rate of Titan of 11 Gyr (Lainey et al. 2020) would imply that the moons could not both be in resonant lock, or that the resonant mode drift is not constant in the rotating frame. Locking to inertial waves which drift together in an inertial reference frame was proposed instead (Lainey et al. 2020), which should result in the same tidal evolution timescale for all moons in resonant lock. The difference between the apparent evolution timescales for Rhea and Titan reported by Lainey et al. (2020) seems to suggest that moons under resonant lock can still have converging orbits, which would challenge the big picture of Lainey et al. (2020), in which moons never cross mutual resonances. On the other hand, if Titan has a much slower orbital evolution consistent with equilibrium tides (as found by Jacobson 2022), then we can be sure that not all moons are locked to resonant modes. Clearly the dynamics of the system is more complex than predicted in any one simple model.

The observed fast evolution of Rhea is at odds with the hypothesis that Dione and Rhea crossed their mutual 5:3 MMR. The Dione–Rhea 5:3 MMR crossing was proposed by Ćuk et al. (2016) as a way of producing the current inclinations of Tethys and Rhea. The inclination of Tethys is particularly significant (i_Θ ≈ 1°) and must have predated Tethys's current resonance with Mimas. As Ćuk et al. (2016) show, the Dione–Rhea 5:3 resonance is usually followed by a secular resonance between Tethys and Dione during which Dione "passes" its eccentricity and inclination to Tethys. Therefore there is a compelling reason to think that a past Dione–Rhea 5:3 MMR crossing did take place.

The orbital evolution of Rhea therefore presents a conundrum. Its rate is much too fast for equilibrium tides (assuming Q/k₂ = 5000), but it is also too fast for resonance locking (assuming that Titan is currently in resonance lock; Lainey et al. 2020). A Saturn-frame mode resonant lock would have $a/\dot{a}\approx 40\,\mathrm{Gyr}$, while a resonance lock with inertial modes would predict $a/\dot{a}=11\,\mathrm{Gyr}$, same as Titan. Furthermore, the apparent past crossing of the Dione–Rhea resonance implies that this fast evolution of Rhea is a very recent phenomenon (i.e., over the past several million years or so). Assuming all of these constraints are correct, the simple solution is that Rhea is not locked to a mode, but is currently crossing a resonant mode. This is a possible situation if a moon (through equilibrium tides) evolves faster than a mode rather than the other way around, making a resonant lock impossible but producing occasional bursts of fast orbital evolution. If we assume that Titan is currently evolving with a timescale of 100 Gyr (consistent with equilibrium tides; Jacobson 2022), then it is possible that resonance locking may not be feasible at all. The argument here is that Titan would have certainly been captured in a resonant lock if resonant modes moving faster than Titan's orbit evolve due to equilibrium tides. In principle, we can also envision resonant modes moving to higher frequencies over time, opposite the direction of tidal evolution. In this case, all accelerated tidal evolution would always be episodic, being driven by the moons overtaking resonant modes.

In the rest of the paper, we will examine how a combination of equilibrium tides and passage through resonant modes can explain other features of the Saturnian system, namely the establishment of the Enceladus–Dione MMR.

The Enceladus–Dione 2:1 MMR contains two first-order subresonances, four more second-order ones, several important three-body resonances, as well as numerous two-body third-order subresonances. The currently occupied e_E subresonance is one of the last to be encountered as the orbits of Enceladus and Dione converge. We ran a number of simulations of the Enceladus–Dione 2:1 MMR encounter assuming equilibrium tides and Q/k₂ = 4000 for Saturn. Typically, Enceladus was captured into the 4:2 i_E ² resonances in a large majority of cases where Enceladus and Dione encounter this resonance under equilibrium tides. The left-hand panels in Figure 5 plot the result of a typical simulation that features capture into the Enceladus–Dione 4:2 i_E ² resonance. Enceladus's inclination grows to about 1°, at which point secondary resonances break the inclination resonance. In this particular simulation all inner moons (from Mimas to Rhea) were assumed to have initially equatorial orbits, but the result is the same when Tethys was assumed to have been already inclined. In all remaining simulations of equilibrium tide Enceladus–Dione 4:2 MMR crossing Enceladus is caught into a more complex inclination resonance, also affecting the inclination of Dione or Tethys. In all cases Enceladus acquires an inclination on the order of a degree, two orders of magnitude above the observed one.

Zoom In Zoom Out Reset image size

The near certainty of capture into one of the inclination subresonances found in our simulations strongly suggests that the Enceladus–Dione 2:1 MMR was not assembled through equilibrium tides. This is in contrast to what we found for the Mimas–Tethys 4:2 resonance (Section 3), which can be accounted for using equilibrium tides (except for the present libration amplitude, which is easily excited by a third body). The high probability of capture into the Enceladus–Dione 4:2 i_E ² resonance, but not Mimas–Tethys 4:2 i_M ² resonance is rather surprising but holds in large numbers of numerical simulations. We suspect that the lack of capture into the Mimas–Tethys 4:2 i_M ² resonance with a low initial i_M is due to the mechanism originally proposed by Sinclair (1972), which states that in cases of small free inclinations and fast orbital precession even a relatively slow resonant encounter may not be adiabatic, i.e., a second-order resonance may not be able to dominate the motion of the ascending node against oblateness-driven precession. Note that this explanation did not hold in numerical simulations using Q ≈ 18,000 for Saturn (Sinclair 1974), but appears to explain the behavior at the ten times faster tidal evolution used in this study.

In the course of our work on the past dynamics of the Saturnian system, we discovered that Enceladus–Tethys 11:8 resonance offers unexpected constraints on the past orbital evolution of those two moons. In the equilibrium tide paradigm, this resonance would have been crossed 10–15 Myr ago. Since Tethys's predicted equilibrium tidal evolution is faster than that of Enceladus (they are the only neighboring pair of major Saturnian moons on diverging orbits, as Tethys's six times larger mass boosts its tidal evolution), this resonance should have been crossed divergently without a chance of capture. Despite the third order of this subresonance, the large eccentricity of Enceladus and inclination of Tethys make this resonant term surprisingly relevant. In particular, we find that the Enceladus–Tethys 11:8 MMR crossing usually both excited the inclination of Enceladus beyond the observed value, and disrupted the Mimas–Tethys 4:2 MMR (Figure 6).

Zoom In Zoom Out Reset image size

These results decisively argue that the current system could not have resulted by evolution through equilibrium tides. Of course, observations of the tidal acceleration on Rhea (Lainey et al. 2017, 2020; Jacobson 2022) already established that Saturn's tidal response is dynamic, and that Rhea and/or Titan may be in a resonance lock. However, as we demonstrated in previous sections, the observed resonances between the inner moons appear to suggest convergent tidal evolution consistent with equilibrium tides. The apparent lack of past crossing of the Enceladus–Tethys 11:8 MMR is the first evidence from mutual resonances that requires frequency-dependent tides to have acted in the past.

In general, the Enceladus–Tethys 11:8 MMR could have been avoided through slower orbital evolution of Tethys or faster migration of Enceladus. We explored the former possibility assuming that the tidal response of Saturn is frequency dependent but does not contain resonant modes. This would be consistent with the findings of Lainey et al. (2020) and Jacobson (2022), who do find a slower orbital evolution of Tethys relative to the other inner moons than expected from equilibrium tides (however, the uncertainties are still large). The current tidal evolution rate of Enceladus is constrained by its tidal heating (Section 2), so a slower evolution of Tethys is necessary to avoid a recent passage through the Enceladus–Tethys 11:8 MMR. However, we have just seen that the equilibrium tides acting on Enceladus are still unable to reproduce the Enceladus–Dione 2:1 MMR (Section 4.1). Therefore it appears that the most likely solutions to the problem posed by the Enceladus–Tethys 11:8 MMR require both that (1) the tidal response at Tethys's frequency is currently weaker than that of Enceladus, and (2) that the migration of Enceladus was much faster in the past. We explore this possibility in the next section.

It is a well-known phenomenon that the probability of MMR capture declines as the resonance is encountered more rapidly (Murray & Dermott 1999). In this subsection we explore the hypothesis that an episode of fast tidal evolution may have enabled Enceladus to cross the inclination resonances with Dione without capture. The results from a number of simulations are plotted in Figure 7. We ran 20 different simulations that had Saturn's Q/k₂ = 4000 for all moons except Enceladus, for which the response was enhanced five times; six resulted in capture into an inclination-type resonance and 14 experienced a kick in inclination. In 20 simulations where Enceladus evolved 10 times faster than expected, none experienced capture. The postpassage inclination of Enceladus was on average lower in the 10× accelerated simulations (typically i_E = 001−002) than in the 5× accelerated ones (typically i_E = 002−004).

Zoom In Zoom Out Reset image size

Our choice of current inclinations and eccentricities for the initial conditions is somewhat arbitrary. While its present inclination is small, Enceladus presumably formed with close to zero inclination and this initial value may affect the postresonant distribution of inclinations. Later simulations (shown in Section 5.3) were started with a lower inclination of Enceladus but the outcomes were very similar. As seen in Figure 7, the three-body resonance with Tethys (middle of the three jumps) is a minor contributor to Enceladus's final inclination, so our assumption that Tethys was already inclined at this epoch is not crucial. The exact inclination of Dione may effect the results somewhat (as the third jump is the i_E i_D subresonance), but we see no strong reason to expect a different inclination for Dione at the time when 2:1 resonance with Enceladus was established.

After the inclination-type subresonances have been crossed, Enceladus encounters a number of eccentricity-type subresonances of the Enceladus–Dione 2:1 MMR. Based on the results shown in Figure 7, we enhanced Enceladus's tidal evolution by a factor of 10, while the tidal evolutions of the other moons correspond to the frequency-independent Q/k₂ = 4000 for Saturn. The choice of initial eccentricities requires making assumptions not only about past resonances but also about the rate of eccentricity damping by the moons.

In our first batch of simulations we used the current eccentricities for most moons, except that Enceladus was given a low initial e_E < 0.001 while Dione was given e_D = 0.004, about twice the current value (to account for the expected resonance kick and subsequent damping). Out of the first group of simulations using these initial conditions, the majority of runs resulted in Enceladus being captured into three-body resonances with Rhea or Tethys (Figure 8, first and third rows), and the rest resulted in the e_D subresonance, with no simulations reaching the current e_E subresonance (which is last in order of encounters).

Zoom In Zoom Out Reset image size

Captures into three-body resonances, which were understandably missed in past semianalytical models (Meyer & Wisdom 2008) demonstrate once again the importance of such resonances in the Saturnian system (Ćuk & El Moutamid 2022). These captures may come as a surprise given the high rate of orbital migration of Enceladus, but we should recognize that these are first-order resonances, unlike the inclination-type resonances shown in Figure 7, which are necessarily of the second order. This preponderance of three-body resonance captures is in part a consequence of us starting the simulations with Tethys and Rhea on very-low-eccentricity orbits. However, these were not arbitrary choices, as these two moons currently do have almost circular orbits. The current low eccentricities of both Tethys and Rhea argue against these three-body resonance captures happening in the recent past, as damping the eccentricities acquired in Figure 8 over realistic timescales may be difficult, especially for Rhea. We therefore prefer these moons having relatively small eccentricities at this epoch and avoiding capture into the three-body resonances, rather than being captured and acquiring large eccentricities (on the order of 0.01 or more).

In the second batch of simulations, we gave both Tethys and Rhea eccentricities of 0.002, which we found (through trial and error) to be largely sufficient to avoid capture into three-body resonances. In the second group of simulations, Enceladus still never reaches the current e_E resonance, but in the majority of cases becomes captured in a second-order resonance which also includes the eccentricity of either Titan or Dione (Figure 8, rows four and five, respectively). While the mixed resonance e_E e_D has been studied before (Meyer & Wisdom 2008), the three-body resonance with Titan is new. In both of these cases we cannot avoid capture into these subresonances by changing the eccentricity of Dione or Rhea. A less-eccentric Dione would make capture into e_E e_D resonance less likely, but would increase the fraction of captures into the first-order e_D subresonance (this fraction is already 20%–30% for initial e_D = 0.004) and may conflict with Dione's current eccentricity. In the case of Titan, its eccentricity does not get affected by the resonance and we must use values very close to the current one.

The simplest way to avoid capture into the abovementioned second-order eccentricity subresonances is to give Enceladus a sizable initial eccentricity. The implication is that at the time of capture into the Enceladus–Dione 2:1 MMR, not only was an accelerated orbital evolution of Enceladus necessary (to avoid capture into the inclination-type resonances), but all inner moons had to have somewhat eccentric orbits (to avoid capture into the eccentricity-type resonances). In the next section we will put together a possible evolution path that satisfies the constraints presented in this and previous sections, and the implications of the moons' previously excited eccentricities will be discussed in Section 6.

In this section we attempt to construct a numerical model of the recent (last 20 Myr) evolution of the Saturnian satellite system. We are using all of the constraints that were established in previous sections. Enceladus and Dione are assumed to experience tidal acceleration with Q/k₂ = 4000 at the present epoch, ³ in agreement with the hypothesis of equilibrium tidal heating of Enceladus (Section 2; Meyer & Wisdom 2007; Howett et al. 2011; Lainey et al. 2012). Mimas is assumed to have the same tidal evolution rate, both because of a lack of other constraints, and as it is consistent with Mimas encountering and being captured into 4:2 resonance with Tethys (Section 3.2). Need to avoid the 11:8 Enceladus–Tethys resonance, as well as the (admittedly uncertain) results of Lainey et al. (2020) and Jacobson (2022), suggested that the tidal acceleration of Tethys is much weaker than that expected from equilibrium tides. Finally, and most radically, we assumed that Enceladus went through a stage of very rapid tidal evolution due to crossing of (but not locking to) a resonant mode. This was done to avoid capture into Enceladus–Dione 4:2 inclination resonances (Section 4.1), as well as to avoid a recent crossing of Enceladus's 1:2 resonance with the horseshoe moons (or their progenitor; Section 3.1). Additional assumptions include low initial tidal dissipation within Enceladus until it experiences strong tidal heating, at which point Enceladus switches to a strong tidal response (this is done to approximate melting of the interior). Due to the nature of the integrator where the tidal response is treated as a constant parameter, we modeled both the transitions in the tidal acceleration of Enceladus and the changes in its own tidal response as abrupt events, implemented manually.

Table 1 shows the initial conditions for two sets of simulations of the integrated evolution of the inner moons (sets differ only by initial i_Θ). In both cases, Enceladus was assumed to be experiencing 10× accelerated evolution in the first 10 Myr of the simulation. After 3 Myr we selected one simulation in each set in which Enceladus was captured into the 2:1 e_E MMR with Dione and acquired a high eccentricity, cloned those simulations, and from that time onward used a much higher tidal response of Enceladus. At 10 Myr we reverted Enceladus to an equilibrium tidal evolution rate and current tidal parameters, evolving the system until 20 Myr or until the Mimas–Tethys resonance reached its present state. We also assumed equilibrium tides on Rhea before 10 Myr and accelerated evolution after than (Section 3.3), and accelerated Titan throughout, but this had little importance as nether Rhea not Titan encountered any resonances.

Figure 9 shows the outcomes of two successful simulations (one from each set), which come close to reproducing the current system. Between the two sets of first-stage simulations (20 total), 10 were "successful," meaning that Enceladus and Dione were in their 2:1 e_E MMR at 3 Myr. The other 10 simulations were caught in one of the other subresonances of the Enceladus–Dione 2:1 MMR (Figure 8), typically exciting the eccentricities of Tethys or Dione well above the observed values, or "fell out" of the resonance due apparent interaction with secondary resonances. We chose two of the successful simulations for cloning, which was done by sharply changing the tidal properties of Enceladus. The tidal Love number was changed from k_2E = 0.01 to k_2E = 1, approximating melting. Different clone simulations had the tidal Q of Enceladus in the 30–39 range. This range was chosen through trial and error, as we found that simulations with A = (Q_S k_2E )/(Q_E k_2S ) < 10 led to breaking of the resonance, while those with A > 10 settled into bound cycle within the 2:1 e_E subresonance (subscript S designates Saturn's properties at Enceladus's frequency). This critical change in the dynamics of resonance with the change in Enceladus's tidal response was first discovered by Meyer & Wisdom (2008), and our value for critical A agrees with theirs. ⁴ This phase of the evolution of Enceladus, despite having seemingly chaotic oscillations, is not stochastic but determined by the tidal properties, so all simulations in which Enceladus is dissipative enough (i.e., A > 10) will stay at the threshold of the Enceladus–Dione 2:1 e_E resonance indefinitely.

Zoom In Zoom Out Reset image size

We assume that Enceladus eventually gets out of the resonant mode and migrates due to the "background" tidal response of Saturn that is not dependent on frequency. In reality this transition would be gradual as resonant modes are expected produce a smooth (if narrow) profile of tidal response over frequency (Fuller et al. 2016). We could not simulate this without extensively modifying the integrator and, more importantly, greatly increasing the complexity of our model. At 10 Myr all of our clone simulations switch to Enceladus experiencing Saturn's tides with Q/k₂ = 4000 and having its own k₂/Q = 0.01, approximately what is expected if the current heating is equal to the observed and is in equilibrium. Out of 20 cloned simulations only three do not complete the transition to a narrow libration within the resonance consistent with observations (in two simulations Enceladus leaves the resonance, in one it is permanently trapped in a large libration state due to secondary resonances). In some simulations the transition is smooth like in the right-hand side of Figure 9, while in the others secondary resonances are encountered during the evolution, leading to temporary increases in the libration amplitude (as in the left-hand side of Figure 9).

The migrations of Mimas, Tethys, and Rhea are largely unaffected by the Enceladus–Dione resonance in our simulations. There is some variation of the eccentricity of Mimas when Enceladus encounters secondary resonances within the 2:1 e_E MMR with Dione, but the effect on Mimas is usually minor. The inclination of Mimas is unaffected by any of the dynamics discussed so far, and passes through a series of kicks due to subresonances of the Mimas–Tethys 2:1 commensurability. Subresonances with arguments 2λ_Θ − λ_M + ϖ_M − 2Ω_M , 2λ_Θ − λ_M + ϖ_M − Ω_M − Ω_Θ (see Figure 4), and 4λ_Θ − 2λ_M − 2Ω_M (the pure i² _M term, Figure 3) are encountered at about 3 Myr, 7 Myr, and 14 Myr, respectively. Since the latter two kicks (which are larger) always happen after the simulations were cloned at 3 Myr, we can consider the clones of the original two simulations to be practically independent when it comes to the inclination of Mimas. Out of 18 clones which stayed in the Enceladus–Dione 2:1 MMR, nine were captured into the current Mimas–Tethys 4:2 i_M i_Θ subresonance (bottom panels in Figure 9). The remaining clones either experienced capture into the i_M ² harmonic or passed the whole forest of subresonances without capture. Therefore we can say that both observed resonances (Mimas–Tethys 4:2 i_M i_Θ and Enceladus–Dione 2:1 e_E ) have high probability (about 50%) when using these initial conditions and assumptions about Enceladus's tidal evolution and internal dissipation. While we can quantify the probabilities of different stochastic outcomes of resonant dynamics, we are currently not able to assess the a priori probability of our initial conditions and assumptions about Saturn's tidal response. Note that our preferred timeline requires Enceladus to settle into the present quiescent state before Mimas acquires a high inclination through its resonance with Dione. The current Mimas–Tethys resonance is very fragile and we often see it broken if it coincides with Enceladus having a very high and/or chaotic eccentricity.

The other parameters of the system at the ends of the simulations in Figure 9 are mostly consistent with the present state. The Mimas–Tethys resonance in every case has a large libration amplitude that is still well short of the observed one (93°), but we expect even very minor subsequent perturbations to modify this quantity (we ignored all moons smaller than Mimas in these runs). The inclinations of Enceladus and Dione are somewhat stochastic but in the range that includes the observed values (i_E = 0008 and i_D = 002, respectively), as is the eccentricity of Dione (probably the most variable quantity in our simulations; currently e_D = 0.002). For both Mimas and Tethys the inclinations are determined by the initial conditions and their mutual resonance, and the eccentricities depend primarily on the initial values and tidal dissipation, without much stochasticity. However, the eccentricity of Rhea is one quantity that our simulation cannot explain, as we have set it much higher than the observed value in order to avoid the three-body argument of the Enceladus–Dione MMR (Figure 8), but for the tidal dissipation we assume for Rhea cannot subsequently modify that eccentricity. This discrepancy may tell us something about the recent dynamics of Rhea, as we discuss in the next section.

Much of our reasoning that resulted in the model of recent evolution presented in Figure 9 is driven by the survival of the low inclination of Enceladus. In order for Enceladus's inclination to stay so low, both the capture into the 4:2 Enceladus–Dione MMR (Sections 4.1 and 5.1) and crossing of the 11:8 Enceladus–Tethys MMR (Section 4.2) must be avoided. Can we be confident that the inclination of Enceladus was not actually higher in the past and then damped by tides?

The simplest way to estimate the tidal damping of inclination within Enceladus is assuming a homogeneous Enceladus that can be described by a single tidal Love number k₂ and a tidal quality factor Q. If we assume that the eccentricity of Enceladus is currently in equilibrium between excitation by resonance with Dione and damping by tides which produce the observed heating (Howett et al. 2011), we get k₂/Q ≈ 0.01 (Meyer & Wisdom 2007; Lainey et al. 2012). This corresponds to an eccentricity-damping timescale of about 0.5 Myr, making Enceladus exceptionally dissipative. The timescale for damping of inclination is longer by the factor $7{(\sin i/\sin \theta )}^{2}$, where θ is the forced obliquity of the moon. The forced obliquity of Enceladus has been modeled by Baland et al. (2016) and they found that θ ≤ 4 × 10⁻⁴°. This would make $(\sin i/\sin \theta )\geqslant 20$, making the timescale for damping of Enceladus's inclination over 1 Gyr, clearly too slow to affect the dynamics in the time frame we consider.

Another mechanism for inclination damping would be the resonant response of a global ocean (Tyler 2008, 2011). Given that Enceladus likely possesses a subsurface ocean, this possibility must be addressed. Chen et al. (2014) analyzed in depth different mechanisms of tidal dissipation, and found that at Enceladus's present obliquity the heating due to obliquity tides in the ocean is more than three orders of magnitude lower than nonresonant obliquity tides due to the whole-body response of Enceladus. This also implies that the inclination-damping timescale due to resonant ocean tides is longer than 10³ Gyr. One important finding by Chen et al. (2014) about resonant tides is that their power depends as a cube of obliquity (as opposed to the square of obliquity in nonresonant tides), meaning that the timescale for inclination damping is inversely proportional to obliquity. However, even if the inclination of Enceladus were about a degree (see Figure 7) and its forced obliquity were θ ≈ 005, resonant-inclination damping would still be an order of magnitude slower than nonresonant obliquity tides. Finally, even if the parameters used by Chen et al. (2014) were not to apply (for some unforeseen reason) to Enceladus on a high-inclination orbit, their calculations clearly demonstrate that the inclination of Enceladus cannot be brought to its present low value by resonant ocean tides.

The above discussion assumed that Enceladus is in Cassini state 1, in which the forced obliquity due to inclinations is low, and that dissipation within other moons does not affect the evolution of Enceladus's inclination. However, it is in principle possible that the obliquity of a moon could be excited by a spin–orbit resonance. This would be equivalent to the excitation of obliquity of Saturn itself by a secular resonance between the axial precession of Saturn and nodal precession of Neptune's orbit (Hamilton & Ward 2004; Ward & Hamilton 2004). Recently Ćuk et al. (2020) suggested that Uranian moon Oberon may have in the past been caught in a similar resonance with the orbit of Umbriel, and that tidal dissipation within Oberon may have damped the inclination of Umbriel. Using the results of Chen et al. (2014), we surveyed the precession frequencies in the Saturnian system and did not find any candidates for such a resonance, although we cannot exclude this possibility due to uncertainties in the moons' shapes and gravity fields. There are however additional arguments against the relevance of such resonances to Enceladus. Most importantly, the dynamics of the Saturnian satellite system is dominated by Saturn's oblateness, and the mutual perturbations by the moons are less important than in the Uranian system. Therefore there is little or no coupling between the precession frequencies of different moons, and the relatively low mass of Enceladus is particularly unlikely to have a noticeable effect on the rotational dynamics of the larger moons. Furthermore, Ćuk et al. (2020) found that damping of inclination through a spin–orbit resonance can operate only until the relevant orbital inclination damps below the level at which resonance can be maintained, and this level is likely to be higher than the very low current i_E = 0008. Therefore we conclude that spin–orbit resonances are very unlikely to have damped the inclination of Enceladus.

Mimas. Our model requires that the eccentricity of Mimas predates the establishment of the current Mimas–Tethys and Enceladus–Dione resonances. The origin of this eccentricity may lie in the past 3:1 Mimas–Dione MMR (see Meyer & Wisdom 2008; Ćuk & El Moutamid 2022), or some other still unidentified, possibly three-body, resonance. The implication is that the damping of Mimas's inclination was limited and therefore that Mimas is unlikely to posses an internal ocean (Rhoden & Walker 2022). A new finding in our paper is that Mimas must have had an ≈01 inclination before encountering the 2:1 MMR with Tethys. The origin of this inclination is harder to explain, as most three-body resonances do not affect inclination, and a 3:1 MMR crossing with Dione produces a smaller inclination "kick" (Ćuk & El Moutamid 2022). One possibility that needs investigation is whether Dione was more inclined at the time of the Mimas–Dione 3:1 MMR crossing, and that would require producing a self-consistent model of the system's evolution further than 20 Myr into the past.

Enceladus. We find that the eccentricity of Enceladus was already excited (e ≤ 0.005) before encountering the 2:1 MMR with Dione, but its inclination was not. The most likely source of this eccentricity is three-body resonances, either isolated or as part of the 5:3 Dione–Rhea MMR crossing. The apparently pristine inclination of Enceladus implies that it did not encounter any major (first or second order) two-body resonances with mid-sized moons prior to the 2:1 MMR with Dione.

Tethys. Our model requires Tethys to have had a somewhat higher eccentricity (e_Θ = 0.002−0.003) 20 Myr ago, and almost its present large inclination. The eccentricity of Tethys poses no challenges, as a past higher eccentricity was suggested on both geophysical (Chen & Nimmo 2008) and dynamical (Ćuk et al. 2016) grounds, and subsequently tidal dissipation is likely to produce Tethys's present low eccentricity. The 1° inclination of Tethys was always known to have to predate the Mimas–Tethys 4:2 resonance (Allan 1969; Sinclair 1972), but its origin was never explained until Ćuk et al. (2016) proposed a Tethys–Dione secular resonance closely following (and dynamically related to) the Dione–Rhea 5:3 MMR crossing. Unless another mechanism of exciting Tethys's inclination is found, our initial conditions effectively require that a passage through the Dione–Rhea 5:3 had happened in the past.

Dione. We start our simulations with Dione that has an excited eccentricity (e_D ≈ 0.005) but very low inclination. This is generally consistent with a past secular resonance with Tethys (Ćuk et al. 2016), in which Dione "passed" its inclination and part of its eccentricity (acquired in a resonance with Rhea) to Tethys. The secular resonance is broken when either the eccentricity or inclination of Dione reach very low values (e.g., i_D < 01 for inclination). Therefore, a substantial eccentricity is to be expected to survive if the inclination is very low, implying that the resonance was broken by depletion of Dione's inclination.

Rhea. At the start of our simulation e_R = 0.002, 10 times in excess of Rhea's current free eccentricity e_R = 2 × 10⁻⁴. ⁵ Rhea's inclination at the start of our simulation is the same as now (i_R = 033) and does not change over the course of it. This substantial inclination is exactly what is expected from a past crossing of the Dione–Rhea 5:3 MMR. This dynamical mechanism is also expected to produce a comparable eccentricity of Rhea, broadly consistent with our initial conditions, but not the observed values. The discrepancy between theoretical expectations and the actual value here is significant, and requires a so-far unknown dynamical mechanism of lowering Rhea's eccentricity. We plan to address this issue in the near future using an integrator that resolves the rotational dynamics of Rhea and search for any additional dynamical features.

Our work is decidedly semiempirical in design, as we acknowledge the importance of the highly variable response of Saturn to the tidal forces of different moons. As we have shown in this paper, both equilibrium tides and resonance-lock-only models that were used so far fail to explain fully the system's dynamics. Apart from the general lack of information on the source of Saturn's dissipation, we tried to keep the number of free parameters to a minimum, leading to our decision to change the evolution rates abruptly. Our decision not to model the interior evolution of Enceladus but to adjust its tidal parameters periodically "by hand" was dictated by our own technical limitations, and we hope that in the future there will be integrated models that fully model both the orbital dynamics and interior evolutions of the moons.

Regardless of the nature and evolution of the resonant modes, it is clear that these are resonant phenomena, and therefore must exhibit a spike-like profile against frequency (possibly Lorentzian or similar). As apparent from Figure 10, we assumed that the width of the normal mode affecting Enceladus to be about 1% of its semimajor axis. We could in principle restrict the very fast evolution to a narrower interval that includes the initial encounter with the Enceladus–Dione 2:1 MMR, but that would require that the resonant mode and the resonance with Dione were encountered at the exact same time due to an unlikely coincidence. Furthermore, if Rhea is only passing though a resonant mode, this mode cannot be too narrow if we were to observe this transient phenomenon. These features in frequency space appear much wider than the resonant modes proposed by Fuller et al. (2016) based on the theory of tidal response in stars and giant planets. We also find that the tidal response at Tethys's frequency must be lower than what appears to be the "background" rate, so this frequency dependence is not restricted to peaks in dissipation, but also has local minima ("troughs"). Apart from being relatively wide, resonant modes cannot be too few and far in between if they were recently encountered by both Enceladus and Rhea as we propose.

Zoom In Zoom Out Reset image size

The determination of a consensus result for the current rate of tidal evolution of Titan is of great importance, as it will give us a major indication if resonance locking is present or not. Of all the moons of Saturn that produce significant tides (Mimas–Titan), Titan is most likely to be in a resonant lock, as it has the slowest (equilibrium) tidal evolution and is the least likely to have been recently reaccreted in some kind of late cataclysm. If Titan is not in a resonance lock, it is possible that the resonant modes are moving inward (in terms of semimajor axis at which they are encountered), making it somewhat less surprising that encounters between moons and modes are common. Still, two moons encountering resonant modes within the last few tens of millions of years appears unlikely unless the modes are very numerous, or their distributions are correlated with those of the moons. The last possibility may indicate that the moons may have reaccreted close to the modes, either due to where the last generation of satellites was before the assumed instability, or for some other reason. The only certainty is that the tidal evolution of the Saturnian system holds yet more surprises for us.

In this work we tried to consider all of the available constraints on the recent (last few tens of millions of years) and current orbital evolution of the major Saturnian satellites. We find that no single mechanism of orbital evolution proposed so far (including frequency-independent equilibrium tides and the evolution through resonance locking) can explain the orbits of these moons. Strong equilibrium tides can explain the existence of the observed resonances (Mimas–Tethys and Enceladus–Dione) and the current heating rate of Enceladus, while resonant modes are necessary to explain the current dynamics of Rhea and the original capture of Enceladus into the resonance with Dione. Additionally, the evolution of Tethys needs to be slower than that of other moons, implying "troughs" as well as "peaks" in the response as a function of frequency.

In order to reproduce the encounter of Enceladus with the 2:1 resonance with Dione successfully we require a passage through a resonant mode, rather than locking to a resonant mode. This is reasonable if the resonant modes in the inner system are evolving more slowly than those at larger distances, as originally predicted Fuller et al. (2016), but would not work in the context of inertial waves (Lainey et al. 2020). Alternatively, if Titan is currently not locked to a resonant mode as the results of Jacobson (2022) suggest, it is also possible that the resonant modes move inward, and in that case the moons could only temporarily cross the modes, rather than become locked to them.

Given the complexity and uncertainties of the tidal evolution rates that not only vary from moon to moon but also over time, it is difficult to reach firm conclusions about the age of the system. However, given the amount of dynamical excitation that the inner moons (especially Mimas and Enceladus) may have experienced in the last 20 Myr, it is difficult to envision that this system of relatively "dynamically cold" satellites evolved this way for hundreds of millions of years, let alone multiple gigayears. We hope that more precise future determinations of the current orbital evolution rates of the Saturnian moons (based on astrometry or spacecraft data) will be able to confirm or falsify our model of their recent evolution.