Orbital evolution of Saturn's mid-sized moons and the tidal heating of Enceladus
aa r X i v : . [ a s t r o - ph . E P ] S e p Orbital evolution of Saturn’s mid-sized moons and the tidal heatingof Enceladus
Ayano Nakajima a, ∗ , Shigeru Ida b , Jun Kimura c , Ramon Brasser b a Department of Earth and Planetary Sciences, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku,Tokyo 152-8551, Japan b Earth-Life Science Institute, Tokyo Institute of Technology, 2-12-1 Ookayama, Meguro-ku, Tokyo 152-8550, Japan c Department of Earth and Space Science, Osaka University, 1-1 Machikaneyama-cho, Toyonaka City, Osaka560-0043, Japan
Abstract
The formation and orbital evolution of Saturn’s inner mid-sized moons – Rhea, Dione, Tethys,Enceladus, and Mimas – are still debated. The most puzzling aspects are 1) how the Tethys-Dione pair and the Mimas-Enceladus pair passed through their strong 3:2 mean-motion resonancesduring the tidal orbital evolution, and 2) the current strong heat flow from Enceladus, which isa few orders of magnitude higher than the tidal energy dissipation caused by the present orbitaleccentricity of Enceladus. Here we perform N-body simulations of the moons’ orbital evolutionfrom various initial conditions – assuming that the moons were formed from Saturn’s hypotheticalmassive ring – and investigate possible paths to solve the above difficulties. If the moons remainon nearly circular orbits and the influence of the rings is neglected, we find that the Tethys-Dionepair cannot avoid becoming trapped in the 2:1 and 3:2 mean-motion resonances as they recede fromSaturn, and that the Tethys-Enceladus pair cannot avoid collisions after the resonance trapping,in case Saturn’s quality factor is smaller than 15 000. These findings are inconsistent with thecurrent orbital configuration. However, taking into account both the eccentricity excitation andthe orbital expansion caused by the ring torque, we find that these resonance captures are avoided.With the relatively high eccentricity pumped up by the torque, Enceladus passes through all themean-motion resonances with Tethys, and the Dione-Tethys pair passes through their 2:1 resonanceand possibly the 3:2 resonance as well. After Enceladus resides beyond the 2:1 resonance with theouter ring edge, the eccentricity can be tidally damped. While this is a promising path of evolution, ∗ Corresponding author
Email address: [email protected] (Ayano Nakajima)
Preprint submitted to Journal of L A TEX Templates September 5, 2018 n most runs, Enceladus collides with Tethys by the excited eccentricity. There is a hint that aring mass decrease (possibly due to Mimas formation) could avoid the collision between Enceladusand Tethys. The parameter survey taking into account detailed ring evolution and Mimas is leftfor future study. The heat that was tidally dissipated due to the eccentricity excitation by the ringtorque in the past is stored in the moons and slowly radiated away through conductive transfer.The stored heat in Enceladus may account for the current anomalously high heat flow.
Keywords:
Saturn, moons, resonance, Enceladus
1. Introduction and tidal evolution
The evolution and origin of Saturn’s mid-sized moons – Mimas, Enceladus, Tethys, Dione, andRhea – remain an enigma. Located closer than Saturn’s massive moon Titan, but farther awaythan Saturn’s famous ring system and a collection of much smaller moons, the classical mid-sizedmoons form a rich dynamical system both now and in the past.The masses and orbital elements of these moons are listed in Table 1. When ignoring mutualgravitational interactions and orbital eccentricities, the relative tidal expansion rate of each moon’ssemi-major axis ( a ) is given by (e.g., Murray & Dermott, 1999; Goldreich & Soter, 1966)1 a dadt = 3 k Q p M m M p (cid:18) R p a (cid:19) Ω , (1)where Ω = p G ( M p + M m ) /a is the orbital frequency of each moon, and G , M and R are thegravitational constant, mass and physical radius, respectively. Q p and k are the quality factorand Love number of the host planet (Saturn). Hereinafter, the subscript “m” and “p” indicatemoon and planet.Saturn’s k is often assumed to be in the range of 0.3-0.4 (Gavrilov & Zharkov, 1977; Helled & Guillot,2013). On the other hand, the quality factor of Saturn is not so well determined. Previously, itwas often thought that Q p ∼
18 000 is a lower bound inferred by assuming that Mimas migratedfrom the outside of the synchronous radius to its current orbit over the past 4.5 Gyr (Peale et al.,1980). Recently, however, Lainey et al. (2012, 2017) proposed a much lower, and controversial,value Q p = 1 700 ± Q p implies a much fastertidal evolution. Figure 1 shows a backward integration of Eq. (1) for each moon starting from theircurrent orbits down to the F-ring. Hereafter, we scale the semi-major axis of the moons by theaverage distance of the F-ring to Saturn ( a F ≃ . × , km ≃ . R p where R p is Saturn’s physicalradius), which is comparable to the planet’s Roche limit ( r Roche = 2 . R p ( ρ p /ρ m ) / , where ρ p and ρ m are the bulk density of Saturn and material, respectively).We used both Q p = 18 000 and Q p = 1 700, assuming k = 0 .
34 (Gavrilov & Zharkov, 1977)for all the moons. We use the same value of k /Q p for all the mid-sized moons in our simulations.For Rhea, Lainey et al. (2017) suggested that k /Q p is ten times smaller than we assumed. In oursimulation, we mostly focus on Enceladus, Tethys and Dione. Our assumption that all the moonshave the same k /Q p does not change the result.The lower value of Q p has profound implications for the formation of Saturn’s inner moons, notleast of which is that they can no longer be primordial. One theory for the formation of the moonswas proposed by Charnoz et al. (2011), who suggested that these mid-sized moons were formedrelatively recently from the spreading of a previous massive ring rather than from an extendedcircumplanetary protosatellite disk (also see Crida & Charnoz, 2012; Salmon & Canup, 2017). Thetheory suggested by Crida & Charnoz (2012) makes predictions about the mass and semi-majoraxis relationship of the moons that accurately fit their observed distribution, and Q p needs to besmall enough for Saturn’s mid-sized moons to be formed from the spreading of massive rings withinthe age of the Solar System. Recently Fuller et al. (2016) proposed that Q p gradually decreasedfrom the large values in an early phase of resonant locking between Saturn’s oscillation mode andthe moon’s orbital frequency in the course of Saturn’s interior evolution and tidal orbital evolution.In this case, the mid-sized moons were formed from an extended circumplanetary protosatellitedisk over 4.5 Gyr ago and the tidal orbital expansion was slow in the early phase until resonantlocking occurred. Further research is needed to distinguish between the two cases, though theyaren’t mutually exclusive.In the case that (1 /a )( da/dt ) of an outer moon is smaller than that of an inner moon, themigration is termed “convergent”, while it is “divergent” otherwise. Two or more moons can becaptured in a mean-motion resonance if the migration is convergent and “adiabatic”, i.e. when thecharacteristic libration timescale of the resonance angle is much shorter than the migration timescaleacross the resonant width (e.g., Murray & Dermott, 1999). Currently, Tethys is just inside of the3emi-major axis ( a F ) mass (10 − M p ) eccentricityRhea 3.77 4.07 0.0010Dione 2.70 1.94 0.0022Tethys 2.11 1.09 0.0000Enceladus 1.70 0.190 0.0045Mimas 1.33 0.067 0.0202 Table 1: Physical parameters of Saturn’s mid-sized moons. The semi-major axis is scaled by the orbital radius ofthe F ring ( a F ≃ . × km) and the masses are scaled by 10 − times Saturn’s mass M p . These values are citedfrom NASA Space Science Data Coordinated Archive. S e m i - m a j o r a x i s ( R F ) Time (Gyrs)
T-D 2:1T-D 3:2
Figure 1: Evolution of semi-major axis of Mimas, Enceladus, Tethys, Dione, and Rhea, backwardly integratingEq. (1) from the current semi-major axes without mutual interactions between the moons. The dashed lines assumed Q p = 18 000 and the solid lines assumed Q p = 1 700 which are normalized by the F-ring radius ( R F = a F = 1).We assumed Saturn’s Love number k = 0 .
34 for all the moons in both lines. The vertical black lines indicate thelocation of most recent 3:2 and 2:1 mean motion resonances between Tethys and Dione. Q p < ∼
15 000. Iftheir eccentricities remained low throughout their migration, it is not yet clear how they avoidedor escaped from these resonances. Even if the time-dependent Q p modeled by Fuller et al. (2016)is considered, the convergent migration near the 3:2 resonance does not change and the problem ofcapture in the 3:2 resonance cannot be easily avoided.Even when considering the formation of these moons, the resonance capture problem remains.Salmon & Canup (2017) performed N-body simulations of collisional growth of the mid-sized moons,based on Charnoz et al. (2011)’s model. While they have succeeded to reproduce the overall mass-distance distribution of the mid-sized moons, the capture of Tethys and Dione in the 3:2 resonancestill looks to be a setback in their results. Zhang & Nimmo (2012) pointed out the possibility that alarge ( ≈ ≈ . e -Dione resonance. If we consider such animpact which may occur during the satellite formation, the resonance can be broken and reproducethe current orbit, although the occurrence probability of such an impact is not clear.For relatively small Q p , Fig. 1 suggests that Enceladus was formed earlier than Tethys andTethys overtook Enceladus during their tidal migration. Note that Fig. 1 does not include the effectof mutual interactions between the moons. In section 3.1, we will show that such an orbital crossingcannot occur because of the resonant interaction between Enceladus and Tethys. We shall furtherpoint out that this orbital crossing is not required if we include the additional orbital expansion byring torque, which is very rapid near the ring’s outer edge. On the other hand, its force vanishesbeyond the 2:1 resonance with the edge of the ring (e.g., Charnoz et al., 2011; Crida & Charnoz,2012). We will also point out that eccentricity excitation by the ring torque (e.g., Goldreich & Sari,2003; Duffell & Chiang, 2015) may need to be taken into account in addition to its effect on theorbital expansion, although the eccentricity excitation was not considered in previous simulations(Charnoz et al., 2011; Crida & Charnoz, 2012; Salmon & Canup, 2017). If the eccentricity is excitedover a threshold value, the capture probability at a resonance is significantly reduced (e.g., Malhotra,1993) and a pair of moons can avoid becoming trapped.The high surface heat flux of Enceladus, which is observed to be ∼ . ± . ∼ . H = 212 Ω GM a k Q m (cid:18) R m a (cid:19) e ∼ . (cid:18) k /Q m − (cid:19) (cid:18) R m R E (cid:19) (cid:18) aa E (cid:19) / (cid:18) ee E (cid:19) GW , (2)where a E , e E (see Table 1) and R E ( ≃
250 km) are semi-major axis, eccentricity ( e ) and physicalradius ( R ) of Enceladus (the subscript “E” represents Enceladus) and we used Saturn’s mass for M p . Even with relatively high tidal damping parameters in Enceladus of k /Q m ∼ − , theestimated tidal heating is only ∼ . ≃ . H ≃ . /Q p ) GW for k /Q m ∼ − , assuming that thesystem is in a steady state. If the moons are in a steady state, e in Eq. (2) is proportional to k /Q p (see Appendix A). Since tidal heat production is inversely proportional to Q p , the heatingis as high as ∼
10 GW for Q p ∼ e than the current value of Enceladus, as shown in Eq. (2), unless the moon damping parametersare extremely high ( k /Q m ∼ − ). Although the extremely dissipative case of Enceladus is notcompletely ruled out (e.g. Ferraz-Mello et al., 2017; Choblet et al., 2017), a very strong dampingwould prevent the moons from passing mean-motion resonances as discussed in section 3.2.3.This inconsistency strongly suggests that the current heat flux from the surface is not equili-brated with the current heat production in Enceladus. Ojakangas & Stevenson (1986) proposedthe idea that Enceladus’ eccentricity increases most of the time with episodical decreases, generat-ing a large amount of heat; the current state is the end of the high heat generation phase with afully damped eccentricity. However, Meyer & Wisdom (2008) argued that such oscillation does notoccur with the Ojakangas-Stevenson model. O’Neill & Nimmo (2010) proposed an episodic heatrelease. On occasion, the generated heat is stored in the interior, and subsequently the stored heatis episodically released. This idea suggests that we are observing the narrow window of a high heatenergy release from Enceladus.Here we consider the possibility that the intense heat production during the past orbital evolu- This formula applies for a moon in a synchronous rotation. If libration is taken into account, tidal dissipationbecomes stronger (Ferraz-Mello et al., 2017, also see the comment below). e and the ring torque easily excites e to values of ∼ Q p ∼ e orbitalphase of Enceladus could be recent enough for Enceladus to keep the stored heat in its interior.Note that a close encounter between two moons can only pump the eccentricity of the smaller bodyup to a value that depends on the surface escape velocity ( v esc ) of the larger one. The correspondingeccentricity is ∼ v esc /v K = [2( M m /M p )( a m /R m )] / ∼ M m /M p ) / ( a m /a F )] / ∼ .
02 where v K is Keplerian velocity and a F is comparable to the Roche radius, a F ∼ r Roche ∼ . ρ p /ρ m ) / R p ∼ . (cid:0) ρ p R /ρR (cid:1) / R m ∼ . M p /M m ) / R m . The past high-eccentricity phase is not the result ofpast close encounters, but was instead caused by resonant interactions between moons or pumpingby the ring torque.Here, we investigate the orbital evolution of the mid-sized moons, based on the model of for-mation from the spreading ring. We employ N-body simulations to compute the evolution anddiscuss resonance capture and tidal heat production due to eccentricity evolution. We take intoaccount the changes in eccentricities and semi-major axes of the moons according to tides in Sat-urn and the moons themselves. We also perform runs where we include orbit torques causedby the ring, including potential eccentricity excitation. The formation models of Charnoz et al.(2011) and Crida & Charnoz (2012) were semi-analytical and did not calculate the gravitationalinteractions between the moons, including resonant configurations and eccentricity evolution. TheN-body simulations by Salmon & Canup (2017) were the first to investigate the formation scenarioof Charnoz et al. (2011) and Crida & Charnoz (2012) in more detail. Although they reproduced theoverall mass distribution of these moons, they did not discuss the details of resonance capture/break-up or eccentricity evolution in individual systems, both of which are important to discuss the currentorbital architecture of the moons.We focus our investigation on resonance passing/capture/break-up and tidal heat productionwith N-body simulations, and therefore we do not include early collisional growth of moonlets fromthe spreading ring.The outline of our paper is as follows. Section 2 describes our numerical model, and Section7 presents our results of the orbital evolution of the moons and discusses the resonance trappingin detail. In Section 4 we estimate the heat energy of Enceladus stored in the course of orbitalevolution. Finally, we present our conclusions and a discussion in the last section.
2. Methods
We simulate the orbital evolution of the system – which mainly consists of Saturn, Enceladus,Tethys, and Dione – to investigate the detailed orbital evolution of strongly interacting moonsstarting from many different initial conditions; in some runs we also added Rhea. Dynamically,Rhea is almost decoupled from other moons. Mimas has the smallest mass and could not affectother moons’ motions significantly.Figure 1 suggests that the Enceladus-Tethys pair undergoes orbital crossing if Q p < ∼
15 000.According to the recently proposed smaller Q p (Lainey et al., 2012, 2017), we adopt constant Q p ∼ C = 10 − , to reduce computation time. This “speed-up factor” has been used inmany other works (e.g. Malhotra & Dermott, 1990; Showman et al., 1997; Meyer & Wisdom, 2008;Zhang & Nimmo, 2009). The ratio between tidal orbital expansion and eccentricity damping ratesis kept the same for different values of the speed-up factor, in order to maintain consistency (seediscussions below). Furthermore, the time in all of following figures represents real time that isobtained by simulation time multiplied by the speed-up factor ( C ). The orbital expansion andeccentricity excitation by the ring torque are implemented in the N-body simulations, and the8cceleration with the same speed-up factor for tides is applied. Some aspects in orbital changes bymutual gravitational interactions, including resonant interactions, cannot be accelerated becausethey are calculated by N-body simulation. We will discuss the effects of the speed-up factor on theprobability of resonant trapping.In our simulations, we find that collisions between moons occur after orbital eccentricities areexcited by resonant perturbations or the ring torque. Consequently, the collision velocity ( v col ) isusually larger than the surface escape velocities ( v esc ) of the moons and the collisions are usually“hit-and-run” collisions (Asphaug et al., 2006). We use the model by Genda et al. (2012) based onSPH simulation results. The critical collision velocity ( v cr ), so that a collision between body 1 and2 (their masses are M and M ) is hit-and-run, is given by v cr v esc = c Γ Θ c + c Γ + c Θ c + c , (3)where Γ = | M − M | / ( M + M ), Θ = 1 − sin θ ( θ is impact angle), and v esc is escape velocity. Thefitting parameters are c = 2 . c = − . c = 1 . c = 1 .
08 and c = 2 .
50 (Genda et al.,2012). If v col < v cr , we assume that the collision results in merging, while it is hit-and-run otherwise.In this study, we consider the most optimal case to reproduce the current orbits, preserving themoons against catastrophic disruption. For simplicity, in the case of hit-and-run collisions, wehave the moons pass through each other in a softened gravitational potential without collisionalenergy dissipation. We find that once moons start orbit crossing, they undergo repeated hit-and-runcollisions, where we assume that they eventually coalesce, regardless of whether or not we includethe energy dissipation during the hit-and-run collision. Tidal deformation of the host planet caused by the moons transfers angular momentum from theplanetary spin to the moon’s orbit. To express the orbital expansion rate in Eq. (1), the tangentialforce per unit mass is added to the equations of motion, which is given by f p ,ψ ≃ a d p GM p adt = 12 a dadt r GM p a = 32 k Q p GM m a (cid:18) R p a (cid:19) ≃ k Q p GM m r (cid:18) R p r (cid:19) , (4)where r means the orbital radius of the moons. In our simulation, we set the radial component ofthe planetary tidal force to be zero, because it does not affect the orbital expansion.9 moon’s deformation caused by the planet dissipates the moon’s orbital kinetic energy, whichresults in a decrease in the moon’s eccentricity ( e ) and semi-major axis ( a ). The eccentricitydamping timescale is given by (e.g., Murray & Dermott, 1999)1 τ e = − e dedt = 21 k Q m M p M m (cid:18) R m a (cid:19) Ω . (5)For the current orbital elements of Enceladus and k m /Q m ∼ − , τ e ∼ × yrs. Using theeccentricity damping timescale, we add the following tidal force per unit mass caused by the moon’sdeformation in a simple form to the equations of motion (Kominami & Ida, 2002), f m = − v − v K τ e , (6)where v is the moon’s velocity and v K is the local circular Keplerian velocity. Adding the semi-majoraxis damping associated with the eccentricity damping, the semi-major axis changes as1 a dadt = 3 k Q p M m M p (cid:18) R p a (cid:19) Ω − k Q m M p M m (cid:18) R m a (cid:19) e Ω . (7)The second term in Eq. (7) corresponds to − e /τ e in Eq. (7). Because the moon’s orbital angularmomentum is conserved during the eccentricity damping, 0 = (1 /L ) dL/t = (1 / a ) da/dt − ( e/ ( e − de/dt ∼ a/ (2 τ a ) − e /τ e . We do not include the second term in f p ,ψ in Eq. (4), because it isautomatically caused by f m (Eq. (6)). The second term is dominant when e > ∼ " k /Q p k /Q m (cid:18) M m M p (cid:19) (cid:18) R p R m (cid:19) / ∼ . (cid:18) k /Q p − (cid:19) / (cid:18) k /Q m − (cid:19) − / (cid:18) M m /M p − (cid:19) / . (8) In one set of runs (SET2), we include the torque from Saturn’s ring following Crida & Charnoz(2012). The timescale is given by1 τ a, ring = 1 a dadt = a F a d (∆ a/a F ) dt = 1627 π M ring M p M m M p (cid:18) ∆ aa F (cid:19) − a F a Ω , (9)where ∆ a = a − a F is the separation from F-ring, M ring is Saturn’s ring mass, which may be compa-rable to the moon’s mass (Crida & Charnoz, 2012) (see discussion below). Note that this formulais an approximate one for computational simplicity and the actual torque decreases in a discretemanner (e.g., Meyer-Vernet & Sicardy, 1987). Because there is no 1st-order Lindblad resonancebeyond 2:1 resonance with the outer edge of the ring, we need to introduce this discreteness at leastbeyond the 2:1 resonance. We set the ring torque to vanish at a > ∼ . a F .10oldreich & Sari (2003) and Duffell & Chiang (2015) argued that the non co-orbital Lindbladtorque excites eccentricity with the timescale as follows,1 τ e, ring = 1 e dedt ≃ π M ring M p M m M p (cid:18) ∆ aa F (cid:19) − Ω , (10)where we used M ring ≃ π Σ a . In our simulation, we include the ring torque effects on semi-majoraxis and eccentricity evolution in similar ways to Eqs. (4) and (6).While the inner/outer Lindblad torques excite e , the co-orbital Lindblad torque and corotationresonance torque damp e . For a uniform surface density distribution without a gap, the damping isstronger. In the case of a gas disk, gas surface density is not zero even at the gap center. Whether e is excited or not depends on how deep the gap is. On the other hand, ring particles are almostcompletely empty outside the ring, so that e of satellites should be excited. Goldreich & Sari (2003)and Duffell & Chiang (2015) considered a gas giant planet that opens up a gap in a protoplanetarydisk. Both inner and outer Lindblad torques excite e . Because only Lindblad torque from an innerdisk exists for the ring torque, we decrease the formula by Goldreich & Sari (2003) by a factor of 2.The e -excitation by the ring torque was not taken into account in previous studies of Saturn’ssystem (Crida & Charnoz, 2012; Salmon & Canup, 2017). As will be shown later, the newly incor-porated e -excitation plays a key role in bypassing strong mean-motion resonances.Note that Eqs. (9) and (10) cannot be valid for ∆ a →
0, because ring particles are stronglyscattered by the moon for ∆ a < ∼ √ r H , where r H is the Hill radius of the moon defined by r H = ( M m / M p ) / a ≃ . M m / − M p ) / a (Ida & Nakazawa, 1989). We should either startintegration from ∆ a > √ r H or start integration from ∆ a = 0 with adding a softening parameter ǫ ∼ √ r H to ∆ a in an approach in which collisions are approximated with a softened gravitationalpotential. Although ǫ ∼ √ r H is physically justified, a softening parameter is also often introducedto secure numerical stability and it corresponds to the numerical resolution. By this reason, wehere adopt ǫ = 10 r H . Because this size is still much smaller than the typical migration distance ofsatellites, this choice does not change the results.Equation (10) suggests that the initial eccentricity of a satellite is important for the subsequenteccentricity evolution. When satellites are formed from the ring, smaller clumps would interactwith one another before coagulation to the satellites and the formed satellites interact with densityfluctuations in the ring edge. The interactions are chaotic and an orbital eccentricity of a few × r H /a would be excited when the orbital separation between the clumps or that between the satellite and11he ring edge is smaller than a few × r H (e.g., Petit & Henon, 1986; Ida, 1990), although N-bodysimulation would be required to prove this argument. Since r H /a ≃ . M m / − M p ) / , we setthe initial values of e as e ∼ . − . M ring ) throughout each run for simplicity. The assumed ringmasses in individual runs are shown in Table 2. The assumption of the constant M ring does notaffect our results, because 1) the decrease in M ring would not be significant during the evolution wesimulated and 2) the a - e evolution path driven by the ring torque is independent of M ring . The ringmass would largely change when new satellites are formed, while it also changes gradually throughviscous diffusion. Here, we do not create a new satellite during the individual runs. In SET2 runs,we start our simulations from the timing of the birth of Enceladus; Dione and Tethys are alreadyformed. The ring mass should decrease at the time of the formation of Mimas, which we do notinclude in our simulations. The a - e evolution path driven by the ring torque is independent of thespeed-up factor C and M ring , because C is multiplied to both Eq. (9) and Eq. (10) in the simulationsand both equations are proportional to M ring . As we will show, the key process to reproduce thecurrent orbital configurations is to avoid capture by mean-motion resonances. The capture can beavoided by sufficiently large e and/or fast migration (non-adiabatic migration) of satellites. The e evolution path as a function of a is independent of M ring and C in the phase dominated by thering torque. The migration speed is proportional to M ring . However, we will show that satellitesare trapped at a resonance that they encounter at the first time even with C = 10 − , whichmeans that the capture is not changed even by order of magnitude larger values of M ring in realcases with C = 1. Thus, the assumption of the constant ring mass and the artificial accelerationwith C are not critical flaws for the purpose of this paper, while simultaneous evolution of orbitsand the ring is important and left for future work.The phase dominated by the ring torque is determined as follows. Comparing Eq. (9) with thefirst term of Eq. (7), we find that the ring torque is dominant for∆ aa F < (∆ a ) crit ,a a F ≡ " π M ring M p Q p k (cid:18) aR p (cid:19) / (11) ≃ . (cid:18) k /Q p − (cid:19) − / (cid:18) M ring /M p − (cid:19) / (cid:18) aa F (cid:19) / , (12)where a F is the F-ring radius ( a F ≃ . R p ). If we use k /Q p ∼ − and M ring /M p ∼ − ,the solution to Eq. (12) with a = a F + (∆ a ) crit ,a is (∆ a ) crit ,a ≃ . a F . Because ∆ a ≃ . a F
12t the 2:1 resonance with the ring edge, the ring torque is dominant in the orbital expansion untilit diminishes at the 2:1 resonance as founded by Crida & Charnoz (2012, section 6.3). Since C ismultiplied to both Eq. (7) and Eq. (9) in the simulations, (∆ a ) crit ,a is independent of the valueof C . It depends on M ring , but only weakly. The eccentricity excitation competes with the tidaleccentricity damping. It dominates over the tidal eccentricity damping, as long as τ e, ring < τ e , thatis, ∆ aa F < (∆ a ) crit ,e a F ≡ " × × π M ring M p (cid:18) M m M p (cid:19) Q m k (cid:18) aR m (cid:19) / (13) ≃ . (cid:18) k /Q m − (cid:19) − / (cid:18) M ring /M p − (cid:19) / (cid:18) M m /M p − (cid:19) / (cid:18) aa F (cid:19) / . (14)For k /Q m ∼ − and M ring /M p ∼ M m /M p ∼ − , the solution to Eq. (14) with a = a F +(∆ a ) crit ,e is (∆ a ) crit ,e ≃ . a F . Again, (∆ a ) crit ,e is independent of the value of C and it dependson M ring only weakly. We performed four sets of simulations (SET1A, SET1B, SET2A, SET2B). Figure 1 suggeststhat Enceladus formed earlier than Tethys and they underwent orbit crossing, if we neglect the ringtorque. Because backward integration is not available back to the state before the orbit crossing,we examine many different initial conditions with a T , < a E , where the subscripts “T” and “E”represent Tethys and Enceladus, and “0” represents the initial values of the simulation. We call thisset of runs “SET1A”. We also carry out “SET1B” where Enceladus and Tethys are formed almostsimultaneously in a horseshoe orbit, a T , ∼ a E , . As we will show in the next section, many of thesesimulations produce a collision between Enceladus and Tethys, and therefore cannot reproduce thecurrent orbital configurations of Enceladus and Tethys.As we already pointed out, the torque from the ring would affect orbital evolution. Figure 2shows backward tidal orbital evolution with Q p = 4 000, taking into account the torque from thering. Because the evolution includes formation of Dione, Tethys, Enceladus and Mimas, we changedthe ring mass for this particular plot, while we use a constant ring mass in N-body simulations. Theinitial mass of the ring is the sum of the masses of Dione, Tethys and 4 × Enceladus, and when onemoon is swept out from the ring, the ring loses some of its own mass. Gravitational interactionsbetween the moons are also neglected in this backward integration. For these parameters, it is13uggested that Enceladus formed later than Tethys and they can avoid the orbit crossing. Becausethe orbital evolution by the ring torque is very rapid near the ring, it is likely that the rapid orbitalevolution of Tethys had already ended when Enceladus began to form. Accordingly, the Enceladus-Tethys pair’s migration is convergent in its early phase until Enceladus migrates beyond the 2:1resonance with the ring edge. Tethys-Dione migration is always convergent. Because the resonantand secular perturbations among these moons –which are not taken into account in Figure 2– arecomplicated, we also need to test various initial conditions with a T , > a E , . We call this set of runs“SET2.” We will perform runs with only semi-major axis expansion by the ring torque (SET2A)and others with both semi-major axis and eccentricity increases by the ring torque (SET2B). Wewill show that only the runs with both semi-major axis and eccentricity increasing by the ringtorque in SET2B potentially reproduce the current orbital configurations of the mid-sized moons. S e m i - m a j o r a x i s ( R F ) Time (Gyrs)
Figure 2: The same as Figure 1 except taking into account the torque from the ring. Q p = 4000 is used. Theorbits of Mimas, Enceladus, Tethys, Dione, and Rhea are represented by violet, green, red, blue and orange lines.The initial mass of the ring is the sum of the masses of Dione, Tethys and 4 × Enceladus, respectively, and when onemoon is swept out from the ring, the ring loses some of its own mass.
Initial masses of the moons are the same as the current masses. In the simulations here, k , k , Q p and Q m are set to be constant with time for all the moons and we adopt k = 0 . k = 10 − , and Q m = 100 in all runs. We use Q p = 1 700 − C = 10 − , depending on runs. 14 . Results Figure 3 shows a typical result of orbital evolution of SET1A. In this case, the ring torque isnot taken into account ( M ring = 0). We adopt Q p = 1 700 and C = 10 . Because interactions ofEnceladus, Tethys and Dione are essential for the orbital evolution in SET1A, we omit Mimas. S e m i - m a j o r a x i s ( R F ) Time (Gyrs)
TethysEnceladusDioneRhea
Figure 3: Evolution of semi-major axis, pericenter and apocenter of Enceladus, Tethys, Dione and Rhea in SET1A.The speed-up factor C for semi-major axis expansion and that for the eccentricity damping are 10 . Enceladusgets trapped in external 2:3 mean-motion resonance with Tethys. After the eccentricity of Enceladus is excited,hit-and-run collisions repeat a few dozens of times. We start simulations when Tethys is formed at a T , ∼ a F = 1. Enceladus was already formedand has migrated to a E , ∼ . a F in agreement with Figure 1. Because Tethys is 5.7 times moremassive (Table 1), it catches up with Enceladus (Eq. 1) and Enceladus gets trapped in outer 2:3mean-motion resonance with Tethys’ orbit at t ≃ .
13 Gyrs. Although the migration is acceleratedby a factor of C = 10 , Enceladus gets trapped in the 1st order mean-motion resonance that it firstencountered. The trapping occurs if the convergent migration speed is low enough (“adiabatic”)(Murray & Dermott, 1999). Thus, the trapping in 2:3 resonance would be robust in the realisticcase with C = 1. 15s Enceladus migrates while being trapped in the resonance, the eccentricity of Enceladus ( e E )secularly increases. The analytical prediction for j + 1 : j mean-motion resonance trapping givenby Eq. (A.3) in the Appendix A based on Malhotra (1995) shows de dt ≃ j + 1 C (cid:18) a T da T dt − a E da E dt (cid:19) ≃ C . j + 1 k Q p M T M p (cid:18) R p a T (cid:19) Ω T , (15)where we used Eq. (1) and assumed M T ≫ M E . The tidal e -damping is given by Eq. (5) as de dt = 2 e τ e = − C E M p R M E a k Q m e . (16)Equilibrating this excitation and damping with a E /a T = ( j + 1 /j ) / , the asymptotic value of e E is estimated as e E ∼ " j + 1) (cid:18) j + 1 j (cid:19) / (cid:18) M E M p (cid:19) / M T M E Q m /k Q p /k / ∼ . (cid:20) ( j + 1 /j ) / /j (cid:21) / (cid:18) M E /M p . × − (cid:19) / (cid:18) M T /M E . (cid:19) / (cid:18) k /Q m − (cid:19) − / (cid:18) k /Q p − (cid:19) / . (17)Although this analytical estimate includes uncertainty for such high e , the numerical simulation inFig. 3 actually shows that the eccentricity of Enceladus secularly increases toward a high value, until e E and e T become ∼ .
18 and ∼ .
04, respectively, and orbital crossing starts between Enceladusand Tethys at ∼ . C cancels out in this analytical estimate,suggesting that the equilibrium eccentricity would be similar in a real system with C = 1. Afterthat, Tethys and Enceladus repeat a few tens of hit-and-run collisions, because the collision velocityis excited by the resonant secular perturbations and is significantly larger than v esc .We performed 50 runs in SET1 (A and B) and a similar evolution was found in all casesexcept one run in which Enceladus was scattered to the inside of Tethys’ orbit without collisions.Initial conditions of semi-major axis of Tethys, Dione and Rhea, speed-up factor, mass of the ring,inclination of Tethys i T and number of simulations are listed on Table 2. In most runs, we adoptedthe same speed-up factor for tides of the moons and the planet. In some runs, we adopted differentvalues for the moons ( C m ) and the planet ( C p ). The varies between the runs are only initial orbitalangle of each moon. Even in the inwardly scattered case, Tethys and Enceladus collide many timesafter the inward scattering of Enceladus and before they become isolated. We also performed runswith non-zero energy dissipation at hit-and-run collisions. In those cases, they still repeat collisionsand eventually they merge because the collision velocity becomes smaller as the collisions repeat.16n some runs in SET1A, we include the ring torque in the initial condition of a E , > a T , . Wefound that the results are similar. Therefore, we conclude that runs in SET1A inevitably end upmerging or disrupting Enceladus and the current orbital configuration of the mid-sized moons isnever reproduced.In SET1B, Enceladus and Tethys migrate together in a horseshoe orbit in the early phase. Weperformed 5 runs of this case and the initial conditions of semi-major axis of Tethys, Dione andRhea, speed-up factor C , mass of the ring and number of simulations are listed on Table 2. However,Enceladus and Tethys eventually start orbit crossing and there are repetitive hit-and-run collisionsas in SET1A runs. Therefore, SET1B does not reproduce the current orbital configuration either. In SET2, we start from the initial conditions with a T , > a E , . We performed 30 runs in thisset: 20 runs considering only semi-major axis expansion by the ring torque (SET2A: subsection3.2.1 and 3.2.2) and 10 runs considering both semi-major axis expansion and eccentricity excitationby the ring torque (SET2B: subsection 3.2.3).The current orbital separation between Enceladus and Tethys is slightly smaller than their 3:2resonance. As Fig. 2 shows, if the ring torque can transfer Enceladus to the orbit inside the 3:2resonance with Tethys, the current orbital configurations may be reproduced. Even if the ring torquesends Enceladus to the orbit beyond 4:3, 5:4, or a higher- j resonance with Tethys, Enceladus is nottrapped at these resonances, because the Tethys-Enceladus migration eventually becomes divergentat a E > .
59. Then Enceladus divergently passes the closer resonances to Tethys and approachesthe 3:2 resonance again. From Eq. (1), in the regions where the ring torque is not effective, theratio of tidal orbital expansion rate of Tethys to that of Enceladus is described as follows(1 /a T )( da T /dt )(1 /a E )( da E /dt ) = M T M E (cid:18) a T a E (cid:19) − . . (18)Because M T /M E ≃ .
7, this ratio is > a T /a E < .
3. Because the3:2 resonance corresponds to a T /a E ≃ .
3, the Tethys-Enceladus pair would not be trapped atresonances deeper (smaller a T /a E ) than 3:2 resonance. Therefore, one of the key points in SET2is whether the ring torque can send Enceladus to an orbit closer to Tethys’s orbit than the 3:2resonance with Tethys. 17 .2.1. SET2A: Results with orbital expansion by ring torque As shown in Eq. (12), the ring torque dominates over the planetary tidal torque for the orbitalexpansion near the ring. Because the rate of orbital expansion by the ring torque is proportionalto (∆ a ) − , with ∆ a the distance from the ring edge (Eq. 9), the expansion proceeds very rapidlynear the ring edge and the Enceladus-Tethys migration is always initially convergent in the settingof SET2 ( a T , > a E , ). From Eq. (9), in the ring torque dominated region,(1 /a T )( da T /dt )(1 /a E )( da E /dt ) = M T M E (cid:18) a T a E (cid:19) . (cid:18) ∆ a T ∆ a E (cid:19) − . (19)Because M T /M E ≃ .
7, the migration becomes divergent when the Enceladus-Tethys pair migratesoutward and ∆ a T / ∆ a E decreases to be < ∼ t ≃ . × − Gyrs . This is the 1st-order mean-motion resonance that Enceladus meets in the first place during its orbital evolution.In this run, C = 10 . The probability of resonance trapping is higher for slower convergence ofthe migration. In the real system with C = 1, where migration is much slower, Enceladus shouldalso get trapped in the 1st-order mean-motion resonance that Enceladus meets at the first place.Enceladus’ eccentricity increases after the resonant trapping, according to Eq. (15).At t ∼ × − Gyrs, ∆ a T / ∆ a E becomes ≃ t = 0.06, 0.20 and 0.60, theEnceladus-Tethys pair pass through the 4:3, 7:5 and 3:2 mean-motion resonances, respectively, andtheir eccentricities are excited. Because Tethys is 5.7 times more massive than Enceladus, theexcited eccentricity is much larger for Enceladus than for Tethys. The maximum eccentricity is ∼ .
08 for Enceladus and ∼ .
02 for Tethys in this run. Enceladus can store much more heat in itsinterior than Tethys (see section 4).Thus, SET2A conditions have the potential to reproduce the current orbital configurations ofthe mid-sized moons and to account for the high thermal activity of Enceladus. We performed 4runs of this simulation and the initial conditions of semi-major axis of Tethys, mass of the ring andnumber of simulations are listed on Table 2. However, runs adding Dione manifest a new problem,as shown below. 18 -5 -4 -3 -2 -1 s e m i - m a j o r a x i s ( R F ) Time (Gyrs)
Tethys Enceladus -5 -4 -3 -2 -1 E cc en t r i c i t y Time (Gyrs)
Enceladus Tethys
Figure 4: Evolution of semi-major axis of Enceladus (green) and Tethys (red). To focus on the initial part ( < . .2.2. SET2A: Interaction with Dione Because Dione is more massive than Tethys/Enceladus and Fig. 2 suggests that Dione may nothave undergone orbit crossing with other moons, we set Dione at the time of birth of Enceladusat the location predicted by the backward integration in Fig. 2. For this setting, we performed 16runs. The initial conditions of semi-major axis of Tethys and Dione, speed-up factor C , mass ofthe ring and number of simulations are listed on Table 2.Figure 5 shows the result of the run in which Dione is added to the initial conditions of the runin Fig. 4. The early orbital evolutions of Enceladus and Tethys are the same as the case in Fig. 4:Enceladus is captured in a 6:5 mean-motion resonance with Tethys and its eccentricity is secularlyincreased. However, in this case, Tethys rapidly gets trapped in a 2:1 resonance with Dione. TheTethys-Dione migration is always convergent because Dione is only 1.8 times more massive thanTethys, while Dione has sufficiently larger a than Tethys. Since Dione slows down Tethys’ migrationthrough the resonance, the migration of Tethys and Enceladus remains convergent. As a result,the multi-resonants state of 6:5 for Enceladus-Tethys and 2:1 for Tethys-Dione is established andthis configuration is stable until the end of the simulation. Due to the resonant migration, theeccentricities of the moons, and that of Tethys in particular, are secularly increased to values thatare extremely high in comparison to the values at present. Currently, orbital separation betweenTethys and Dione is smaller than 3:2 resonance, which is inconsistent with the trapping in 2:1resonance obtained by the simulation. The other 9 runs show similar results. Therefore, how tobreak up the Tethys-Dione’s mean-motion resonance is a critical issue. For the resonant capture, adiabatic convergent migration is required. In addition to that, for thecapture, the eccentricity of a moon that encounters a mean-motion resonance with a more massivemoon of mass M m must be smaller than a critical value ( M m /M p ∼ × − for Dione), given by(Malhotra, 1993) e crit ≃ . (cid:20) j ( j + 1) M m M p (cid:21) / ∼ . (cid:20) j/ ( j + 1) . (cid:21) / (cid:18) M m /M p − (cid:19) / , (20)For e > e crit , the capture probability for the j + 1 : j resonance abruptly decays.We simulate the orbital evolution of the moons by varying the initial eccentricity of Enceladus e E , between 0.01 and 0.03, while the other moons’ eccentricities are set to 0. We performed 1120 S e m i - m a j o r a x i s ( R F ) Time (Gyrs) E cc en t r i c i t y Time (Gyrs)
Figure 5: Evolution of semi-major axis of Enceladus (green), Tethys (red), and Dione (blue). Enceladus is swept outfrom the ring soon after Tethys. Enceladus is captured in a 6:5 mean-motion resonance with Tethys (same as thecase without Dione). Tethys is captured in a 2:1 resonance with Dione. The resonances are kept to the end. Q p , mass of the ring and number of runs are listed on Table 2.In Fig. 6, we include the eccentricity excitation due to the ring torque with M ring = 4 M E and e E , = 0 .
01 for Enceladus. In the beginning of Fig. 6 (on the left of this figure), we acceleratedthe tidal orbital evolution with C = 10 , and adopted k /Q m as 10 − . Equation (14) predictsthat the eccentricity increases until a reaches (∆ a ) crit ,e ≃ . e E is already excited up to ∼ .
04 at the timingof Tethys-Dione 2:1 resonance passing. Because e E is well excited beyond e crit , Enceladus passesthrough 3:2, 4:3 and closer to 1st order resonances with Tethys, although the Enceladus-Tethysmigration is convergent until Enceladus reaches a E ∼ .
59 and the ring torque decays.Secular perturbation from Enceladus with the relatively high eccentricity enhances Tethys’ ec-centricity up to e T ∼ . ∼ e crit . Tethys encounters a 2:1 resonance with Dione at t ≃ .
04 Gyrs.At the resonance passage, a relatively large amount of angular momentum is exchanged betweenTethys and Dione and also between Tethys and Enceladus. Owing to e T ∼ e crit , Tethys successfullyavoided getting trapped in a 2:1 resonance with Dione.In Fig. 6, the tidal orbital evolution is accelerated from C = 10 (the left panels) to C = 10 (the right panels) to follow the whole orbital evolution. After a E exceeds 1.59 at 0.6 Gyr, theEnceladus-Tethys migration becomes divergent and it becomes impossible for them to be trappedin 1st-order resonances between them. The calculation of the following evolution shows that thefinal a T and a D are consistent with their current semi-major axes ( a T = 2 .
11 and a D = 2 . e T and e D , respectively) are enhanced too muchby the resonance in the simulation. Indeed, the Tethys-Dione pair gets trapped in a 5:3 resonanceat t ≃ .
45 Gyrs and the resonance configuration remains stable, while the current Tethys-Dioneseparation is slightly smaller than the 3:2 resonance. Immediately after that, Enceladus gets trappedin a 2nd-order 9:7 mean-motion resonance with Tethys. These resonances are not consistent withthe current orbit.In this case, a mechanism is needed to kick the Tethys-Dione pair out of their resonance, suchas the impact that created the Odysseus crater on Tethys (Zhang & Nimmo, 2012). If we considerthe growth of the moons by merging as they migrate outward, it is, in principle, possible thatthe moons avoid the resonance capture. However, N-body simulations of collisional growth of themid-sized moons by Salmon & Canup (2017) showed that capturing at Tethys-Dione’s 2:1 and 3:222 s e m i - m a j o r a x i s ( R F ) Time (Gyrs) s e m i - m a j o r a x i s ( R F ) Time (Gyrs) E cc en t r i c i t y Time (Gyrs) E cc en t r i c i t y Time (Gyrs)
Figure 6: The evolution of semi-major axis (upper) and eccentricity (lower) of Enceladus (green), Tethys (red), andDione (blue) with e E , = 0 .
01. In the left panels, the speed-up factor of C = 10 is used. Tethys passes 2:1 resonancewith Dione at t ≃ .
04 Gyrs. The results in the right panels start from the end state of the evolution in the leftpanels and use C = 10 . Note that the oscillation patterns of the eccentricities are modified by the change in C .Tethys is captured in a 5:3 resonance with Dione at t ≃ .
45 Gyrs and immediately afterward, Enceladus is trappedinto a 9:7 resonance with Tethys (the right panels). S e m i - m a j o r a x i s ( R F ) E cc en t r i c i t y Time (Gyrs)
Figure 7: Same as Fig. 6, except e E , = 0 .
03. In the left panels, Tethys passes 2:1 resonance with Dione at t ≃ . C = 10 . In this case, Tethys and Enceladus collide with each other at ∼ .
18 Gyr. However, when the eccentricityof Enceladus becomes smaller, they never collide and it can reproduce the orbital configuration (the right panels).At the end of the right panels ( ∼ . e E , = 0 .
03 and other parameters were the same as those in the results ofFig. 6. With the three times higher e E , , the eccentricity of Enceladus increases up to e E ∼ . t ≃ .
04 Gyrs, as e T is enhanced well beyond e crit , Tethys passes through the Tethys-Dione 2:1 resonance and angular momentum is transportedbetween the moons (the left panels). The eccentricities of all the moons are excited and oscillatesubstantially from t ≃ .
04 Gyrs to the point of a E ≃ .
59, where the excitation of eccentricity andsemi-major axis by the ring torque for Enceladus decays.In the middle panels, the following orbital evolution is calculated with C = 10 . As e D is excitedto ∼ . − .
15, Enceladus eventually collides with Tethys at t = 0 .
18 Gyrs. We performed 11 runswith orbital expansion and eccentricity excitation by the ring torque. In 2 runs with e E , < ∼ . e E , > ∼ . e E is excited so much that Enceladus undergoes collisions or closescattering with Tethys. Therefore, it is difficult for Tethys to successfully pass the 3:2 resonancewith Dione and simultaneously avoid a collision with Enceladus, as long as only Enceladus, Tethysand Dione are integrated.We consider a hypothetical case in which the eccentricity of Enceladus is smaller than in themiddle panels. In the right panels in Fig. 7, we artificially decreased e E by a factor of 2 from the finalstate of the left panels and calculated the following orbital evolution. Although we do not specifythe cause of the decrease, the ring mass decrease due to the birth of Mimas could be responsiblefor it. In this case, Enceladus and Tethys never collide. But, the eccentricity of Tethys excited bysecular perturbations from Enceladus is still large enough for Tethys to avoid getting trapped in the3:2 resonance with Dione. After passing the 3:2 resonance, both e T and e D are damped. Enceladusis eventually trapped into the 2:1 resonance with Dione, which is the current resonance relation,because both e D and e E are damped below e crit . After the trapping, it is predicted that e E increasesto an equilibrium value, e E ∼ .
04 (Appendix A), which is 10 times larger than the current value.Hence, Enceladus must have recently become trapped in the 2:1 resonance with Dione.So far, we have neglected Mimas because it is the smallest mid-sized moon. However, if thering is still massive enough after Mimas’ formation, the torque from the ring can be transportedto Enceladus, Tethys and Dione through Mimas, because Mimas is currently located within a 2:125esonance with the ring edge and should have suffered the ring torque throughout its entire orbitalevolution. Figure 1 suggests that Mimas-Enceladus encounters their 3:2 resonances at a similartime as the trapping of the Tethys-Dione pair at their 3:2 resonance. It is very likely that the ringtorque pumps up Mimas’ eccentricity to a value larger than e crit and the Mimas-Enceladus pairavoids the trapping at their 3:2 resonance. However, interactions among the four moons with thering torque are complicated. The ring mass should also change at the formation of Mimas. Becausethese investigations require much more parameter surveys, we leave them to future study.
4. Heat Flux
As we have shown, the moons would have undergone a high eccentricity phase in the past duringorbital evolution. As we show below, the heat generated during the high eccentricity phase can bestored in the interior and the current high heat flux can reflect the stored heat (the current heatgeneration is not balanced with the surface heat flux). From the numerical simulations, here wecalculate the stored heat energy for each moon, E ∼ Z Hdt, (21)where H is given by Eq. (2).Although Enceladus would have a subsurface ocean at present, most parts of bulk Enceladuswould be in a solid phase. For simplicity, we assume that conduction is a major heat transfermechanism in the interior of the mid-sized moons. We estimate the conduction timescale veryroughly. The thermal conductivity and the specific heat capacity of solid ice are 2 W/m K and 2000J/kg, and those of rock are 3 W/m K and 900 J/kg, respectively. Assuming the densities of ice androck are 1000 kg m − and ∼ − , we can estimate the volume fraction, x , of rock by bulkdensity of moons; ρ (g cm − ) ∼ (1 − x ) + 3 x = 1 + 2 x . Using the obtained x , the mean thermalconductivity and heat capacity are λ ∼ [2(1 − x ) + 3 x ] W / m K and ρc ∼ { (1 + 2 x ) × [2000(1 − x ) + 900 x ] } Jm − / K, respectively. The thermal diffusion coefficient of Enceladus is then κ = λρc ∼ ζ × − m / s , (22)where ζ = (2 + x ) / [(1 + 2 x )(2 − . x )] ∼
1. The thermal conduction timescale for a moon with aphysical surface radius of R is τ cond ∼ R κ ∼ . ζ − (cid:18) R
250 km (cid:19) Gyrs . (23)26 able 2: Initial conditions SET a E , ( a F ) a T , ( a F ) a D , ( a F ) a R , ( a F ) C m C p i T (degree) M ring Q p runs1.55 1.0 2.5 3.7 10 M E M E M E M E M E M E M E M E M E M E M E M E M E M E M E M E M E M E M E M E . ×
12 3.1Dione – 0.19 0.44 0.14
Table 3: Calculated heat flux of Enceladus, Tethys, and Dione based on our simulations of SET2. These values areshown in GW units.
Because τ cond may be longer than the age of Enceladus in the low Q p model, we simply assume thatmost of the heat energy generated during the high eccentricity phase is still stored in the interior.Then, the total heat flux is given by L ∼ Eτ cond . (24)For e max being a typical e of the high e phase, H max being the heat generation rate for e max , and∆ t being the duration of the high e phase, L ∼ H max (∆ t/τ cond ). From Eq. (2), for example, if e max ∼ . , k m /Q m ∼ − and (∆ t/τ cond ) ∼ .
1, then L ∼
16 GW. In other words, e ∼ . (cid:18) k /Q m − (cid:19) − (cid:18) R m (cid:19) − (cid:18) a m . × km (cid:19) − / (cid:18) ∆ t/τ cond . (cid:19) − (25)is need to reproduce the current heat flux.We integrate H with t to obtain the energy E from the data of the numerical simulations andthe total heat flux L for individual simulations are listed in Table 3. In the case correspondingto Fig. 4, the tidal heat is consistent with the current observed value. However, in the results inFigs. 5 and 6, Tethys is captured by a resonance with Dione and subsequently has higher tidal heatgeneration than Enceladus. In the case of Fig. 7, tidal heat of Enceladus is not inconsistent withthe observed value which is suggested by Howett et al. (2011); Spencer & Nimmo (2013). Whilethe heat flux of other moons is smaller than Enceladus, the heat flux of Tethys predicted by oursimulation is not sufficiently small, which could be consistent with the geological features of Tethys(Giese et al., 2007; Chen & Nimmo, 2008). The case of Fig. 7 may be the most preferable case notonly in final orbital configurations but also in the heat flux.28 . Conclusion and discussion Through N-body simulations, we have numerically investigated the orbital evolution of Saturn’smid-sized moons (mainly Dione, Tethys and Enceladus), under the influence of Saturn’s tidal force,tidal dissipation in the moons, and the torque exerted by its ring. Our work was based on the modelof the mid-sized moons having formed relatively recently from the spreading out of a massive ring,a theory that was proposed by Charnoz et al. (2011). We have performed 80 runs in total withvarious initial conditions at two different settings.If the ring torque is sufficiently weak, Enceladus must be formed prior to Tethys and scatteredinward across Tethys’ orbit. In this set of runs (SET1), we found that Enceladus is always trappedin the outer 1st-order mean-motion resonance with Tethys due to the rapid migration of Tethyswhich is more massive than Enceladus. Such resonance trapping is inevitable because of the tidalmigration timescale of these moons. The eccentricity of Enceladus is secularly increased by theresonant migration. When the eccentricity becomes large enough, Enceladus undergoes collisionswith Tethys until the end of simulation. Therefore, it is impossible for Enceladus to cross Tethys’orbit and become isolated from Tethys, and the current orbital configuration of Enceladus andTethys is never reproduced.If the ring mass is comparable to the mass of a forming moon, which is a reasonable assumption,the torque from the ring is strong enough that Enceladus can be formed after Tethys, and Enceladusneed not cross Tethys’ orbit. We also performed many simulations with many different initialconditions in the setting in which Enceladus is initially located inside Tethys’ orbit (SET2). Becausethe ring torque is very strong near the outer edge of the ring, Enceladus-Tethys migration isconvergent as long as Enceladus’ orbital radius is smaller than the 2:1 resonance with the ringedge. As a result, Enceladus is always trapped in the 1st-order mean-motion resonance with Tethysthat Enceladus meets in the first place. After the trapping, Enceladus’ eccentricity is secularlyincreased to undergo repeated hit-and-run collisions as in SET1 runs. However, for some range ofinitial conditions, the migration turns into a divergent one before the hit-and-run collisions start.After that, Enceladus is never trapped at 1st-order resonances with Tethys and orbital eccentricitydecays via the tidal damping, which results in the final orbits of Enceladus and Tethys beingconsistent with the current ones.However, Tethys is trapped in the 2:1 or 3:2 mean-motion resonance with Dione in this setof simulations, because their migration is always convergent and adiabatic. It is very difficult to29eproduce their current orbital separation closer than the 3:2 resonance relation between Dione andTethys.Fuller et al. (2016) considered time-dependent Q p , in which Q p > ∼ Q p ∼ O (1000) after resonant locking between the orbital frequency andSaturn’s oscillation mode. In this case, moons can be formed in an extended circumplanetary diskand orbit crossing of the moons does not occur. Even in this model, the migration between Dioneand Tethys is usually convergent and adiabatic and it is difficult for them to avoid becoming trappedin the 3:2 resonance.We found that, if eccentricity excitation by the ring torque is effective, it has the potential tosolve the problem. This excitation is effective only in the regions close to the outer edge of thering. Owing to the excited eccentricity, Enceladus can easily pass through the resonances withTethys. Enceladus’ eccentricity could be mostly enhanced by the ring torque, but not by theresonant perturbations. The modest eccentricity of Tethys raised by the secular perturbation fromEnceladus with relatively high eccentricity breaks the 2:1 resonance between Tethys and Dione. Inpart of the range of Enceladus’ eccentricity, Tethys can also pass the 3:2 resonance with Dione. Asthe distance between Enceladus and the ring increases, tidal eccentricity damping dominates overthe excitation by ring torque, and the moons’ eccentricities decay, which enables the Enceladus-Dione pair to get trapped in the 2:1 resonance.This orbital evolution path is promising to reproduce the current orbits of Enceladus, Tethys,and Dione. However, in the calculations consisting only of Enceladus, Tethys and Dione with aconstant ring mass, we only found the orbital evolution path if Enceladus’ eccentricity is artificiallydecreased after the Tethys-Dione 2:1 resonance is passed; this is because Enceladus’ eccentricitynecessary to bypass the Tethys-Dione 3:2 resonance is so large that it results in collisions betweenEnceladus and Tethys afterward. Although we assumed that the ring mass is constant to focusourselves on orbital evolution, it must evolve with time. The ring mass decrease due to Mimasformation could lower Enceladus’ eccentricity to avoid the collision between Enceladus and Tethys.sThe dynamical effect of Mimas, the smallest moon among the mid-sized moons, could also play animportant role, although we have not explored the effect. The assumption of constant k /Q m forall the moons is also too simple. The parameter survey taking these effects into account is left fora separate paper.We also estimated the heat flux of individual moons, which is caused by thermal energy stored30n past periods of high eccentricity as a result of resonant interactions and ring torque. UnlessEnceladus and Tethys are captured by the mean-motion resonance, the heat generation is thehighest for Enceladus. The orbital evolution with the eccentricity excitation by the ring torque canproduce the heat flux that is comparable to or slightly smaller than the observationally inferredvalue of Enceladus. Tidal heating due to the high eccentricity events may make the moons moreconvective and dissipative, which may significantly increase k /Q m . We will also address theseissues in a separate paper.In conclusion, if we take into account of the orbital expansion and eccentricity excitation byring torque, there will be one possible pass to solve the problem of the current heat budged onEnceladus and the resonance capture from the birth to the current orbit. Acknowledgments
We thank anonymous referees for their helpful and very detailed comments. This work wassupported by JSPS Kakenhi grant 15H02065 and 17K05635. We thank Daigo Shoji for helpfulcomments.
Appendix A. Equilibrium eccentricity of a resonant pair
We consider a system of an outwardly migrating planet in a circular orbit and a test particletrapped at the exterior j : j + 1 mean-motion resonance with the planet. According to the outwardmigration of the planet, the trapped test particle also migrates outward and its eccentricity issecularly increased. The test particle’s eccentricity increase rate is given by (Malhotra, 1995) de dt ≃ j + 1 1 a dadt , (A.1)where a is the semi-major axis of the test particle. From this relation, it is suggested that for aconvergent resonant pair of bodies with mass M i , semi-major axis a i and eccentricity e i ( i = 1 , M /a / > M /a / ), their eccentricityincrease rates are de dt ≃ M M + M j + 1 (cid:18) a da dt − a da dt (cid:19) , (A.2) de dt ≃ M M + M j + 1 (cid:18) a da dt − a da dt (cid:19) , (A.3)31here 1 a da dt − a da dt = 3 k Q p R M p (cid:18) M a Ω − M a Ω (cid:19) > . (A.4)We neglected the second term in Eq. (7) for simplicity, because it is smaller than the first term for e ∼ .
03 that we consider here (Eq. (8)).The eccentricity damping rates by tide are de dt ≃ − k , Q m , M p M R m , a Ω e , (A.5) de dt ≃ − k , Q m , M p M R m , a Ω e . (A.6)By balancing (A.2) and (A.3) with (A.5) and (A.6), e has equilibrium values, e ∼ M M + M j + 1) k /Q p k , /Q m , (cid:18) M M p (cid:19) (cid:18) R p R m , (cid:19) " − M M (cid:18) a a (cid:19) / ∼ M M + M j + 1) k /Q p k , /Q m , R m , R p " − M M (cid:18) jj + 1 (cid:19) / , (A.7) e ∼ M M + M j + 1) k /Q p k , /Q m , R m , R p " M M (cid:18) j + 1 j (cid:19) / − . (A.8)For the Enceladus and Dione pair trapped at 2:1 resonance, j = 1 and body 1 and 2 are Enceladusand Dione, respectively. Substituting M /M = M E /M D ∼ . R E /R p ∼ . R D /R p ∼ . e E ∼ . (cid:18) k /Q p − (cid:19) / (cid:18) k /Q m − (cid:19) − / , (A.9) e D ∼ . (cid:18) k /Q p − (cid:19) / (cid:18) k /Q m − (cid:19) − / , (A.10)where we used k ∼ . Q p ∼ k ∼ − and Q m ∼ . These eccentricities are one orderhigher than the current values, one suggestion is raised that their eccentricities are now on the wayto the equilibrium.Note that the derivation for the equilibrium eccentricity here is simplified. More rigorous deriva-tions with Lagrange equations and detailed resonant properties are found in the past literatures(e.g., Meyer & Wisdom, 2008; Zhang & Nimmo, 2009). While the numerical factors differ fromthe rigorous treatment by a factor of up to a few, the dependence on k , p /Q p and k , m /Q m arereproduced by the simple derivation here and the difference in the numerical factors does not affectthe discussions in this paper. 32 eferencesReferences Asphaug, E., Agnor, C. B., & Williams, Q. (2006). Hit-and-run planetary collisions.
Nature , ,155–160. doi: .Charnoz, S., Crida, A., Castillo-Rogez, J. C., Lainey, V., Dones, L., Karatekin, ¨O., Tobie, G.,Mathis, S., Le Poncin-Lafitte, C., & Salmon, J. (2011). Accretion of Saturn’s mid-sized moonsduring the viscous spreading of young massive rings: Solving the paradox of silicate-poorrings versus silicate-rich moons. Icarus , , 535–550. doi: . arXiv:1109.3360 .Chen, E. M. A., & Nimmo, F. (2008). Implications from Ithaca Chasma for the thermal and orbitalhistory of Tethys. grl , , L19203. doi: .Choblet, G., Tobie, G., Sotin, C., Bˇehounkov´a, M., ˇCadek, O., Postberg, F., & Souˇcek, O. (2017).Powering prolonged hydrothermal activity inside Enceladus. Nature Astronomy , , 841–847.doi: .Crida, a., & Charnoz, S. (2012). Formation of regular satellites from ancient mas-sive rings in the solar system. Science (New York, N.Y.) , , 1196–9. URL: . doi: . arXiv:1301.3808 .Duffell, P. C., & Chiang, E. (2015). Eccentric Jupiters Via Disk-Planet Interactions. Astro-physical Journal , , 1DUMMY. URL: http://dx.doi.org/10.1088/0004-637X/812/2/94 .doi: . arXiv:1507.08667 .Duncan, M. J., Levison, H. F., & Lee, M. H. (1998). A Multiple Time Step Symplectic Algo-rithm for Integrating Close Encounters. The Astronomical Journal , , 2067–2077. URL: http://adsabs.harvard.edu/cgi-bin/nph-data{_}query?bibcode=1998AJ....116.2067D{&}link{_}type=ABSTRACT{%}5Cnpapers3://publication/doi/10.1086/300541 .doi: .Ferraz-Mello, S., Folonier, H. A., & Andrade-Ines, E. (2017). Tidal synchronization of close-insatellites and exoplanets. III. Tidal dissipation revisited and application to Enceladus. ArXive-prints , . arXiv:1707.09229 . 33uller, J., Luan, J., & Quataert, E. (2016). Resonance locking as the source of rapid tidal migrationin the Jupiter and Saturn moon systems.
Monthly Notices of the Royal Astronomical Society , , 3867–3879. doi: . arXiv:1601.05804 .Gavrilov, S. V., & Zharkov, V. N. (1977). Love numbers of the giant planets. Icarus , , 443–449.doi: .Genda, H., Kokubo, E., & Ida, S. (2012). Merging Criteria for Gi-ant Impacts of Protoplanets. The Astrophysical Journal , , 137. URL: http://stacks.iop.org/0004-637X/744/i=2/a=137?key=crossref.13660ad625848d1eee511669b9db133a .doi: . arXiv:arXiv:1109.4330v1 .Giese, B., Wagner, R., Neukum, G., Helfenstein, P., & Thomas, P. C. (2007). Tethys: Lithosphericthickness and heat flux from flexurally supported topography at Ithaca Chasma. grl , , L21203.doi: .Goldreich, P., & Sari, R. (2003). Eccentricity Evolution for Planets in Gaseous Disks. The Astro-physical Journal , , 1024–1037. URL: http://stacks.iop.org/0004-637X/585/i=2/a=1024 .doi: . arXiv:0202462 .Goldreich, P., & Soter, S. (1966). Q in the solar system. Icarus , , 375–389. URL: http://linkinghub.elsevier.com/retrieve/pii/0019103566900510 .doi: .Helled, R., & Guillot, T. (2013). Interior models of saturn: Including the uncertaintiesin shape and rotation. Astrophysical Journal , . doi: . arXiv:arXiv:1302.6690v1 .Howett, C. J. A., Spencer, J. R., Pearl, J., & Segura, M. (2011). High heat flow from Enceladus’south polar region measured using 10-600 cm-1 Cassini/CIRS data. Journal of GeophysicalResearch E: Planets , , 1–15. doi: .Ida, S. (1990). Stirring and dynamical friction rates of planetesimals in the solar gravitational field. Icarus , , 129–145. doi: .Ida, S., & Nakazawa, K. (1989). Collisional probability of planetesimals revolving in the solargravitational field. III. Astronomy and Astrophysics (ISSN 0004-6361) , , 303–315. URL:34 ttp://adsabs.harvard.edu/abs/1989A{&}A...224..303I{%}5Cnpapers3://publication/uuid/6C8ECD9B-6E4A-4F17-8D88-CF14E571F7C3 .doi: . arXiv:arXiv:1011.1669v3 .Kominami, J., & Ida, S. (2002). The Effect of Tidal Interaction with aGas Disk on Formation of Terrestrial Planets. Icarus , , 43–56. URL: http://linkinghub.elsevier.com/retrieve/doi/10.1006/icar.2001.6811 .doi: .Lainey, V., Jacobson, R. A., Tajeddine, R., Cooper, N. J., Murray, C., Robert, V., Tobie, G., Guil-lot, T., Mathis, S., Remus, F., Desmars, J., Arlot, J. E., De Cuyper, J. P., Dehant, V., Pascu,D., Thuillot, W., Poncin-Lafitte, C. L., & Zahn, J. P. (2017). New constraints on Saturn’s inte-rior from Cassini astrometric data. Icarus , , 286–296. doi: . arXiv:1510.05870 .Lainey, V., Karatekin, ¨O., Desmars, J., Charnoz, S., Arlot, J.-E., Emelyanov,N., Le Poncin-Lafitte, C., Mathis, S., Remus, F., Tobie, G., & Zahn, J.-P. (2012). Strong Tidal Dissipation in Saturn and Constraints on Enceladus’Thermal State From Astrometry. The Astrophysical Journal , , 14. URL: http://stacks.iop.org/0004-637X/752/i=1/a=14?key=crossref.d68a903d62213eb4f281cc2f4669349f .doi: . arXiv:1204.0895 .Malhotra, R. (1993). Orbital resonances in the solar nebula - Strengths and weaknesses. Icarus , , 264. doi: .Malhotra, R. (1995). The Origin of Pluto’s Orbit: Implications for the Solar System BeyondNeptune. Astronomical Journal , , 420. doi: . arXiv:9504036 .Malhotra, R., & Dermott, S. F. (1990). The role of secondary resonances in the orbital history ofMiranda. Icarus , , 444–480. doi: .Meyer, J., & Wisdom, J. (2007). Tidal heating in Enceladus. Icarus , , 535–539.doi: .Meyer, J., & Wisdom, J. (2008). Episodic volcanism on Enceladus: Application of the OjakangasStevenson model. Icarus , , 178–180. doi: .35eyer, J., & Wisdom, J. (2008). Tidal evolution of Mimas, Enceladus, and Dione. Icarus , ,213–223. doi: .Meyer-Vernet, N., & Sicardy, B. (1987). On the physics of resonant disk-satellite interaction. Icarus , , 157–175. doi: .Murray, C. D., & Dermott, S. F. (1999). Solar system dynamics .Ojakangas, G. W., & Stevenson, D. J. (1986). Episodic volcanism of tidally heated satellites withapplication to Io.
Icarus , , 341–358. doi: .O’Neill, C., & Nimmo, F. (2010). The role of episodic overturn in generating thesurface geology and heat flow on Enceladus. Nature Geoscience , , 88–91. URL: http://dx.doi.org/10.1038/ngeo731 . doi: .Peale, S. J., Cassen, P., & Reynolds, R. (1980). Tidal dissipation, orbital evo-lution, and the nature of Saturn’s inner satellites. Icarus , , 1196–1205.URL: .doi: .Petit, J.-M., & Henon, M. (1986). Satellite encounters. Icarus , , 536–555.doi: .Porco, C. C., Helfenstein, P., Thomas, P. C., Ingersoll, A. P., Wisdom, J., West, R., Neukum,G., Denk, T., Wagner, R., Roatsch, T., Kieffer, S., Turtle, E., McEwen, A., Johnson, T. V.,Rathbun, J., Veverka, J., Wilson, D., Perry, J., Spitale, J., Brahic, A., Burns, J. A., DelGenio,A. D., Dones, L., Murray, C. D., & Squyres, S. (2006). Cassini observes the active south pole ofenceladus. Science , , 1393–1401. doi: .Salmon, J., & Canup, R. M. (2017). Accretion of Saturn’s Inner Mid-sized Moons from a Massive Primordial Ice Ring, . (pp. 1–28). URL: http://arxiv.org/abs/1702.04385{%}0Ahttp://dx.doi.org/10.3847/1538-4357/836/1/109 .doi: . arXiv:1702.04385 .Showman, A. P., Stevenson, D. J., & Malhotra, R. (1997). CoupledOrbital and Thermal Evolution of Ganymede. Icarus , , 367–383.36RL: http://linkinghub.elsevier.com/retrieve/pii/S001910359795778X .doi: .Spencer, J. R., & Nimmo, F. (2013). Enceladus: An Active Ice World inthe Saturn System. Annual Review of Earth and Planetary Sciences , , 693–717. URL: .doi: .Zhang, K., & Nimmo, F. (2009). Recent orbital evolution and the internal structures of Enceladusand Dione. Icarus , , 597–609. URL: http://dx.doi.org/10.1016/j.icarus.2009.07.007 .doi: .Zhang, K., & Nimmo, F. (2012). Late-stage impacts and the orbital and thermal evolution ofTethys. Icarus , , 348–355. URL: http://dx.doi.org/10.1016/j.icarus.2011.12.013 .doi:10.1016/j.icarus.2011.12.013