{'fraction_non_alphanumeric': 0.037374577500929054, 'fraction_numerical': 0.01985524429736856, 'mean_word_length': 4.383442888444318, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Orbital evolution of Saturn's mid-sized moons and the tidal heating of Enceladus

The formation and orbital evolution of Saturn's inner mid-sized moons 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 resonances during the tidal orbital evolution, and 2) the current strong heat flow from Enceladus, which is a few orders of magnitude higher than the tidal energy dissipation caused by the present orbital eccentricity of Enceladus. Here we perform N-body simulations of the moons' orbital evolution from various initial conditions -- assuming that the moons were formed from Saturn's hypothetical massive ring -- and investigate possible paths to solve the above difficulties. If the influence of the rings is neglected, we find that the Tethys-Dione pair cannot avoid becoming trapped in the mean-motion resonances as they recede from Saturn, 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 the current orbital configuration. However, taking into account both the eccentricity excitation and the orbital expansion caused by the ring torque, we find that these resonance captures are avoided. With the high eccentricity pumped up by the torque, Enceladus passes through the resonances with Tethys, and the Dione-Tethys pair passes through their 2:1 and possibly the 3:2 resonance as well. After Enceladus resides beyond the 2:1 resonance with the outer ring edge, the eccentricity can be tidally damped. While this is a promising path of evolution, in most runs, Enceladus collides with Tethys by the excited eccentricity. There is a hint that a ring mass decrease could avoid the collision. The tidal heat due to the eccentricity excitation by the ring torque is stored in the moons and slowly radiated away. It may account for the current high heat flow.

a few orders of magnitude higher than the tidal energy dissipation caused by the present orbital eccentricity of Enceladus. Here we perform N-body simulations of the moons' orbital evolution from various initial conditions -assuming that the moons were formed from Saturn's hypothetical massive ring -and investigate possible paths to solve the above difficulties. If the moons remain on nearly circular orbits and the influence of the rings is neglected, we find that the Tethys-Dione pair cannot avoid becoming trapped in the 2:1 and 3:2 mean-motion resonances as they recede from Saturn, 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 the current orbital configuration. However, taking into account both the eccentricity excitation and the 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 the mean-motion resonances with Tethys, and the Dione-Tethys pair passes through their 2:1 resonance and possibly the 3:2 resonance as well. After Enceladus resides beyond the 2:1 resonance with the outer ring edge, the eccentricity can be tidally damped. While this is a promising path of evolution,

Introduction and tidal evolution
The evolution and origin of Saturn's mid-sized moons -Mimas, Enceladus, Tethys, Dione, and Rhea -remain an enigma. Located closer than Saturn's massive moon Titan, but farther away than Saturn's famous ring system and a collection of much smaller moons, the classical mid-sized moons 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 mutual gravitational interactions and orbital eccentricities, the relative tidal expansion rate of each moon's semi-major axis (a) is given by (e.g., Murray & Dermott, 1999;Goldreich & Soter, 1966) 1 a where Ω = G(M p + M m )/a 3 is the orbital frequency of each moon, and G, M and R are the gravitational constant, mass and physical radius, respectively. Q p and k 2p are the quality factor and Love number of the host planet (Saturn). Hereinafter, the subscript "m" and "p" indicate moon and planet.
Saturn's k 2p 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, it was often thought that Q p ∼ 18 000 is a lower bound inferred by assuming that Mimas migrated from 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. (2012Lainey et al. ( , 2017 proposed a much lower, and controversial, value Q p = 1 700 ± 500, which was deduced from the detailed analysis of astrometric observational data of the tidal evolution of Saturn's major moons. This lower value of Q p implies a much faster tidal evolution. Figure 1 shows a backward integration of Eq. (1) for each moon starting from their current orbits down to the F-ring. Hereafter, we scale the semi-major axis of the moons by the average distance of the F-ring to Saturn (a F ≃ 1.40×10 5 , km ≃ 2.4R p where R p is Saturn's physical radius), which is comparable to the planet's Roche limit (r Roche = 2.4R p (ρ p /ρ m ) 1/3 , 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 2p = 0.34 (Gavrilov & Zharkov, 1977) for all the moons. We use the same value of k 2p /Q p for all the mid-sized moons in our simulations.
For Rhea, Lainey et al. (2017) suggested that k 2p /Q p is ten times smaller than we assumed. In our simulation, we mostly focus on Enceladus, Tethys and Dione. Our assumption that all the moons have the same k 2p /Q p does not change the result.
The lower value of Q p has profound implications for the formation of Saturn's inner moons, not least of which is that they can no longer be primordial. One theory for the formation of the moons was proposed by Charnoz et al. (2011), who suggested that these mid-sized moons were formed relatively recently from the spreading of a previous massive ring rather than from an extended circumplanetary protosatellite disk (also see Crida & Charnoz, 2012;Salmon & Canup, 2017). The theory suggested by Crida & Charnoz (2012) makes predictions about the mass and semi-major axis relationship of the moons that accurately fit their observed distribution, and Q p needs to be small enough for Saturn's mid-sized moons to be formed from the spreading of massive rings within the age of the Solar System. Recently Fuller et al. (2016) proposed that Q p gradually decreased from the large values in an early phase of resonant locking between Saturn's oscillation mode and the 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 protosatellite disk over 4.5 Gyr ago and the tidal orbital expansion was slow in the early phase until resonant locking occurred. Further research is needed to distinguish between the two cases, though they aren't mutually exclusive.
In the case that (1/a)(da/dt) of an outer moon is smaller than that of an inner moon, the migration is termed "convergent", while it is "divergent" otherwise. Two or more moons can be captured in a mean-motion resonance if the migration is convergent and "adiabatic", i.e. when the characteristic libration timescale of the resonance angle is much shorter than the migration timescale across the resonant width (e.g., Murray & Dermott, 1999). Currently, Tethys is just inside of the  (1) from the current semi-major axes without mutual interactions between the moons. The dashed lines assumed Qp = 18 000 and the solid lines assumed Qp = 1 700 which are normalized by the F-ring radius (R F = a F = 1).
We assumed Saturn's Love number k 2p = 0.34 for all the moons in both lines. The vertical black lines indicate the location of most recent 3:2 and 2:1 mean motion resonances between Tethys and Dione.
3:2 resonance with Dione and their migration is convergent. Figure 1 shows that the Tethys-Dione pair should have passed a strong 3:2 resonance and also their 2:1 resonance if Q p < ∼ 15 000. If their eccentricities remained low throughout their migration, it is not yet clear how they avoided or 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 of capture 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)  impact which may occur during the satellite formation, the resonance can be broken and reproduce the 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 and Tethys overtook Enceladus during their tidal migration. Note that Fig. 1 does not include the effect of mutual interactions between the moons. In section 3.1, we will show that such an orbital crossing cannot occur because of the resonant interaction between Enceladus and Tethys. We shall further point out that this orbital crossing is not required if we include the additional orbital expansion by ring torque, which is very rapid near the ring's outer edge. On the other hand, its force vanishes beyond 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 the orbital 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 excited over 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 ∼ 15.9 ± 3.1 GW in the South Polar Terrain by the Cassini CIRS instrument (Howett et al., 2011), is also a big mystery.
If Enceladus is in a thermal equilibrium, the heat flux must be equal to the heat production in its interior. However, it is unlikely that Enceladus is in thermal equilibrium. The present radiogenic heating within the rocky core is estimated to be only ∼ 0.3 GW (Porco et al., 2006). The tidal energy dissipation rate is (Murray & Dermott, 1999) where a E , e E (see Table 1) and R E (≃ 250 km) are semi-major axis, eccentricity (e) and physical radius (R) of Enceladus (the subscript "E" represents Enceladus) and we used Saturn's mass for This inconsistency strongly suggests that the current heat flux from the surface is not equilibrated with the current heat production in Enceladus. Ojakangas & Stevenson (1986) proposed the idea that Enceladus' eccentricity increases most of the time with episodical decreases, generating a large amount of heat; the current state is the end of the high heat generation phase with a fully damped eccentricity. However,  argued that such oscillation does not occur with the Ojakangas-Stevenson model. O'Neill & Nimmo (2010) proposed an episodic heat release. On occasion, the generated heat is stored in the interior, and subsequently the stored heat is episodically released. This idea suggests that we are observing the narrow window of a high heat energy release from Enceladus.
Here we consider the possibility that the intense heat production during the past orbital evolu-tion of Enceladus is the result of eccentricity excitation due to the ring torque, and that this heat is still stored in Enceladus' interior and is slowly being released. Enceladus' potential past high eccentricity caused by the suggested ring torque is damped after it migrates beyond the 2:1 resonance with the ring edge, as we will clearly show in Section 3.2.3. Since the tidal heat rate is proportional to e 2 and the ring torque easily excites e to values of ∼0.05-0.1, the stored heat can be large enough to account for the current high heat flux. With the recently proposed Q p ∼ 2 000, the high e orbital phase 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 body up to a value that depends on the surface escape velocity (v esc ) of the larger one. The corresponding is Keplerian velocity and a F is comparable to the Roche radius, The past high-eccentricity phase is not the result of past close encounters, but was instead caused by resonant interactions between moons or pumping by the ring torque.
Here, we investigate the orbital evolution of the mid-sized moons, based on the model of formation from the spreading ring. We employ N-body simulations to compute the evolution and discuss resonance capture and tidal heat production due to eccentricity evolution. We take into account the changes in eccentricities and semi-major axes of the moons according to tides in Saturn and the moons themselves. We also perform runs where we include orbit torques caused by 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 gravitational interactions between the moons, including resonant configurations and eccentricity evolution. The N-body simulations by Salmon & Canup (2017) were the first to investigate the formation scenario of Charnoz et al. (2011) and Crida & Charnoz (2012) in more detail. Although they reproduced the overall mass distribution of these moons, they did not discuss the details of resonance capture/breakup or eccentricity evolution in individual systems, both of which are important to discuss the current orbital architecture of the moons.
We focus our investigation on resonance passing/capture/break-up and tidal heat production with N-body simulations, and therefore we do not include early collisional growth of moonlets from the spreading ring.
The outline of our paper is as follows. Section 2 describes our numerical model, and Section 3 presents our results of the orbital evolution of the moons and discusses the resonance trapping in detail. In Section 4 we estimate the heat energy of Enceladus stored in the course of orbital evolution. Finally, we present our conclusions and a discussion in the last section.

Numerical model
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 moons starting 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 affect other 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(Lainey et al., , 2017, we adopt constant Q p ∼ 2 000-4 000 in our simulations. Because the backward integration does not include gravitational perturbations -such as the resonant perturbations-between the moons, the backward integration cannot be available in the entire space from the birth places of the moons to their current positions. Therefore, we need to perform forward N-body simulations with various initial conditions to investigate which initial conditions result in the current orbital configuration by making individual simulations simple. We will explain the initial conditions in more detail in section 2.4. We have incorporated the semi-major axis expansion and the eccentricity damping due to tidal interaction into the N-body code SyMBA (Duncan et al., 1998). In our simulation, we accelerate tidal evolution by increasing the Love numbers of both the planet and the moons by the same factor, C = 10 3 − 10 4 , to reduce computation time. This "speed-up factor" has been used in many other works (e.g. Malhotra & Dermott, 1990;Showman et al., 1997;Zhang & Nimmo, 2009). The ratio between tidal orbital expansion and eccentricity damping rates is kept the same for different values of the speed-up factor, in order to maintain consistency (see discussions below). Furthermore, the time in all of following figures represents real time that is obtained by simulation time multiplied by the speed-up factor (C). The orbital expansion and eccentricity excitation by the ring torque are implemented in the N-body simulations, and the acceleration with the same speed-up factor for tides is applied. Some aspects in orbital changes by mutual gravitational interactions, including resonant interactions, cannot be accelerated because they are calculated by N-body simulation. We will discuss the effects of the speed-up factor on the probability of resonant trapping.
In our simulations, we find that collisions between moons occur after orbital eccentricities are excited by resonant perturbations or the ring torque. Consequently, the collision velocity (v col ) is usually 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 on SPH simulation results. The critical collision velocity (v cr ), so that a collision between body 1 and 2 (their masses are M 1 and M 2 ) is hit-and-run, is given by and v esc is escape velocity. The fitting parameters are c 1 = 2.43, c 2 = −0.0408, c 3 = 1.86, c 4 = 1.08 and c 5 = 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 the moons against catastrophic disruption. For simplicity, in the case of hit-and-run collisions, we have the moons pass through each other in a softened gravitational potential without collisional energy dissipation. We find that once moons start orbit crossing, they undergo repeated hit-and-run collisions, where we assume that they eventually coalesce, regardless of whether or not we include the energy dissipation during the hit-and-run collision.

Tidal forces
Tidal deformation of the host planet caused by the moons transfers angular momentum from the planetary spin to the moon's orbit. To express the orbital expansion rate in Eq. (1), the tangential force per unit mass is added to the equations of motion, which is given by where r means the orbital radius of the moons. In our simulation, we set the radial component of the planetary tidal force to be zero, because it does not affect the orbital expansion.
A moon's deformation caused by the planet dissipates the moon's orbital kinetic energy, which results in a decrease in the moon's eccentricity (e) and semi-major axis (a). The eccentricity damping timescale is given by (e.g., Murray & Dermott, 1999) 1 For the current orbital elements of Enceladus and k 2m /Q m ∼ 10 −5 , τ e ∼ 8 × 10 8 yrs. Using the eccentricity damping timescale, we add the following tidal force per unit mass caused by the moon's deformation in a simple form to the equations of motion (Kominami & Ida, 2002), where v is the moon's velocity and v K is the local circular Keplerian velocity. Adding the semi-major axis damping associated with the eccentricity damping, the semi-major axis changes as The second term in Eq. (7) corresponds to −2e 2 /τ e in Eq. (7). Because the moon's orbital angular momentum is conserved during the eccentricity damping, automatically caused by f m (Eq. (6)). The second term is dominant when

Saturn's ring torque
In one set of runs (SET2), we include the torque from Saturn's ring following Crida & Charnoz (2012). The timescale is given by where ∆a = a − a F is the separation from F-ring, M ring is Saturn's ring mass, which may be comparable to the moon's mass (Crida & Charnoz, 2012) (see discussion below). Note that this formula is an approximate one for computational simplicity and the actual torque decreases in a discrete manner (e.g., Meyer-Vernet & Sicardy, 1987). Because there is no 1st-order Lindblad resonance beyond 2:1 resonance with the outer edge of the ring, we need to introduce this discreteness at least beyond the 2:1 resonance. We set the ring torque to vanish at a > ∼ 1.59a F . Goldreich & Sari (2003) and Duffell & Chiang (2015) argued that the non co-orbital Lindblad torque excites eccentricity with the timescale as follows, where we used M ring ≃ πΣa 2 . In our simulation, we include the ring torque effects on semi-major axis 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 corotation resonance torque damp e. For a uniform surface density distribution without a gap, the damping is stronger. 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 almost completely 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 protoplanetary disk. Both inner and outer Lindblad torques excite e. Because only Lindblad torque from an inner disk 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's system (Crida & Charnoz, 2012;Salmon & Canup, 2017). As will be shown later, the newly incorporated 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 strongly scattered by the moon for ∆a < ∼ 2 √ 3 r H , where r H is the Hill radius of the moon defined by Ida & Nakazawa, 1989). We should either start integration from ∆a > 2 √ 3 r H or start integration from ∆a = 0 with adding a softening parameter ǫ ∼ 2 √ 3 r H to ∆a in an approach in which collisions are approximated with a softened gravitational potential. Although ǫ ∼ 2 √ the ring edge is smaller than a few × r H (e.g., Petit & Henon, 1986;Ida, 1990), although N-body simulation would be required to prove this argument. Since r H /a ≃ 0.007(M m /10 −6 M p ) 1/3 , we set the initial values of e as e 0 ∼ 0.01 − 0.03.
We assume a constant ring mass (M ring ) throughout each run for simplicity. The assumed ring masses in individual runs are shown in Table 2. The assumption of the constant M ring does not affect our results, because 1) the decrease in M ring would not be significant during the evolution we simulated and 2) the a-e evolution path driven by the ring torque is independent of M ring . The ring mass would largely change when new satellites are formed, while it also changes gradually through viscous 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 already formed. The ring mass should decrease at the time of the formation of Mimas, which we do not include in our simulations. The a-e evolution path driven by the ring torque is independent of the speed-up factor C and M ring , because C is multiplied to both Eq. (9) and Eq. (10) in the simulations and both equations are proportional to M ring . As we will show, the key process to reproduce the current orbital configurations is to avoid capture by mean-motion resonances. The capture can be avoided 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 the ring torque. The migration speed is proportional to M ring . However, we will show that satellites are trapped at a resonance that they encounter at the first time even with C = 10 3 − 10 4 , which means that the capture is not changed even by order of magnitude larger values of M ring in real cases with C = 1. Thus, the assumption of the constant ring mass and the artificial acceleration with C are not critical flaws for the purpose of this paper, while simultaneous evolution of orbits and the ring is important and left for future work.
The phase dominated by the ring torque is determined as follows. Comparing Eq. (9) with the first term of Eq. (7), we find that the ring torque is dominant for where a F is the F-ring radius (a F ≃ 2.4R p ). If we use k 2p /Q p ∼ 10 −4 and M ring /M p ∼ 10 −6 , the solution to Eq. (12) with a = a F + (∆a) crit,a is (∆a) crit,a ≃ 0.82 a F . Because ∆a ≃ 0.59 a F at the 2:1 resonance with the ring edge, the ring torque is dominant in the orbital expansion until it diminishes at the 2:1 resonance as founded by Crida & Charnoz (2012, section 6.3). Since C is multiplied to both Eq. (7) and Eq. (9) in the simulations, (∆a) crit,a is independent of the value of C. It depends on M ring , but only weakly. The eccentricity excitation competes with the tidal eccentricity damping. It dominates over the tidal eccentricity damping, as long as τ e,ring < τ e , that is,

Initial conditions
We performed four sets of simulations (SET1A, SET1B, SET2A, SET2B). Figure 1 suggests that Enceladus formed earlier than Tethys and they underwent orbit crossing, if we neglect the ring torque. Because backward integration is not available back to the state before the orbit crossing, we examine many different initial conditions with a T,0 < a E,0 where the subscripts "T" and "E" represent Tethys and Enceladus, and "0" represents the initial values of the simulation. We call this set of runs "SET1A". We also carry out "SET1B" where Enceladus and Tethys are formed almost simultaneously in a horseshoe orbit, a T,0 ∼ a E,0 . As we will show in the next section, many of these simulations produce a collision between Enceladus and Tethys, and therefore cannot reproduce the current orbital configurations of Enceladus and Tethys.
As we already pointed out, the torque from the ring would affect orbital evolution. Figure 2 shows backward tidal orbital evolution with Q p = 4 000, taking into account the torque from the ring. Because the evolution includes formation of Dione, Tethys, Enceladus and Mimas, we changed the ring mass for this particular plot, while we use a constant ring mass in N-body simulations. The initial mass of the ring is the sum of the masses of Dione, Tethys and 4 × Enceladus, and when one moon is swept out from the ring, the ring loses some of its own mass. Gravitational interactions between the moons are also neglected in this backward integration. For these parameters, it is suggested that Enceladus formed later than Tethys and they can avoid the orbit crossing. Because the orbital evolution by the ring torque is very rapid near the ring, it is likely that the rapid orbital evolution 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:1 resonance with the ring edge. Tethys-Dione migration is always convergent. Because the resonant and secular perturbations among these moons -which are not taken into account in Figure 2-are complicated, we also need to test various initial conditions with a T,0 > a E,0 . 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). We will show that only the runs with both semi-major axis and eccentricity increasing by the ring torque in SET2B potentially reproduce the current orbital configurations of the mid-sized moons. The initial mass of the ring is the sum of the masses of Dione, Tethys and 4×Enceladus, respectively, and when one moon 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 2p , k 2m , Q p and Q m are set to be constant with time for all the moons and we adopt k 2p = 0.34, k 2m = 10 −3 , and Q m = 100 in all runs. We use Q p = 1 700 − 4 000 and the speed-up parameter of C = 10 3 − 10 4 , depending on runs. Figure 3 shows a typical result of orbital evolution of SET1A. In this case, the ring torque is not taken into account (M ring = 0). We adopt Q p = 1 700 and C = 10 4 . Because interactions of Enceladus, Tethys and Dione are essential for the orbital evolution in SET1A, we omit Mimas. The speed-up factor C for semi-major axis expansion and that for the eccentricity damping are 10 4 . Enceladus gets 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.

SET1: Enceladus forms no later than Tethys
We start simulations when Tethys is formed at a T,0 ∼ a F = 1. Enceladus was already formed and has migrated to a E,0 ∼ 1.5a F in agreement with Figure 1. Because Tethys is 5.7 times more massive (Table 1), it catches up with Enceladus (Eq. 1) and Enceladus gets trapped in outer 2:3 mean-motion resonance with Tethys' orbit at t ≃ 0.13 Gyrs. Although the migration is accelerated by a factor of C = 10 4 , Enceladus gets trapped in the 1st order mean-motion resonance that it first encountered. 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 realistic case with C = 1.
As 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 given by Eq. (A.3) in the Appendix A based on Malhotra (1995) shows where we used Eq. (1) and assumed M T ≫ M E . The tidal e-damping is given by Eq. (5) as Equilibrating this excitation and damping with a E /a T = (j + 1/j) 2/3 , the asymptotic value of e E is estimated as Although this analytical estimate includes uncertainty for such high e, the numerical simulation in Fig. 3 actually shows that the eccentricity of Enceladus secularly increases toward a high value, until e E and e T become ∼ 0.18 and ∼ 0.04, respectively, and orbital crossing starts between Enceladus and Tethys at ∼ 0.6 Gyr. Note that the speed-up factor C cancels out in this analytical estimate, suggesting that the equilibrium eccentricity would be similar in a real system with C = 1. After that, Tethys and Enceladus repeat a few tens of hit-and-run collisions, because the collision velocity is 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 cases except 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 adopted the same speed-up factor for tides of the moons and the planet. In some runs, we adopted different values for the moons (C m ) and the planet (C p ). The varies between the runs are only initial orbital angle of each moon. Even in the inwardly scattered case, Tethys and Enceladus collide many times after the inward scattering of Enceladus and before they become isolated. We also performed runs with non-zero energy dissipation at hit-and-run collisions. In those cases, they still repeat collisions and eventually they merge because the collision velocity becomes smaller as the collisions repeat.
In some runs in SET1A, we include the ring torque in the initial condition of a E,0 > a T,0 . We found that the results are similar. Therefore, we conclude that runs in SET1A inevitably end up merging or disrupting Enceladus and the current orbital configuration of the mid-sized moons is never reproduced.
In SET1B, Enceladus and Tethys migrate together in a horseshoe orbit in the early phase. We performed 5 runs of this case and the initial conditions of semi-major axis of Tethys, Dione and Rhea, 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 collisions as in SET1A runs. Therefore, SET1B does not reproduce the current orbital configuration either.

SET2: Enceladus posterior to Tethys
In SET2, we start from the initial conditions with a T,0 > a E,0 . We performed 30 runs in this set: 20 runs considering only semi-major axis expansion by the ring torque (SET2A: subsection 3.2.1 and 3.2.2) and 10 runs considering both semi-major axis expansion and eccentricity excitation by the ring torque (SET2B: subsection 3.2.3).
The current orbital separation between Enceladus and Tethys is slightly smaller than their 3:2 resonance. As Fig. 2 shows, if the ring torque can transfer Enceladus to the orbit inside the 3:2 resonance with Tethys, the current orbital configurations may be reproduced. Even if the ring torque sends Enceladus to the orbit beyond 4:3, 5:4, or a higher-j resonance with Tethys, Enceladus is not trapped at these resonances, because the Tethys-Enceladus migration eventually becomes divergent at a E > 1.59. Then Enceladus divergently passes the closer resonances to Tethys and approaches the 3:2 resonance again. From Eq. (1), in the regions where the ring torque is not effective, the ratio of tidal orbital expansion rate of Tethys to that of Enceladus is described as follows Because M T /M E ≃ 5.7, this ratio is > 1 (divergent migration) for a T /a E < 1.3. Because the 3:2 resonance corresponds to a T /a E ≃ 1.3, the Tethys-Enceladus pair would not be trapped at resonances deeper (smaller a T /a E ) than 3:2 resonance. Therefore, one of the key points in SET2 is whether the ring torque can send Enceladus to an orbit closer to Tethys's orbit than the 3:2 resonance with Tethys.

SET2A: Results with orbital expansion by ring torque
As shown in Eq. (12), the ring torque dominates over the planetary tidal torque for the orbital expansion near the ring. Because the rate of orbital expansion by the ring torque is proportional to (∆a) −3 , with ∆a the distance from the ring edge (Eq. 9), the expansion proceeds very rapidly near the ring edge and the Enceladus-Tethys migration is always initially convergent in the setting of SET2 (a T,0 > a E,0 ). From Eq. (9), in the ring torque dominated region, Because M T /M E ≃ 5.7, the migration becomes divergent when the Enceladus-Tethys pair migrates outward and ∆a T /∆a E decreases to be < ∼ 2. In order to highlight the resonant interaction between Enceladus and Tethys, Figure 4 shows the orbital evolution in the case of only Enceladus and Tethys with the effect of orbital expansion by the ring torque. Enceladus gets trapped in 6:5 resonance at t ≃ 2.0 × 10 −4 Gyrs. This is the 1storder mean-motion resonance that Enceladus meets in the first place during its orbital evolution.
In this run, C = 10 4 . The probability of resonance trapping is higher for slower convergence of the migration. In the real system with C = 1, where migration is much slower, Enceladus should also get trapped in the 1st-order mean-motion resonance that Enceladus meets at the first place.
At t ∼ 3 × 10 −2 Gyrs, ∆a T /∆a E becomes ≃ 2 and the migration becomes divergent. Enceladus leaves the 6:5 resonance with Tethys and its eccentricity decays. At t = 0.06, 0.20 and 0.60, the Enceladus-Tethys pair pass through the 4:3, 7:5 and 3:2 mean-motion resonances, respectively, and their eccentricities are excited. Because Tethys is 5.7 times more massive than Enceladus, the excited eccentricity is much larger for Enceladus than for Tethys. The maximum eccentricity is ∼ 0.08 for Enceladus and ∼ 0.02 for Tethys in this run. Enceladus can store much more heat in its interior than Tethys (see section 4).
Thus, SET2A conditions have the potential to reproduce the current orbital configurations of the mid-sized moons and to account for the high thermal activity of Enceladus. We performed 4 runs of this simulation and the initial conditions of semi-major axis of Tethys, mass of the ring and number of simulations are listed on Table 2. However, runs adding Dione manifest a new problem, as shown below. Gyr), we took x-axis as logarithm scale. Enceladus is swept out from the ring soon after Tethys. Enceladus is captured in a 6:5 mean-motion resonance with Tethys and they migrate outward together until 0.03 Gyrs.

SET2A: Interaction with Dione
Because Dione is more massive than Tethys/Enceladus and Fig. 2 suggests that Dione may not have undergone orbit crossing with other moons, we set Dione at the time of birth of Enceladus at the location predicted by the backward integration in Fig. 2. For this setting, we performed 16 runs. The initial conditions of semi-major axis of Tethys and Dione, speed-up factor C, mass of the 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 run in 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 secularly increased. However, in this case, Tethys rapidly gets trapped in a 2:1 resonance with Dione. The Tethys-Dione migration is always convergent because Dione is only 1.8 times more massive than Tethys, while Dione has sufficiently larger a than Tethys. Since Dione slows down Tethys' migration through 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 and this configuration is stable until the end of the simulation. Due to the resonant migration, the eccentricities of the moons, and that of Tethys in particular, are secularly increased to values that are extremely high in comparison to the values at present. Currently, orbital separation between Tethys and Dione is smaller than 3:2 resonance, which is inconsistent with the trapping in 2:1 resonance obtained by the simulation. The other 9 runs show similar results. Therefore, how to break up the Tethys-Dione's mean-motion resonance is a critical issue.

SET2B: Results with both orbital expansion and eccentricity excitation by ring torque
For the resonant capture, adiabatic convergent migration is required. In addition to that, for the capture, the eccentricity of a moon that encounters a mean-motion resonance with a more massive moon of mass M m must be smaller than a critical value (M m /M p ∼ 2 × 10 −6 for Dione), given by (Malhotra, 1993) e crit ≃ 1.58 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 runs of simulations in this set and the initial conditions of semi-major axis of Tethys and Dione, dissipation factor of Saturn 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 = 4M E and e E,0 = 0.01 for Enceladus. In the beginning of Fig. 6 (on the left of this figure), we accelerated the tidal orbital evolution with C = 10 3 , and adopted k 2m /Q m as 10 −5 . Equation (14) predicts that the eccentricity increases until a reaches (∆a) crit,e ≃ 0.3 for Enceladus, which is consistent with the numerical result. In the numerical result, e E is already excited up to ∼ 0.04 at the timing of Tethys-Dione 2:1 resonance passing. Because e E is well excited beyond e crit , Enceladus passes through 3:2, 4:3 and closer to 1st order resonances with Tethys, although the Enceladus-Tethys migration is convergent until Enceladus reaches a E ∼ 1.59 and the ring torque decays.
Secular perturbation from Enceladus with the relatively high eccentricity enhances Tethys' eccentricity up to e T ∼ 0.01 ∼ e crit . Tethys encounters a 2:1 resonance with Dione at t ≃ 0.04 Gyrs.
At the resonance passage, a relatively large amount of angular momentum is exchanged between Tethys and Dione and also between Tethys and Enceladus. Owing to e T ∼ e crit , Tethys successfully avoided getting trapped in a 2:1 resonance with Dione.
In Fig. 6, the tidal orbital evolution is accelerated from C = 10 3 (the left panels) to C = 10 4 (the right panels) to follow the whole orbital evolution. After a E exceeds 1.59 at 0.6 Gyr, the Enceladus-Tethys migration becomes divergent and it becomes impossible for them to be trapped in 1st-order resonances between them. The calculation of the following evolution shows that the final a T and a D are consistent with their current semi-major axes (a T = 2.11 and a D = 2.70).
However, both Tethys' and Dione's eccentricity (e T and e D , respectively) are enhanced too much by the resonance in the simulation. Indeed, the Tethys-Dione pair gets trapped in a 5:3 resonance at t ≃ 0.45 Gyrs and the resonance configuration remains stable, while the current Tethys-Dione separation is slightly smaller than the 3:2 resonance. Immediately after that, Enceladus gets trapped in a 2nd-order 9:7 mean-motion resonance with Tethys. These resonances are not consistent with the current orbit.
In this case, a mechanism is needed to kick the Tethys-Dione pair out of their resonance, such as the impact that created the Odysseus crater on Tethys (Zhang & Nimmo, 2012). If we consider the growth of the moons by merging as they migrate outward, it is, in principle, possible that the moons avoid the resonance capture. However, N-body simulations of collisional growth of the mid-sized moons by Salmon & Canup (2017) showed that capturing at Tethys-Dione's 2:1 and 3:2 Dione (blue) with e E,0 = 0.01. In the left panels, the speed-up factor of C = 10 3 is used. Tethys passes 2:1 resonance with Dione at t ≃ 0.04 Gyrs. The results in the right panels start from the end state of the evolution in the left panels and use C = 10 4 . 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 ≃ 0.45 Gyrs and immediately afterward, Enceladus is trapped into a 9:7 resonance with Tethys (the right panels). Gyrs, as in the results of Fig. 6. The middle panels start from the end state of the evolution in the left panels and use C = 10 4 . In this case, Tethys and Enceladus collide with each other at ∼ 0.18 Gyr. However, when the eccentricity of Enceladus becomes smaller, they never collide and it can reproduce the orbital configuration (the right panels).
At the end of the right panels (∼ 1.1 Gyr), Enceladus and Dione get trapped in 2:1 resonance, which is the current state.
resonances is still a great difficulty (their Figs. 11 and 12), although they did not discuss this issue.
In Fig. 7, we set e E,0 = 0.03 and other parameters were the same as those in the results of Fig. 6. With the three times higher e E,0 , the eccentricity of Enceladus increases up to e E ∼ 0.15.
In addition, Tethys' eccentricity is also enhanced, because the effect of the secular perturbation is proportional to the perturber's eccentricity. At t ≃ 0.04 Gyrs, as e T is enhanced well beyond e crit , Tethys passes through the Tethys-Dione 2:1 resonance and angular momentum is transported between the moons (the left panels). The eccentricities of all the moons are excited and oscillate substantially from t ≃ 0.04 Gyrs to the point of a E ≃ 1.59, where the excitation of eccentricity and semi-major axis by the ring torque for Enceladus decays.
In the middle panels, the following orbital evolution is calculated with C = 10 4 . As e D is excited to ∼ 0.1 − 0.15, Enceladus eventually collides with Tethys at t = 0.18 Gyrs. We performed 11 runs with orbital expansion and eccentricity excitation by the ring torque. In 2 runs with e E,0 < ∼ 0.03, the orbital evolution is similar to Fig. 6 and Tethys is captured in the Tethys-Dione 3:2 resonance.
In other 8 runs with e E,0 > ∼ 0.04, e E is excited so much that Enceladus undergoes collisions or close scattering with Tethys. Therefore, it is difficult for Tethys to successfully pass the 3:2 resonance with Dione and simultaneously avoid a collision with Enceladus, as long as only Enceladus, Tethys and Dione are integrated.
We consider a hypothetical case in which the eccentricity of Enceladus is smaller than in the middle panels. In the right panels in Fig. 7, we artificially decreased e E by a factor of 2 from the final state of the left panels and calculated the following orbital evolution. Although we do not specify the cause of the decrease, the ring mass decrease due to the birth of Mimas could be responsible for it. In this case, Enceladus and Tethys never collide. But, the eccentricity of Tethys excited by secular perturbations from Enceladus is still large enough for Tethys to avoid getting trapped in the 3:2 resonance with Dione. After passing the 3:2 resonance, both e T and e D are damped. Enceladus is 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 increases to an equilibrium value, e E ∼ 0.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 the ring is still massive enough after Mimas' formation, the torque from the ring can be transported to Enceladus, Tethys and Dione through Mimas, because Mimas is currently located within a 2:1 resonance with the ring edge and should have suffered the ring torque throughout its entire orbital evolution. Figure 1 suggests that Mimas-Enceladus encounters their 3:2 resonances at a similar time as the trapping of the Tethys-Dione pair at their 3:2 resonance. It is very likely that the ring torque pumps up Mimas' eccentricity to a value larger than e crit and the Mimas-Enceladus pair avoids the trapping at their 3:2 resonance. However, interactions among the four moons with the ring torque are complicated. The ring mass should also change at the formation of Mimas. Because these investigations require much more parameter surveys, we leave them to future study.

Heat Flux
As we have shown, the moons would have undergone a high eccentricity phase in the past during orbital evolution. As we show below, the heat generated during the high eccentricity phase can be stored in the interior and the current high heat flux can reflect the stored heat (the current heat generation is not balanced with the surface heat flux). From the numerical simulations, here we calculate the stored heat energy for each moon, where H is given by Eq. (2).
Although Enceladus would have a subsurface ocean at present, most parts of bulk Enceladus would be in a solid phase. For simplicity, we assume that conduction is a major heat transfer mechanism in the interior of the mid-sized moons. We estimate the conduction timescale very x) + 900x]} Jm −3 /K, respectively. The thermal diffusion coefficient of Enceladus is then where ζ = (2 + x)/[(1 + 2x)(2 − 1.1x)] ∼ 1. The thermal conduction timescale for a moon with a physical surface radius of R is Gyrs.
(23)    Because τ cond may be longer than the age of Enceladus in the low Q p model, we simply assume that most of the heat energy generated during the high eccentricity phase is still stored in the interior.
Then, the total heat flux is given by 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) 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 and the total heat flux L for individual simulations are listed in Table 3. In the case corresponding to Fig. 4, the tidal heat is consistent with the current observed value. However, in the results in Figs. 5 and 6, Tethys is captured by a resonance with Dione and subsequently has higher tidal heat generation than Enceladus. In the case of Fig. 7, tidal heat of Enceladus is not inconsistent with the observed value which is suggested by Howett et al. (2011);Spencer & Nimmo (2013). While the heat flux of other moons is smaller than Enceladus, the heat flux of Tethys predicted by our simulation 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 not only in final orbital configurations but also in the heat flux.

Conclusion and discussion
Through N-body simulations, we have numerically investigated the orbital evolution of Saturn's mid-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 model of 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 with various initial conditions at two different settings.
If the ring torque is sufficiently weak, Enceladus must be formed prior to Tethys and scattered inward across Tethys' orbit. In this set of runs (SET1), we found that Enceladus is always trapped in the outer 1st-order mean-motion resonance with Tethys due to the rapid migration of Tethys which is more massive than Enceladus. Such resonance trapping is inevitable because of the tidal migration timescale of these moons. The eccentricity of Enceladus is secularly increased by the resonant migration. When the eccentricity becomes large enough, Enceladus undergoes collisions with 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 and Tethys 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 Enceladus need not cross Tethys' orbit. We also performed many simulations with many different initial conditions in the setting in which Enceladus is initially located inside Tethys' orbit (SET2). Because the ring torque is very strong near the outer edge of the ring, Enceladus-Tethys migration is convergent as long as Enceladus' orbital radius is smaller than the 2:1 resonance with the ring edge. As a result, Enceladus is always trapped in the 1st-order mean-motion resonance with Tethys that Enceladus meets in the first place. After the trapping, Enceladus' eccentricity is secularly increased to undergo repeated hit-and-run collisions as in SET1 runs. However, for some range of initial 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 eccentricity decays via the tidal damping, which results in the final orbits of Enceladus and Tethys being consistent with the current ones.
However, Tethys is trapped in the 2:1 or 3:2 mean-motion resonance with Dione in this set of simulations, because their migration is always convergent and adiabatic. It is very difficult to reproduce their current orbital separation closer than the 3:2 resonance relation between Dione and Tethys. Fuller et al. (2016) considered time-dependent Q p , in which Q p > ∼ 10000 in the initial phase and it decreases down to Q p ∼ O(1000) after resonant locking between the orbital frequency and Saturn's oscillation mode. In this case, moons can be formed in an extended circumplanetary disk and orbit crossing of the moons does not occur. Even in this model, the migration between Dione and Tethys is usually convergent and adiabatic and it is difficult for them to avoid becoming trapped in the 3:2 resonance.
We found that, if eccentricity excitation by the ring torque is effective, it has the potential to solve the problem. This excitation is effective only in the regions close to the outer edge of the ring. Owing to the excited eccentricity, Enceladus can easily pass through the resonances with Tethys. Enceladus' eccentricity could be mostly enhanced by the ring torque, but not by the resonant perturbations. The modest eccentricity of Tethys raised by the secular perturbation from Enceladus with relatively high eccentricity breaks the 2:1 resonance between Tethys and Dione. In part of the range of Enceladus' eccentricity, Tethys can also pass the 3:2 resonance with Dione. As the distance between Enceladus and the ring increases, tidal eccentricity damping dominates over the 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 a constant ring mass, we only found the orbital evolution path if Enceladus' eccentricity is artificially decreased after the Tethys-Dione 2:1 resonance is passed; this is because Enceladus' eccentricity necessary to bypass the Tethys-Dione 3:2 resonance is so large that it results in collisions between Enceladus and Tethys afterward. Although we assumed that the ring mass is constant to focus ourselves on orbital evolution, it must evolve with time. The ring mass decrease due to Mimas formation could lower Enceladus' eccentricity to avoid the collision between Enceladus and Tethys.s The dynamical effect of Mimas, the smallest moon among the mid-sized moons, could also play an important role, although we have not explored the effect. The assumption of constant k 2m /Q m for all the moons is also too simple. The parameter survey taking these effects into account is left for a separate paper.
We also estimated the heat flux of individual moons, which is caused by thermal energy stored in past periods of high eccentricity as a result of resonant interactions and ring torque. Unless Enceladus and Tethys are captured by the mean-motion resonance, the heat generation is the highest for Enceladus. The orbital evolution with the eccentricity excitation by the ring torque can produce the heat flux that is comparable to or slightly smaller than the observationally inferred value of Enceladus. Tidal heating due to the high eccentricity events may make the moons more convective and dissipative, which may significantly increase k 2m /Q m . We will also address these issues in a separate paper.
In conclusion, if we take into account of the orbital expansion and eccentricity excitation by ring torque, there will be one possible pass to solve the problem of the current heat budged on Enceladus and the resonance capture from the birth to the current orbit. where we used k 2p ∼ 0.3, Q p ∼ 3000, k 2m ∼ 10 −3 and Q m ∼ 10 2 . These eccentricities are one order higher than the current values, one suggestion is raised that their eccentricities are now on the way to the equilibrium.
Note that the derivation for the equilibrium eccentricity here is simplified. More rigorous derivations with Lagrange equations and detailed resonant properties are found in the past literatures (e.g., Zhang & Nimmo, 2009). While the numerical factors differ from the rigorous treatment by a factor of up to a few, the dependence on k 2,p /Q p and k 2,m /Q m are reproduced by the simple derivation here and the difference in the numerical factors does not affect the discussions in this paper.