© ESO 2020

1. Introduction

The orbital configuration of Saturn’s mid-sized moons, Mimas, Enceladus, Tethys, Dione, and Rhea (from inner to outer orbits), is puzzling: they are trapped in mean-motion resonances for every-other pairs (Mimas-Tethys 4:2 and Enceladus-Dione 2:1), but not for adjacent pairs. The observed current fast tidal orbital expansion rate (Lainey et al. 2012, 2017) and observations of the rings by Cassini suggest late formation of the satellites (see, e.g., Ida 2019) such as the formation model from a hypothetical ancient massive rings (Charnoz et al. 2011; Crida & Charnoz 2012), which we refer to as the disk. We note, however, that a possible slower tidal orbital expansion rate in the past, before resonance locking between the orbital frequency and the planetary internal oscillation mode, could allow the satellite formation in the circumplanetary disk 4.5 G years ago (Lainey et al. 2020).

In the model of the formation from the disk, satellites were formed one after another at the disk outer edge and the outward orbital migrations of adjacent pairs of satellites due to planetary tides are generally convergent. In that case, the satellite pairs are usually captured into a mutual first-order mean-motion resonance, which is inconsistent with the current orbital configuration (e.g., Nakajima et al. 2019). The avoidance of this resonance capture requires moderate orbital eccentricity of the satellites or fast orbital migration on a timescale shorter than the resonant libration period (Malhotra 1996).

Recent Cassini observations determined the current ring mass to be M_disk = (1.54 ± 0.49)×10¹⁹ ≃ 0.4× Mimas mass (Iess et al. 2019). The rings are still undergoing viscous spreading and were probably much more massive in the past (Salmon et al. 2010). Crida & Charnoz (2012) suggested that the satellite-disk (rings) interaction is more effective for the orbital migration than Saturn’s tide until the satellite reaches 2:1 resonance with the disk’s outer edge, beyond which the disk torque would quickly decay. They applied the theoretical model for a planet in a gap of a protoplanetary disk (e.g., Lin & Papaloizou 1986) to estimate the migration rate of the satellite. The satellite-disk interactions can also excite the orbital eccentricity of the satellite (Goldreich & Sari 2003; Duffell & Chiang 2015). The eccentricity excitation may be much faster than the eccentricity damping by the satellite’s tide near the edge of the disk, as we show in Sect. 2. If the eccentricity is excited beyond a critical value, the satellites can avoid the resonance capture to reach the current orbital configuration (Nakajima et al. 2019).

Previous studies of planet-disk interactions usually assumed protoplanetary gas disks that are stable against self-gravitational instability, while the rings are often in a marginally unstable state (e.g., Salo 1995; Daisaka et al. 2001). Thus, it is important to investigate the interactions of a satellite and a marginally unstable self-gravitating particle disk (rings) via high-resolution N-body simulations.

Hyodo et al. (2015) performed high-resolution N-body simulations (N = 3 × 10⁴ − 5 × 10⁴) of the formation of satellites from a disk with mass M_disk ≃ (0.01–0.06) M_p (M_p is the planet mass) to find the dependence of the forming satellite mass on the initial disk mass. On the other hand, they were not concerned with the orbital evolution of satellites and disk structures.

Here we focus on the detailed evolution of the disk structures and the satellite’s orbit. In our study we perform high-resolution (N ∼ 10⁵) N-body simulations of particle disk evolution due to the mutual gravitational interactions and inelastic collisions of the particles and the disk’s interactions with a satellite in an orbit exterior to the disk. Our aim is to investigate the orbital evolution of the satellite.

2. Methods

We use a new N-body simulation code called GPLUM (Ishigaki, in prep.) that adopts the particle-particle particle-tree scheme (P³T; Oshino et al. 2011) for planetary formation. The P³T scheme uses the fourth-order Hermite integrator to calculate gravitational interactions between particles within a cutoff radius and the Barnes-Hut tree method for gravity from particles beyond the cutoff (Iwasawa et al. 2016), which guarantees higher-order integrations for close interactions and fast integrations for perturbations from a large number of distant particles simultaneously. GPLUM adopts an individual cutoff radius scheme for individual particles, depending on their mass and distance from the central star, resulting in a significant increase in the speed of calculations, while keeping the accuracy.

We follow the orbital evolution of the satellite and the disk particles by the gravitational interactions, inelastic collisions between the particles, and the accretion of the particles onto the satellite and the host planet. When physical sizes overlap, we regard that a collision occurs. For collisions between the disk particles, we apply inelastic collisions with the normal restitution coefficient ϵ_n = 0.1 and the tangential restitution coefficient ϵ_t = 1 (free-slip condition). When a particle collides with the satellite at an orbital radius larger than the Roche limit radius (Eq. (1)), the collision results in gravitational binding of the particle and the satellite. Because the satellite does not reenter the Roche limit in the results we show here, we make a merged body from the satellite and the particle, keeping their total mass and momentum.

The particles are initially distributed from the physical radius of the planet (R_p) to the Roche limit radius (denoted by a_R), which is given by

(1)

where ρ and ρ_p are the bulk densities of the disk particles and the planet, respectively. In this paper we assume ρ = 0.9 g cm⁻³ and ρ_p = 0.7 g cm⁻³, so that a_R ≃ 2.26 R_p. The initial surface density of the particles follows Σ(r)∝r⁻³, and their total mass is ∼(10⁻³–10⁻²) M_p. We use 8 × 10⁴ − 1.2 × 10⁵ particles with equal mass of M ∼ (10⁻⁸–10⁻⁷) M_p for the disk. Initially, the particles have circular orbits with low inclinations, following a normal distribution of ⟨e²⟩^1/2 = 2⟨i²⟩^1/2 ∼ R/r ∼ (2 − 4)×10⁻³, where R is the particle physical radius. Because of the inelastic collisions and self-gravity of the disk particles, they are quickly relaxed to quasi-equilibrium values, which are also ∼R/r (e is a few times larger than R/r probably due to the scattering by the satellite). The satellite with a mass M_s ∼ 10⁻³M_p is placed outside the Roche limit.

We do not include the outward orbital migration due to the planetary tide and eccentricity damping due to the satellite tide because they are negligible compared with the migration due to the satellite-disk interactions at the radius inside the 2:1 resonance with the disk’s outer edge. The tidal e-damping and a-expansion timescales are (Q_s/k_2s)(M_s/M_p)(a_s/R_m)⁵Ω⁻¹ and (Q_s/k_2s)](M_p/M_s)^1/3τ_e, tide (e.g., Charnoz et al. 2011). For the tidal parameters for the satellite Q_s/k_2s ∼ 10⁵, M_p/M_s ∼ 10⁶, and the satellite orbital radius a_s ∼ a_R ≃ 2.26R_p, years. For the planet tidal parameter Q_p/k_2p ∼ 10³–10⁵, –10⁹ years. As we show in Sect. 3, the a-expansion timescale due to the satellite-disk interactions is τ_a, disk ∼ (π²/51)(M_p/M_disk)³(M_s/M_p)(a_R/a_s)^1/2Ω⁻¹ ∼ 7 × 10³Ω⁻¹ ∼ 2 years for a realistic case with M_disk/M_p ∼ 3 × 10⁻⁴ and M_s/M_p ∼ 10⁻⁶ (see Sect. 4). Since near the disk outer edge, the assumption to neglect the tidal force is justified in our simulation.

Table 1 shows the parameter sets of the runs with the different initial disk mass (M_disk) and satellite mass (M_p). The disk particles have equal masses. We use the satellite masses that are much higher than Saturn’s mid-size moons. We derive semi-analytical formulas from the results of N-body simulations to clarify intrinsic physics in this system. Applying the derived mass scaling law for the realistic masses of Saturn’s mid-size moons, we discuss the possibility to avoid the resonance trapping.

3. Simulation results

3.1. Ring structures

Figure 1c is a snapshot at t = 1.47 × 10³T_Kep of RUN 1, where T_Kep is the Keplerian period at r = a_R. The figure shows that two distinct structures are superposed: (1) the dense wake structure with small wavelengths caused by the disk self-gravity (e.g., Salo 1995; Daisaka et al. 2001; Takeda & Ida 2001) and (2) the m = 2 spiral arms produced by the Lindblad resonance torque from the satellite. The radial wavelength and wavenumber of the self-gravity wakes are estimated as (Takeda & Ida 2001)

(2)

(3)

Fig. 1. Time evolution of the system of RUN1: (a) t = 0, (b) t = 29.4 TKep, (c) t = 1.47 × 103 TKep, and (d) t = 2.94 × 104 TKep, where TKep is the Keplerian period at r = aR. The green point is the outer satellite and the purple dots show the disk particles. The inner and outer circles in black solid lines represent the planetary surface and the Roche limit (r = aR).

where M_disk is the total disk mass and we assumed the pitch angle ∼π/4 to estimate m_self. These estimates are consistent with the result in Fig. 1c. The Lindblad torque exerted by the satellite makes the spiral arms on the disk (e.g., Goldreich & Tremaine 1982). The wavenumber of the spiral arms is given by

(4)

where m_res ≠ 1, and Ω_s and Ω are the orbital frequencies of the satellite and the disk. Figure 1c clearly shows m_self = 2 spiral arms. With a_s ∼ 3R_p and Ω at ∼2R_p, Eq. (4) predicts m_res ∼ 2^−3/2/(2^−3/2 − 3^−3/2)∼2, which agrees with Fig. 1c.

Our N-body simulation simultaneously show the short-wavelength wakes due to the self-gravity, which were often shown in the local shearing sheet simulations (e.g., Salo 1995; Daisaka et al. 2001), and the m_res = 2 global spiral arms structure. While Hyodo et al. (2015) also performed N-body simulation with N = 3 × 10⁴ − 5 × 10⁴, they did not clearly show these two distinct structures, because they used ten times higher M_disk (accordingly, ten times fewer m_res) and because their simulations often had multiple clumps.

3.2. Time evolution of disk structures and the satellite’s orbit

Figures 1a–d show the time evolution of the disk structures. We initially set the satellite near the disk’s outer edge (Fig. 1a). In the disk, the dense wakes quickly emerge due to the combined effect of self-gravity and inelastic collisions. The satellite creates m_res = 3 spiral arms by the Lindblad torque and scatters or accretes nearby disk particles to open a gap with the half width ∼3r_Hill, where r_Hill is the satellite Hill’s radius defined by r_Hill = (M_s/3M_p)^1/3a_s (Fig. 1b). However, the number of particles scattered outside of the Roche radius is only ∼2000 at this time. The disk mass loss is mostly caused by accretion onto the planet due to the viscous spreading. At the initial satellite location (a_s ≃ a_R), Eq. (4) at the disk outer edge (r ∼ a_R − 3r_Hill) is

(5)

This is consistent with the m_res = 3 spiral mode in Fig. 1b.

After the gap opening, the satellite migrates outward, and m_res decreases from 3 to 2 (Fig. 1c). When the satellite orbit expands beyond the 2:1 resonance with the disk’s outer edge, m_res becomes much smaller than 2 and the spiral arms disappears (Fig. 1d). Because the theoretically predicted effective viscosity is ∝Σ² (Eq. (6)), the r-gradient of Σ is quickly flattened. The self-gravity wakes become fainter (Eq. (2)) through the loss of M_disk and through the spiral arms decay (Fig. 1d). Consequently, the satellite’s orbital migration slows down. In the next subsection, we quantitatively discuss the satellite migration rate.

3.3. Satellite orbital evolution

Figure 2 shows the time evolution of the satellite’s semimajor axis (a_s) and the total disk mass (M_disk) obtained by our N-body simulation (RUN 1). The changes in a_s and M_disk become slower with t, which are theoretically explained as follows. Through N-body simulations, Daisaka et al. (2001) found that the effective viscosity of a self-gravitating disk is given by

(6)

Fig. 2. Time evolution of the satellite’s semimajor axis (as; upper panel) and the disk mass (Mdisk; lower panel). The units of the semimajor axis and the time are the Roche limit radius aR and the Keplerian period at aR, respectively. The red curves are the results of RUN1 and the blue curve in the lower panel is the analytical estimation given by Eq. (8).

where the subscript “_R” represents the values at r ≃ a_R, , and represents the effect of finite physical size of the particles. Because we are concerned with the outer disk region, we adopted r ≃ 0.8a_R. The rate of disk accretion onto the planet is

(7)

Integrating this equation, we predict the explicit time evolution of the disk mass,

(8)

which reproduces the N-body simulation result (lower panel of Fig. 2).

Using Eq. (8), we show that the migration is regulated by the self-gravity wakes, but not by the Lindblad resonance torque (the spiral arms) induced by the satellite. The migration rate due to the (one-sided) Lindblad resonance is (e.g., Lin & Papaloizou 1986; Crida & Charnoz 2012)

(9)

where Δa_s = a_s − a_R. The migration rate by the self-gravity wakes is evaluated as follows. When the disk viscous spreading beyond the Roche limit is prevented by the satellite’s perturbations, the angular momentum flux in the disk (∼3πΣν r²Ω) is transferred from the disk outer edge to the satellite’s orbit. In this case,

(10)

Substituting Eq. (6) and into Eq. (10), we obtain

(11)

In Fig. 3, da_s/dt from each run of our N-body simulations is compared with analytical estimations. In the analytical estimations, M_s and M_disk obtained by the N-body simulations at each a_s are substituted to Eqs. (9) and (11) to calculate (da_s/dt)_res and (da_s/dt)_self. This figure shows that the results of N-body simulations fit (da_s/dt)_self. In the vicinity of the disk’s edge, (da_s/dt)_res dominates over (da_s/dt)_self. The theoretical prediction for (da_s/(dt)_res assumes a non-self-gravitating disk with modest viscosity. The spiral arms may be weakened by the relatively strong diffusion due to self-gravity wakes. Because the self-gravity wake torque is independent of the satellite locations, it dominates over the Lindblad torque, which is very sensitive to the distance from the disk outer edge.

Fig. 3. Orbital expansion rate, das/dt, as a function of as in RUN1 to RUN5. The red curves are the results of individual N-body simulations. The blue and green curves are the theoretical predictions of (das/dt)res and (das/dt)self given respectively by Eqs. (9) and (11) using Ms and Mdisk obtained by the N-body simulations at individual as in each run.

In these runs we adopted M_s ∼ 10⁻³M_p, while the masses of the actual mid-sized moons are M_s ∼ (10⁻⁷–10⁻⁶)M_p. Hyodo et al. (2015) showed through N-body simulations that M_s/M_p ∼ 10(M_disk(0)/M_p)² for M_s ∼ 10⁻³M_p. Although Crida & Charnoz (2012) proposed M_s/M_p ∝ (M_disk(0)/M_p)³ for the smaller value of M_s/M_p, generated clumps would quickly coagulate. Here we use the Hyodo et al. (2015) relation to discuss the cases of M_s ∼ (10⁻⁷–10⁻⁶)M_p. As is shown in Sect. 4, the high migration rate induced by the self-gravity wake torque, which has been overlooked in the past studies, plays an important role in avoidance of the mean-motion resonant capture between adjacent satellites.

Figure 4 shows the eccentricity evolution of the satellite in the individual runs. The eccentricity is excited only in the early phase, t ≲ (10⁴–10⁵) T_Kep, when the Lindblad torque may not be negligible compared to the self-gravity wake torque. As we discuss in Sect. 4, the excited eccentricity in the early phase (e ∼ 0.01) is marginal for the condition to avoid a mean-resonance trapping.

Fig. 4. Eccentricity evolution of the massive satellite during the migration in each run of our simulations.

4. Resonance capture probability

The probability of the mean-motion resonance capture depends on the eccentricity and migration rate of the satellites (Dermott et al. 1988; Malhotra 1993). Here, we consider the avoidance of Tethys-Dione 3:2 resonance tapping as an example. Dione’s orbit is currently located beyond that of Tethys with the separation within 3:2 mean-motion resonance. Because tidal migration rate is a strong function of orbital radius, while the mass difference between Dione and Tethys is within a factor of 2 (M_Dione/M_Saturn ≃ 1.94 × 10⁻⁶, M_Tethys/M_Saturn ≃ 1.09 × 10⁻⁶), their tidal migration are convergent unless the planetary tidal parameter Q_p is much lower (much more dissipative) for Dione. To avoid the trapping into their 3:2 mean-motion resonance, a large enough eccentricity and/or fast enough convergence of their orbits is required.

The critical eccentricity, beyond which the j + 1 : j resonance trapping is inhibited, is given by (Malhotra 1996)

(12)

Nakajima et al. (2019) pointed out the possibility to avoid Tethys-Dione 3:2 resonance trapping by Enceladus’ eccentricity excitation. If Enceladus’ eccentricity is excited enough by the Lindblad torque from the disk, Tethys’ eccentricity can also be excited to ≳e_crit by the secular perturbations from Enceladus.

The other possibility is fast orbital migration on a timescale shorter than the resonant libration timescale. As we have shown, the migration is significantly faster if we consider the satellite-disk interactions. According to Ogihara & Kobayashi (2013), resonance trapping is avoided for Tethys-Dione 3:2 resonance (j = 2) if

(13)

where a_s is the semimajor axis of the inner satelliite (Tethys), α is the semimajor axis ratio (≃0.763), and f(α)∼ − 1.55/α (see Murray & Dermott 1999, Table 8.5). With M_s/M_p ∼ 10⁻⁶. τ_a, crit ∼ 8.0 × 10⁶T_Kep.

If we consider the Tethys-mass satellite (M_s ∼ 10⁻⁶M_p), the disk mass may be M_disk ∼ 3 × 10⁻⁴M_p (Hyodo et al. 2015). The orbital migration timescale of the self-gravity wake is estimated to be τ_a, self ∼ 1.2 × 10³(M_disk/3 × 10⁻⁴M_p)⁻³T_Kep (Eq. (11)), which is shorter than t_a,crit by more than three orders of magnitude. Although the migration by the self-gravity wake torque could be weaker along with the decrease in M_disk, the fast migration potentially prevents the resonance capture of the Dione-Tethys pair.

5. Conclusions

In order to investigate the gravitational interactions between Saturn’s mid-sized moons and hypothetical ancient massive rings and the associated orbital evolution of the moons, we performed global high-resolution N-body simulations (N ∼ 10⁵) of a self-gravitating particle disk interacting with a single satellite, taking account of gravitational forces among all the disk particles and the satellite and inelastic collisions between the particles. Our simulations show that the dense short-wavelength wake structure and m = 2 or 3 global spiral arms simultaneously develop in the disk and are produced by the disk self-gravity and Lindblad torque from the satellite, respectively. These structures transfer the angular momentum of the disk to the satellite and regulate the early phase of orbital evolution of the satellite. We found that the orbital migrations of the satellites are determined by the self-gravity wakes. The literature assumes that the Lindblad torque regulates the migrations because they considered the self-gravity wakes only for the source of the disk diffusion.

In this paper, we focused on investigating the detailed dynamics of gravitational interactions between a circumplanetary particle disk and a satellite, and derived the semi-analytical formulas for the satellite’s migrating rate. While the simulations used a much more massive satellite than Saturn’s mid-sized moons due to the simulation limitation, we extrapolated the formulas to realistic satellite masses to find that the migration is fast enough to avoid the resonance capture of adjacent moons on the way to the current orbital configuration of Saturn’s mid-size moons. To confirm this conclusion, simulations with much higher resolution and with multiple satellites are needed, which is left for future study.

Acknowledgments

We thank Takaaki Takeda for helpful and useful comments. This research was supported by JSPS Grants-in-Aid for Scientific Research (# JP19J12542) and MEXT “Exploratory Challenge on Post-K computer” hp190143. Numerical computations were in part carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.

References

All Tables

All Figures

	Fig. 1. Time evolution of the system of RUN1: (a) t = 0, (b) t = 29.4 TKep, (c) t = 1.47 × 103 TKep, and (d) t = 2.94 × 104 TKep, where TKep is the Keplerian period at r = aR. The green point is the outer satellite and the purple dots show the disk particles. The inner and outer circles in black solid lines represent the planetary surface and the Roche limit (r = aR).
In the text

	Fig. 2. Time evolution of the satellite’s semimajor axis (as; upper panel) and the disk mass (Mdisk; lower panel). The units of the semimajor axis and the time are the Roche limit radius aR and the Keplerian period at aR, respectively. The red curves are the results of RUN1 and the blue curve in the lower panel is the analytical estimation given by Eq. (8).
In the text

	Fig. 3. Orbital expansion rate, das/dt, as a function of as in RUN1 to RUN5. The red curves are the results of individual N-body simulations. The blue and green curves are the theoretical predictions of (das/dt)res and (das/dt)self given respectively by Eqs. (9) and (11) using Ms and Mdisk obtained by the N-body simulations at individual as in each run.
In the text

	Fig. 4. Eccentricity evolution of the massive satellite during the migration in each run of our simulations.
In the text