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/00046361/202038743  
Published online  14 August 2020 
Letter to the Editor
Orbital evolution of Saturn’s satellites due to the interaction between the moons and the massive rings
^{1}
Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Tokyo 1528551, Japan
email: nakajima.a.ah@m.titech.ac.jp
^{2}
EarthLife Science Institute, Tokyo Institute of Technology, Tokyo, Japan
email: ida@elsi.ac.jp
^{3}
Department of Earth and Planetary Sciences, University of Tokyo, Tokyo, Japan
^{4}
ISAS, JAXA, Kanagawa, Japan
email: y.ishigaki@stp.isas.jaxa.jp
Received:
25
June
2020
Accepted:
20
July
2020
Context. Saturn’s midsized moons (satellites) have a puzzling orbital configuration with trapping in meanmotion resonances with everyother pairs (MimasTethys 4:2 and EnceladusDione 2:1). To reproduce their current orbital configuration on the basis of a recent model of satellite formation from a hypothetical ancient massive ring, adjacent pairs must pass firstorder meanmotion resonances without being trapped.
Aims. The trapping could be avoided by fast orbital migration and/or excitation of the satellite’s eccentricity caused by gravitational interactions between the satellites and the rings (the disk), which are still unknown. In our research we investigate the satellite orbital evolution due to interactions with the disk through full Nbody simulations.
Methods. We performed global highresolution Nbody simulations of a selfgravitating particle disk interacting with a single satellite. We used N ∼ 10^{5} particles for the disk. Gravitational forces of all the particles and their inelastic collisions are taken into account.
Results. Dense shortwavelength wake structure is created by the disk selfgravity and a few global spiral arms are induced by the satellite. The selfgravity wakes regulate the orbital evolution of the satellite, which has been considered as a disk spreading mechanism, but not as a driver for the orbital evolution.
Conclusions. The selfgravity wake torque to the satellite is so effective that the satellite migration is much faster than was predicted with the spiral arm torque. It provides a possible model to avoid the resonance capture of adjacent satellite pairs and establish the current orbital configuration of Saturn’s midsized satellites.
Key words: planets and satellites: dynamical evolution and stability / planets and satellites: rings / comets: individual: Saturn
© ESO 2020
1. Introduction
The orbital configuration of Saturn’s midsized moons, Mimas, Enceladus, Tethys, Dione, and Rhea (from inner to outer orbits), is puzzling: they are trapped in meanmotion resonances for everyother pairs (MimasTethys 4:2 and EnceladusDione 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 firstorder meanmotion 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^{19} ≃ 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 satellitedisk (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 satellitedisk 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 planetdisk interactions usually assumed protoplanetary gas disks that are stable against selfgravitational 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 selfgravitating particle disk (rings) via highresolution Nbody simulations.
Hyodo et al. (2015) performed highresolution Nbody simulations (N = 3 × 10^{4} − 5 × 10^{4}) 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 highresolution (N ∼ 10^{5}) Nbody 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 Nbody simulation code called GPLUM (Ishigaki, in prep.) that adopts the particleparticle particletree scheme (P^{3}T; Oshino et al. 2011) for planetary formation. The P^{3}T scheme uses the fourthorder Hermite integrator to calculate gravitational interactions between particles within a cutoff radius and the BarnesHut tree method for gravity from particles beyond the cutoff (Iwasawa et al. 2016), which guarantees higherorder 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 (freeslip 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
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 a_{R} ≃ 2.26 R_{p}. The initial surface density of the particles follows Σ(r)∝r^{−3}, and their total mass is ∼(10^{−3}–10^{−2}) M_{p}. We use 8 × 10^{4} − 1.2 × 10^{5} particles with equal mass of M ∼ (10^{−8}–10^{−7}) M_{p} for the disk. Initially, the particles have circular orbits with low inclinations, following a normal distribution of ⟨e^{2}⟩^{1/2} = 2⟨i^{2}⟩^{1/2} ∼ R/r ∼ (2 − 4)×10^{−3}, where R is the particle physical radius. Because of the inelastic collisions and selfgravity of the disk particles, they are quickly relaxed to quasiequilibrium 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^{−3}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 satellitedisk interactions at the radius inside the 2:1 resonance with the disk’s outer edge. The tidal edamping and aexpansion timescales are (Q_{s}/k_{2s})(M_{s}/M_{p})(a_{s}/R_{m})^{5}Ω^{−1} 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^{5}, M_{p}/M_{s} ∼ 10^{6}, and the satellite orbital radius a_{s} ∼ a_{R} ≃ 2.26R_{p}, years. For the planet tidal parameter Q_{p}/k_{2p} ∼ 10^{3}–10^{5}, –10^{9} years. As we show in Sect. 3, the aexpansion timescale due to the satellitedisk interactions is τ_{a, disk} ∼ (π^{2}/51)(M_{p}/M_{disk})^{3}(M_{s}/M_{p})(a_{R}/a_{s})^{1/2}Ω^{−1} ∼ 7 × 10^{3}Ω^{−1} ∼ 2 years for a realistic case with M_{disk}/M_{p} ∼ 3 × 10^{−4} and M_{s}/M_{p} ∼ 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 (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 midsize moons. We derive semianalytical formulas from the results of Nbody simulations to clarify intrinsic physics in this system. Applying the derived mass scaling law for the realistic masses of Saturn’s midsize moons, we discuss the possibility to avoid the resonance trapping.
Parameter sets of our simulations.
3. Simulation results
3.1. Ring structures
Figure 1c is a snapshot at t = 1.47 × 10^{3}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 selfgravity (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 selfgravity wakes are estimated as (Takeda & Ida 2001)
Fig. 1. Time evolution of the system of RUN1: (a) t = 0, (b) t = 29.4 T_{Kep}, (c) t = 1.47 × 10^{3} T_{Kep}, and (d) t = 2.94 × 10^{4} T_{Kep}, where T_{Kep} is the Keplerian period at r = a_{R}. 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 = a_{R}). 
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
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 Nbody simulation simultaneously show the shortwavelength wakes due to the selfgravity, 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 Nbody simulation with N = 3 × 10^{4} − 5 × 10^{4}, 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 selfgravity 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/3}a_{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
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 ∝Σ^{2} (Eq. (6)), the rgradient of Σ is quickly flattened. The selfgravity 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 Nbody simulation (RUN 1). The changes in a_{s} and M_{disk} become slower with t, which are theoretically explained as follows. Through Nbody simulations, Daisaka et al. (2001) found that the effective viscosity of a selfgravitating disk is given by
Fig. 2. Time evolution of the satellite’s semimajor axis (a_{s}; upper panel) and the disk mass (M_{disk}; lower panel). The units of the semimajor axis and the time are the Roche limit radius a_{R} and the Keplerian period at a_{R}, 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
Integrating this equation, we predict the explicit time evolution of the disk mass,
which reproduces the Nbody simulation result (lower panel of Fig. 2).
Using Eq. (8), we show that the migration is regulated by the selfgravity wakes, but not by the Lindblad resonance torque (the spiral arms) induced by the satellite. The migration rate due to the (onesided) Lindblad resonance is (e.g., Lin & Papaloizou 1986; Crida & Charnoz 2012)
where Δa_{s} = a_{s} − a_{R}. The migration rate by the selfgravity 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^{2}Ω) is transferred from the disk outer edge to the satellite’s orbit. In this case,
Substituting Eq. (6) and into Eq. (10), we obtain
In Fig. 3, da_{s}/dt from each run of our Nbody simulations is compared with analytical estimations. In the analytical estimations, M_{s} and M_{disk} obtained by the Nbody 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 Nbody 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 nonselfgravitating disk with modest viscosity. The spiral arms may be weakened by the relatively strong diffusion due to selfgravity wakes. Because the selfgravity 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, da_{s}/dt, as a function of a_{s} in RUN1 to RUN5. The red curves are the results of individual Nbody simulations. The blue and green curves are the theoretical predictions of (da_{s}/dt)_{res} and (da_{s}/dt)_{self} given respectively by Eqs. (9) and (11) using M_{s} and M_{disk} obtained by the Nbody simulations at individual a_{s} in each run. 
In these runs we adopted M_{s} ∼ 10^{−3}M_{p}, while the masses of the actual midsized moons are M_{s} ∼ (10^{−7}–10^{−6})M_{p}. Hyodo et al. (2015) showed through Nbody simulations that M_{s}/M_{p} ∼ 10(M_{disk}(0)/M_{p})^{2} for M_{s} ∼ 10^{−3}M_{p}. Although Crida & Charnoz (2012) proposed M_{s}/M_{p} ∝ (M_{disk}(0)/M_{p})^{3} 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^{−7}–10^{−6})M_{p}. As is shown in Sect. 4, the high migration rate induced by the selfgravity wake torque, which has been overlooked in the past studies, plays an important role in avoidance of the meanmotion 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^{4}–10^{5}) T_{Kep}, when the Lindblad torque may not be negligible compared to the selfgravity 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 meanresonance 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 meanmotion resonance capture depends on the eccentricity and migration rate of the satellites (Dermott et al. 1988; Malhotra 1993). Here, we consider the avoidance of TethysDione 3:2 resonance tapping as an example. Dione’s orbit is currently located beyond that of Tethys with the separation within 3:2 meanmotion 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^{−6}, M_{Tethys}/M_{Saturn} ≃ 1.09 × 10^{−6}), 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 meanmotion 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)
Nakajima et al. (2019) pointed out the possibility to avoid TethysDione 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 satellitedisk interactions. According to Ogihara & Kobayashi (2013), resonance trapping is avoided for TethysDione 3:2 resonance (j = 2) if
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^{−6}. τ_{a, crit} ∼ 8.0 × 10^{6}T_{Kep}.
If we consider the Tethysmass satellite (M_{s} ∼ 10^{−6}M_{p}), the disk mass may be M_{disk} ∼ 3 × 10^{−4}M_{p} (Hyodo et al. 2015). The orbital migration timescale of the selfgravity wake is estimated to be τ_{a, self} ∼ 1.2 × 10^{3}(M_{disk}/3 × 10^{−4}M_{p})^{−3}T_{Kep} (Eq. (11)), which is shorter than t_{a,crit} by more than three orders of magnitude. Although the migration by the selfgravity wake torque could be weaker along with the decrease in M_{disk}, the fast migration potentially prevents the resonance capture of the DioneTethys pair.
5. Conclusions
In order to investigate the gravitational interactions between Saturn’s midsized moons and hypothetical ancient massive rings and the associated orbital evolution of the moons, we performed global highresolution Nbody simulations (N ∼ 10^{5}) of a selfgravitating 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 shortwavelength wake structure and m = 2 or 3 global spiral arms simultaneously develop in the disk and are produced by the disk selfgravity 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 selfgravity wakes. The literature assumes that the Lindblad torque regulates the migrations because they considered the selfgravity 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 semianalytical formulas for the satellite’s migrating rate. While the simulations used a much more massive satellite than Saturn’s midsized 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 midsize 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 GrantsinAid for Scientific Research (# JP19J12542) and MEXT “Exploratory Challenge on PostK computer” hp190143. Numerical computations were in part carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan.
References
 Charnoz, S., Crida, A., CastilloRogez, J. C., et al. 2011, Icarus, 216, 535 [Google Scholar]
 Crida, A., & Charnoz, S. 2012, Science, 338, 1196 [NASA ADS] [CrossRef] [Google Scholar]
 Daisaka, H., Tanaka, H., & Ida, S. 2001, Icarus, 154, 296 [NASA ADS] [CrossRef] [Google Scholar]
 Dermott, S. F., Malhotra, R., & Murray, C. D. 1988, Icarus, 76, 295 [NASA ADS] [CrossRef] [Google Scholar]
 Duffell, P. C., & Chiang, E. 2015, ApJ, 812, 94 [Google Scholar]
 Goldreich, P., & Sari, R. 2003, ApJ, 585, 1024 [Google Scholar]
 Goldreich, P., & Tremaine, S. 1982, ARA&A, 20, 249 [NASA ADS] [CrossRef] [Google Scholar]
 Hyodo, R., Ohtsuki, K., & Takeda, T. 2015, ApJ, 799, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Ida, S. 2019, Science, 364, 1028 [NASA ADS] [CrossRef] [Google Scholar]
 Iess, L., Militzer, B., Kaspi, Y., et al. 2019, Science, 364, aat2965 [NASA ADS] [CrossRef] [Google Scholar]
 Iwasawa, M., Tanikawa, A., Hosono, N., et al. 2016, PASJ, 68, 54 [NASA ADS] [CrossRef] [Google Scholar]
 Lainey, V., Karatekin, Ö., Desmars, J., et al. 2012, ApJ, 752, 14 [Google Scholar]
 Lainey, V., Jacobson, R. A., Tajeddine, R., et al. 2017, Icarus, 281, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Lainey, V., Casajus, L. G., Fuller, J., et al. 2020, Nat. Astron., in press, [arXiv:2006.06854] [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 [NASA ADS] [CrossRef] [Google Scholar]
 Malhotra, R. 1993, Icarus, 106, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Malhotra, R. 1996, AJ, 111, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Nakajima, A., Ida, S., Kimura, J., & Brasser, R. 2019, Icarus, 317, 570 [CrossRef] [Google Scholar]
 Ogihara, M., & Kobayashi, H. 2013, ApJ, 775, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Oshino, S., Funato, Y., & Makino, J. 2011, PASJ, 63, 881 [NASA ADS] [Google Scholar]
 Salmon, J., Charnoz, S., Crida, A., & Brahic, A. 2010, Icarus, 209, 771 [CrossRef] [Google Scholar]
 Salo, H. 1995, Icarus, 117, 287 [NASA ADS] [CrossRef] [Google Scholar]
 Takeda, T., & Ida, S. 2001, ApJ, 560, 514 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1. Time evolution of the system of RUN1: (a) t = 0, (b) t = 29.4 T_{Kep}, (c) t = 1.47 × 10^{3} T_{Kep}, and (d) t = 2.94 × 10^{4} T_{Kep}, where T_{Kep} is the Keplerian period at r = a_{R}. 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 = a_{R}). 

In the text 
Fig. 2. Time evolution of the satellite’s semimajor axis (a_{s}; upper panel) and the disk mass (M_{disk}; lower panel). The units of the semimajor axis and the time are the Roche limit radius a_{R} and the Keplerian period at a_{R}, 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, da_{s}/dt, as a function of a_{s} in RUN1 to RUN5. The red curves are the results of individual Nbody simulations. The blue and green curves are the theoretical predictions of (da_{s}/dt)_{res} and (da_{s}/dt)_{self} given respectively by Eqs. (9) and (11) using M_{s} and M_{disk} obtained by the Nbody simulations at individual a_{s} 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 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.