Free Access
Issue
A&A
Volume 640, August 2020
Article Number L15
Number of page(s) 6
Section Letters to the Editor
DOI https://doi.org/10.1051/0004-6361/202038743
Published online 14 August 2020

© 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 Mdisk = (1.54 ± 0.49)×1019 ≃ 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 × 104 − 5 × 104) of the formation of satellites from a disk with mass Mdisk ≃ (0.01–0.06) Mp (Mp 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 ∼ 105) 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 (P3T; Oshino et al. 2011) for planetary formation. The P3T 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 (Rp) to the Roche limit radius (denoted by aR), 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−3 and ρp = 0.7 g cm−3, so that aR ≃ 2.26 Rp. The initial surface density of the particles follows Σ(r)∝r−3, and their total mass is ∼(10−3–10−2) Mp. We use 8 × 104 − 1.2 × 105 particles with equal mass of M ∼ (10−8–10−7) Mp for the disk. Initially, the particles have circular orbits with low inclinations, following a normal distribution of ⟨e21/2 = 2⟨i21/2 ∼ R/r ∼ (2 − 4)×10−3, 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 Ms ∼ 10−3Mp 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 (Qs/k2s)(Ms/Mp)(as/Rm)5Ω−1 and (Qs/k2s)](Mp/Ms)1/3τe, tide (e.g., Charnoz et al. 2011). For the tidal parameters for the satellite Qs/k2s ∼ 105, Mp/Ms ∼ 106, and the satellite orbital radius as ∼ aR ≃ 2.26Rp, years. For the planet tidal parameter Qp/k2p ∼ 103–105, –109 years. As we show in Sect. 3, the a-expansion timescale due to the satellite-disk interactions is τa, disk ∼ (π2/51)(Mp/Mdisk)3(Ms/Mp)(aR/as)1/2Ω−1 ∼ 7 × 103Ω−1 ∼ 2 years for a realistic case with Mdisk/Mp ∼ 3 × 10−4 and Ms/Mp ∼ 10−6 (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 (Mdisk) and satellite mass (Mp). 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.

Table 1.

Parameter sets of our simulations.

3. Simulation results

3.1. Ring structures

Figure 1c is a snapshot at t = 1.47 × 103TKep of RUN 1, where TKep is the Keplerian period at r = aR. 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)

thumbnail Fig. 1.

Time evolution of the system of RUN1: (a) t = 0, (b) t = 29.4 TKep, (c) t = 1.47 × 103TKep, and (d) t = 2.94 × 104TKep, 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 Mdisk is the total disk mass and we assumed the pitch angle ∼π/4 to estimate mself. 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 mres ≠ 1, and Ωs and Ω are the orbital frequencies of the satellite and the disk. Figure 1c clearly shows mself = 2 spiral arms. With as ∼ 3Rp and Ω at ∼2Rp, Eq. (4) predicts mres ∼ 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 mres = 2 global spiral arms structure. While Hyodo et al. (2015) also performed N-body simulation with N = 3 × 104 − 5 × 104, they did not clearly show these two distinct structures, because they used ten times higher Mdisk (accordingly, ten times fewer mres) 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 mres = 3 spiral arms by the Lindblad torque and scatters or accretes nearby disk particles to open a gap with the half width ∼3rHill, where rHill is the satellite Hill’s radius defined by rHill = (Ms/3Mp)1/3as (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 (as ≃ aR), Eq. (4) at the disk outer edge (r ∼ aR − 3rHill) is

(5)

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

After the gap opening, the satellite migrates outward, and mres decreases from 3 to 2 (Fig. 1c). When the satellite orbit expands beyond the 2:1 resonance with the disk’s outer edge, mres becomes much smaller than 2 and the spiral arms disappears (Fig. 1d). Because the theoretically predicted effective viscosity is ∝Σ2 (Eq. (6)), the r-gradient of Σ is quickly flattened. The self-gravity wakes become fainter (Eq. (2)) through the loss of Mdisk 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 (as) and the total disk mass (Mdisk) obtained by our N-body simulation (RUN 1). The changes in as and Mdisk 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)

thumbnail 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 ≃ aR, , and represents the effect of finite physical size of the particles. Because we are concerned with the outer disk region, we adopted r ≃ 0.8aR. 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 Δas = as − aR. 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πΣνr2Ω) 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, das/dt from each run of our N-body simulations is compared with analytical estimations. In the analytical estimations, Ms and Mdisk obtained by the N-body simulations at each as are substituted to Eqs. (9) and (11) to calculate (das/dt)res and (das/dt)self. This figure shows that the results of N-body simulations fit (das/dt)self. In the vicinity of the disk’s edge, (das/dt)res dominates over (das/dt)self. The theoretical prediction for (das/(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.

thumbnail 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 Ms ∼ 10−3Mp, while the masses of the actual mid-sized moons are Ms ∼ (10−7–10−6)Mp. Hyodo et al. (2015) showed through N-body simulations that Ms/Mp ∼ 10(Mdisk(0)/Mp)2 for Ms ∼ 10−3Mp. Although Crida & Charnoz (2012) proposed Ms/Mp ∝ (Mdisk(0)/Mp)3 for the smaller value of Ms/Mp, generated clumps would quickly coagulate. Here we use the Hyodo et al. (2015) relation to discuss the cases of Ms ∼ (10−7–10−6)Mp. 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 ≲ (104–105) TKep, 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.

thumbnail 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 (MDione/MSaturn ≃ 1.94 × 10−6, MTethys/MSaturn ≃ 1.09 × 10−6), their tidal migration are convergent unless the planetary tidal parameter Qp 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 ≳ecrit 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 as 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 Ms/Mp ∼ 10−6. τa, crit ∼ 8.0 × 106TKep.

If we consider the Tethys-mass satellite (Ms ∼ 10−6Mp), the disk mass may be Mdisk ∼ 3 × 10−4Mp (Hyodo et al. 2015). The orbital migration timescale of the self-gravity wake is estimated to be τa, self ∼ 1.2 × 103(Mdisk/3 × 10−4Mp)−3TKep (Eq. (11)), which is shorter than ta,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 Mdisk, 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 ∼ 105) 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

  1. Charnoz, S., Crida, A., Castillo-Rogez, J. C., et al. 2011, Icarus, 216, 535 [Google Scholar]
  2. Crida, A., & Charnoz, S. 2012, Science, 338, 1196 [NASA ADS] [CrossRef] [Google Scholar]
  3. Daisaka, H., Tanaka, H., & Ida, S. 2001, Icarus, 154, 296 [Google Scholar]
  4. Dermott, S. F., Malhotra, R., & Murray, C. D. 1988, Icarus, 76, 295 [NASA ADS] [CrossRef] [Google Scholar]
  5. Duffell, P. C., & Chiang, E. 2015, ApJ, 812, 94 [Google Scholar]
  6. Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024 [Google Scholar]
  7. Goldreich, P., & Tremaine, S. 1982, ARA&A, 20, 249 [NASA ADS] [CrossRef] [Google Scholar]
  8. Hyodo, R., Ohtsuki, K., & Takeda, T. 2015, ApJ, 799, 40 [NASA ADS] [CrossRef] [Google Scholar]
  9. Ida, S. 2019, Science, 364, 1028 [Google Scholar]
  10. Iess, L., Militzer, B., Kaspi, Y., et al. 2019, Science, 364, aat2965 [NASA ADS] [CrossRef] [Google Scholar]
  11. Iwasawa, M., Tanikawa, A., Hosono, N., et al. 2016, PASJ, 68, 54 [Google Scholar]
  12. Lainey, V., Karatekin, Ö., Desmars, J., et al. 2012, ApJ, 752, 14 [Google Scholar]
  13. Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286 [NASA ADS] [CrossRef] [Google Scholar]
  14. Lainey, V., Casajus, L. G., Fuller, J., et al. 2020, Nat. Astron., in press, [arXiv:2006.06854] [Google Scholar]
  15. Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 [NASA ADS] [CrossRef] [Google Scholar]
  16. Malhotra, R. 1993, Icarus, 106, 264 [Google Scholar]
  17. Malhotra, R. 1996, AJ, 111, 504 [NASA ADS] [CrossRef] [Google Scholar]
  18. Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (Cambridge, UK: Cambridge University Press) [Google Scholar]
  19. Nakajima, A., Ida, S., Kimura, J., & Brasser, R. 2019, Icarus, 317, 570 [Google Scholar]
  20. Ogihara, M., & Kobayashi, H. 2013, ApJ, 775, 34 [NASA ADS] [CrossRef] [Google Scholar]
  21. Oshino, S., Funato, Y., & Makino, J. 2011, PASJ, 63, 881 [NASA ADS] [Google Scholar]
  22. Salmon, J., Charnoz, S., Crida, A., & Brahic, A. 2010, Icarus, 209, 771 [CrossRef] [Google Scholar]
  23. Salo, H. 1995, Icarus, 117, 287 [Google Scholar]
  24. Takeda, T., & Ida, S. 2001, ApJ, 560, 514 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1.

Parameter sets of our simulations.

All Figures

thumbnail Fig. 1.

Time evolution of the system of RUN1: (a) t = 0, (b) t = 29.4 TKep, (c) t = 1.47 × 103TKep, and (d) t = 2.94 × 104TKep, 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
thumbnail 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
thumbnail 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
thumbnail Fig. 4.

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

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.