# Dynamical Evolution of the Debris Disk after a Satellite Catastrophic Disruption around Saturn

Ryuki Hyodo, Sébastien Charnoz

► To cite this version:

Ryuki Hyodo, Sébastien Charnoz. Dynamical Evolution of the Debris Disk after a Satellite Catastrophic Disruption around Saturn. The Astronomical Journal, 2017, 154 (1), pp.34. 10.3847/1538-3881/aa74c9. insu-02914711

HAL Id: insu-02914711

https://insu.hal.science/insu-02914711v1

Submitted on 2 Sep 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Dynamical Evolution of the Debris Disk after a Satellite Catastrophic Disruption around Saturn

Ryuki Hyodo<fnref n="1" />,<fnref n="2" /> and Sébastien Charnoz<fnref n="2" />

Received 2017 March 22; revised 2017 May 13; accepted 2017 May 21; published 2017 June 30

## Abstract

The hypothesis of the recent origin of Saturn's rings and its midsized moons is actively debated. It was suggested that a proto-Rhea and a proto-Dione might have collided recently, giving birth to the modern system of midsized moons. It has also been suggested that the rapid viscous spreading of the debris may have implanted mass inside Saturn's Roche limit, giving birth to its modern ring system. However, this scenario has only been investigated in a very simplified way for the moment. This paper investigates it in detail to assess its plausibility by using *N*-body simulations and analytical arguments. When the debris disk is dominated by its largest remnant, *N*-body simulations show that the system quickly reaccretes into a single satellite without significant spreading. On the other hand, if the disk is composed of small particles, analytical arguments suggest that the disk experiences dynamical evolutions in three steps. The disk starts significantly excited after the impact and collisional damping dominates over the viscous spreading. After the system flattens, the system can become gravitationally unstable when particles are smaller than ~100 m. However, the particles grow faster than spreading. Then, the system becomes gravitationally stable again and accretion continues at a slower pace, but spreading is inhibited. Therefore, the debris is expected to reaccrete into several large bodies. In conclusion, our results show that such a scenario may not form today's ring system. In contrast, our results suggest that today's midsized moons are likely reaccreted from such a catastrophic event.

**Key words:** planets and satellites: dynamical evolution and stability – planets and satellites: individual (Tethys, Dione, Rhea, Titan) – planets and satellites: rings

## 1. Introduction

Origin, age, and dynamical evolution of icy Saturn's rings and satellites are still debated. Canup (2010) has proposed that Saturn's rings formed by tidal disruption of a Titan-sized body that migrates inward through the interaction with the circum-planetary gas disk about 4.5 Gyr ago. On the other hand, Hyodo et al. (2017) showed that tidal disruption of a passing Pluto-sized Kuiper Belt object can form ancient massive rings around not only Saturn but also other giant planets during the late heavy bombardment (LHB) about 3.8 Gyr ago. Then, the inner regular satellite systems around Saturn, Uranus, and Neptune are generally thought to be formed by spreading of such ancient massive rings (Charnoz et al. 2010; Crida & Charnoz 2012; Hyodo et al. 2015; Hyodo & Ohtsuki 2015).

The pure icy rings would continuously darken over the age of solar system due to micrometeorid bombardment (e.g., Cuzzi & Estrada 1998). Thus, the rings might have formed more recently than previously thought. However, note that they might be older if they are more massive (Elliott & Esposito 2011; Esposito et al. 2012).

Recently, Cuk et al. (2016) has investigated the past orbital evolutions of Saturn's midsized moons (Tethys, Dione and Rhea) and found that Tethys–Dione 3:2 orbital resonance is not likely to have occurred whereas the Dione–Rhea 5:3 resonance may have occurred. Then, they conclude that the midsized moons are not primordial and propose that the moons reaccreted from the debris disk that formed by a catastrophic collision between primordial Rhea-sized moons about 100 Myr ago (Cuk et al. 2016). They also propose that the debris disk may spread inward rapidly (due to fast gravitational instability) and feed the Roche limit to form today's rings. In addition they propose that outward spreading may form and push outward a

population of small moons (with a mass of $m = 4 \times 10^{20}$ kg) that would excite Titan's current eccentricity through the resonant interaction.

The aim of the present paper is to test this scenario by using direct simulations and detailed analytical arguments. In Section 2, we first use smoothed-particle hydrodynamics (SPH) simulations to investigate the outcome of the collision between two proto-Rhea-sized objects at impact velocity $3 \text{ km s}^{-1}$ (Cuk et al. 2016). In Section 3, using *N*-body simulations, we investigate the long-term evolution of the debris, starting from the impact simulation and assuming that debris is not collisionally disrupted. In Section 4, using analytical arguments, we estimate the fate of the disk of small particles as an extreme case of collisional evolution. In Section 5, we discuss the plausibility of this scenario to form today's rings and moons.

## 2. Catastrophic Collision between Rhea-sized Bodies

### 2.1. SPH Methods and Models

Using SPH simulations, we model the collision between Rhea-sized objects ($M_{\text{body}} = 10^{21}$ kg) in free space. The silicate mass fraction of Saturn's icy moons are diverse (Charnoz et al. 2011). Thus, we assume 60wt% silicate core for one object and 40wt% silicate core for the other with both covered by an icy mantel. Following Cuk et al. (2016) arguments, impact velocity is set to be about six times the mutual escape velocity, which is about $v_{\text{imp}} = 3 \text{ km s}^{-1}$. The impact angle is set to be either $\theta = 0^{\circ}, 20^{\circ}, 45^{\circ}, 60^{\circ}$, and $80^{\circ}$. The total mass of the two colliding objects is $M_{\text{tot}} = 2 \times 10^{21}$ kg and the total number of SPH particles is $N = 2 \times 10^{5}$. We simulated about 3.88 hr, which is much shorter than the orbital period at the distance of

<footnote n="1">¹Earth-Life Science Institute/Tokyo Institute of Technology, 2-12-1 Tokyo, Japan</footnote>

<footnote n="2">²Institut de Physique du Globe, 75005 Paris, France</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Largest two remnants after the collision as a function of different impact angles obtained from SPH simulations. Filled circles and squares represent the largest and the second largest fragments, respectively.</figcaption></figure>

Rhea (4.5 days). Our numerical code is the same as that used in Hyodo et al. (2016, 2017), which was developed in Genda et al. (2012).

## 2.2. Results of SPH Simulations

SPH simulations show that the collision is energetic enough to catastrophically destroy colliding objects (Figure 1) as suggested by Cuk et al. (2016). However, after the collision, in most cases, two large fragments remain as direct leftovers of the cores covered by water ice of the original two colliding objects. In the case of $\theta = 45^{\circ}$, the largest remnants consist of masses of $M = 7.8 \times 10^{20}$ kg and $M = 7.2 \times 10^{20}$ kg, which are both about 40% of the total mass of the two objects. Figure 2 shows the orbital elements of the debris after the impact in the case of $\theta = 45^{\circ}$, assuming the impact occurs at semimajor axis $a = 5 \times 10^5$ km (as in Cuk et al. (2016) and used as initial condition for N-body simulations; Section 3). Initial dispersion of the semimajor axes and eccentricities are about $3.5 \times 10^5$ km and 0.35, respectively, which are consistent with what we can derive from the first-order approximation as

$$ \Delta a_{\text{ini}} \sim 2\Delta v/\Omega \qquad (1) $$

$$ \Delta e_{\text{ini}} \sim \Delta v/(a\Omega) \qquad (2) $$

where $\Delta v \sim v_{\text{imp}}$ and $\Omega$ are the velocity dispersion and orbital frequency, respectively. In the next section, we investigate the longer-term evolution of the debris.

# 3. Dynamical Evolution of Debris with Large Fragments

## 3.1. N-body Methods and Models

Orbits of the debris are integrated by using a fourth-order Hermite method (Makino & Aarseth 1992). The collisions between particles are solved as a hard-sphere model with the normal and tangential coefficient of restitutions $\epsilon_n = 0.1$ and $\epsilon_t = 1$, respectively. However, following the argument of Kokubo et al. (2000) and Canup & Esposito (1995), we allow accretion only when the following two conditions are satisfied.

First, the Jacobi energy of two particles after the collision $E_J$ has to be negative as

$$ E_J = \frac{1}{2}v^2\epsilon_{\text{eff}}^2 - \frac{3}{2}x^2\Omega^2 + \frac{1}{2}z^2\Omega^2 - \frac{G(m_1 + m_2)}{r} + \frac{9}{2}r_{\text{H}}^2\Omega^2 < 0, \qquad (3) $$

where $x$, $y$, and $z$ are the relative positions, $r^2 = x^2 + y^2 + z^2$, $m_{1,2}$ are the masses of particles, and $\epsilon_{\text{eff}}$ is an effective coefficient of restitution written as

$$ \epsilon_{\text{eff}} = [(\epsilon_n^2 v_n^2 + \epsilon_t^2 v_t^2)/(v_n^2 + v_t^2)]^{1/2}, \qquad (4) $$

where $v_n$ and $v_t$ are the normal and tangential components of the relative velocity between particles. In addition, the sum of the radii of two particles should be smaller than the Hill radius as

$$ r_1 + r_2 = 3^{1/3} \frac{1 + \mu^{1/3}}{(1 + \mu)^{1/3}} \left( \frac{\rho_{\text{par}}}{\rho_{\text{pla}}} \right)^{-1/3} \frac{R_{\text{pla}}}{a} r_{\text{H}} \le r_{\text{H}}, \qquad (5) $$

where $\rho_{\text{pla}}$ and $\rho_{\text{par}}$ are the densities of planet and particles, respectively. $\mu$ is the mass ration $m_2/m_1$. $R_{\text{pla}}$ is the radius of the planet, and $r_H$ is the Hill radius defined as

$$ r_H = \left( \frac{m_1 + m_2}{3M_{\text{pla}}} \right)^{1/3} a. \qquad (6) $$

We use the tree-method for the gravity calculations and collisional detections (Rein & Liu 2012; Hyodo et al. 2015). The numerical code is the same as that used in Hyodo et al. (2015).

## 3.2. Initial Conditions

Positions and velocities of particles obtained from the SPH simulation with $\theta = 45^{\circ}$ (see Figure 2) are passed to N-body simulations, assuming the collision takes place in the equatorial plane of Saturn and the center of mass of the two colliding objects' orbits around Saturn with a semimajor axis $a = 5.0 \times 10^5$ km and the eccentricity $e = 0$. We also include Titan with the current semimajor axis $a_{\text{Titan}} = 1.2 \times 10^6$ km, eccentricity $e_{\text{Titan}} = 0.0288$, and inclination $i_{\text{Titan}} = 0.34^{\circ}$. Due to the computational power limitation, we randomly select 20,000 particles from 200,000 particles used in SPH simulations. We run five different simulations by changing the random choice of particles. Initially, each particle has the same mass of $m = M_{\text{tot}}/N$ ($m = 1 \times 10^{17}$ kg) and they are either silicate or icy particles. We assume silicate particles have density $\rho_{\text{sil}} = 3000$ kg m$^{-3}$ and icy particles have $\rho_{\text{icy}} = 900$ kg m$^{-3}$. During the calculation, we track the density change when two particles merge into a new particle. Just after the calculations start, numerous particles merge into single particles because they are initially the constituent particles of large remnants.

## 3.3. Results of N-body Simulations

Figure 3 shows the time evolution of the system. Just after the impact, most of the mass is contained in the two largest remnants (Figure 3, panel (a)). Since the two remnants have large eccentricities ($e \sim 0.2$), their orbits cross. Thus, after several periods, they collide and merge into a single large body
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Orbital elements of the debris after the impact in the case of $\theta = 45^{\circ}$ obtained from SPH simulations, assuming the collision takes place at $a = 5 \times 10^5$ km. The left and right panels show eccentricities and inclinations of particles against their semimajor axes. The blue dots represent water ice particles and red dots represent silicate particles. Cyan area in the left panel corresponds to the region inside the Roche limit, where $a(1 - e) < 135,000$ km.</figcaption></figure>

with a mass of $m \sim 1.5 \times 10^{21}$ kg with small eccentricity (Figure 3, panels (b) and (c)).

Masses of most field particles are $m_p = 10^{17}$ kg and their escape velocity is $v_{esc} \sim 20$ m s$^{-1}$. In order for accretion between such particles to take place, relative velocities should be smaller than their escape velocity. Thus, in order to accrete, the eccentricity of field particles should be smaller than $e_{cri} \sim 2.5 \times 10^{-3}$. Left panel of Figure 4 shows the time evolution of the root mean square (rms) of eccentricities $\langle e^2 \rangle^{1/2}$. Since the field particles have much larger eccentricities than $e_{cri}$, accretion between field particles is initially difficult. However, collisional damping is effective and the rms eccentricity decreases with time (Figure 4, left panel). Because the largest remnant is much larger than the field particles, field particles whose eccentricities are below $e \sim 0.06$ can accrete onto the largest remnant rather than between themselves. We confirm by N-body simulations that the remnant keeps growing by eating field particles and the number of particles in the system keeps decreasing (Figure 4, right panel).

At the end of our N-body simulations, we have less than 1000 particles without significant spreading of the system (Figure 3, panel (d)). At this time, the largest remnant (satellite) has accreted most of the field particles whose orbits cross that of the satellite and it has a mass of $m \sim 1.9 \times 10^{21}$ kg (~95% of the total system mass) with small eccentricity $e \sim 10^{-2}$. The size of the Hill sphere of this largest remnant is about 5000 km. The typical separation between two bodies is 10 Hill radius (Kokubo & Ida 1998). Thus, the remaining field particles would accrete onto the largest remnant and the system is expected to reaccrete into a single large object.

## 4. Dynamical Evolution of Debris of Small Particles

In the previous section, using N-body simulations, we investigated the long-term evolution of the debris disk within which initially the two large fragments are embedded as a result of catastrophic collision. However, we neglected the

effect of fragmentation and the debris particles initially have large eccentricities, thus collisional grinding may occur in the real system. Here, we analytically estimate the fate of the debris, initially consisting of same-sized small particles (radius $r_p$). The velocity dispersion $c$ of the system is controlled by the following equation.

$$ \frac{dc^2}{dt} = f_{\text{trans}} c^2 + f_{\text{col}} (r_p \Omega)^2 + f_{\text{vs}} \frac{v_{\text{esc}}^4}{c^2} - f_{\text{damp}} c^2, \quad (7) $$

where the first two terms are the contribution of viscous heating: the first term is due to velocity shear sampled by the random motion of particles (Goldreich & Tremaine 1978) and the second term is due to physical collisions (Araki & Tremaine 1986). The third term is due to gravitational scattering described by Chandrasekhar's relaxation time (Ida 1990; Michikoshi & Kokubo 2016) and the last term is due to collisional damping (Goldreich & Tremaine 1978). The coefficients are written as

$$ f_{\text{trans}} = c_1 \times \frac{9}{4} \frac{\tau}{1 + \tau^2} \Omega \quad (8) $$

$$ f_{\text{col}} = \frac{9}{4} \tau \Omega \quad (9) $$

$$ f_{\text{vs}} = \frac{\Omega \tau \ln \Lambda}{4} \quad (10) $$

$$ f_{\text{damp}} = c_2 \times \Omega \tau (1 - \epsilon^2) \quad (11) $$

where $\tau$ is the optical depth and is written with the assumption that all particles have the same radius $r_p$ as

$$ \tau = \frac{N_{\text{tot}} \pi r_p^2}{S} \sim 1.1 \times \left( \frac{\rho_p}{1200 \text{ kg m}^{-3}} \right)^{-1} \left( \frac{S}{1.1 \times 10^{18} \text{ m}^2} \right)^{-1} \times \left( \frac{M_{\text{tot}}}{2 \times 10^{21} \text{ kg}} \right) \left( \frac{r_p}{1.0 \text{ m}} \right)^{-1} \quad (12) $$
<figure><img src="image_5.png" /><figcaption>Figure 3. Time evolution of a debris disk on the a-e plane. The dots represent particles and their color and size represent their mass. The black filled big dot on the bottom right represents Titan. Cyan area on the left panel corresponds to the region inside the Roche limit, where a(1 - e) < 135,000 km. Two black lines in panel (d) represent the orbital elements that cross the orbits of the largest remnant at either pericenter or apocenter.</figcaption></figure>

where N<sub>tot</sub> is the total number of particles, ρ<sub>p</sub> is the particle density, and S is the surface area, respectively. Assuming ρ<sub>p</sub> = 1200 kg m<sup>-3</sup> and S = 2πaΔa = 2π × (5 × 10<sup>8</sup> m) × (3.5 × 10<sup>8</sup> m) ~ 1.1 × 10<sup>18</sup> m<sup>2</sup>, τ is in the range between τ = 10<sup>-5</sup> – 10<sup>3</sup>, depending on the size of particles between r<sub>p</sub> = 10<sup>-3</sup> – 10<sup>5</sup> m. ϵ is the coefficient of restitution and is in the range between 0–1 depending on the material properties and we use ϵ = 0.1 for our calculation. In Λ is in the range between 1–10, respectively. The coefficients c<sub>1</sub> and c<sub>2</sub> are of the order of unity and depend on τ (Goldreich & Tremaine 1978) and/or the spin state of particles (Morishima & Salo 2006). The dynamical evolution of the debris can be divided into three stages that we will discuss in detail in the following subsections. At each stage, we compare timescales of accretion, damping, and spreading.

## 4.1. Collisional Damping of the Initial Hot Debris

### 4.1.1. Collisional Damping Timescale

After the giant impact, the velocity dispersion of particles is much larger than their escape velocity and their shear velocity. Initially, the accretion is prohibited. Instead, collisional damping is effective and velocity dispersion gradually decreases. In the particle-in-a-box approximation, the collision timescale is written as

$$T_{\text{col}} = \frac{1}{n\sigma_{\text{col}}v_{\text{rel}}} \qquad (13)$$

where n is the number density of particles, σ<sub>col</sub> is the collisional cross-section, and v<sub>rel</sub> ~ c is the relative velocity. The cross-section is written as

$$\sigma_{\text{col}} = \pi r_{\text{p}}^2 (1 + v_{\text{esc}}^2 / v_{\text{rel}}^2). \qquad (14)$$

Considering the particles are distributed toroidally after the impact, the volume of this toroid can be expressed as V = (2πa) ⋅ (πa⟨e⟩a⟨I⟩) = 2π<sup>2</sup>a<sup>3</sup>⟨e⟩⟨I⟩, assuming radial and vertical widths are a⟨e⟩ and a⟨I⟩, respectively, where ⟨e⟩ and ⟨I⟩ are the mean eccentricity and inclination, respectively. Thus, the number density is written as

$$n = N/V = \frac{N}{2\pi^2 a^3 \langle e \rangle \langle I \rangle}, \qquad (15)$$

where we assume that ⟨e⟩ ~ √(5/3)⟨v<sub>rel</sub>/v<sub>K</sub>⟩ and ⟨I⟩ ~ √(1/3)⟨v<sub>rel</sub>/v<sub>K</sub>⟩, where v<sub>K</sub> is the Keplerian velocity (see also Jackson & Wyatt 2012). Figure 5 shows the collision timescale as a function of r<sub>p</sub> and velocity dispersion v<sub>rel</sub> ~ c. The timescale varies significantly depending on the size of the particle and relative velocity. We will compare this timescale to the viscous spreading timescale in the next subsection.

### 4.1.2. Spreading Timescale without Gravitational Instability

As the velocity dispersion decreases, the system may viscously spread. The timescale of viscous spreading can be written as T<sub>spr</sub> = Δa<sup>2</sup>/ν, where Δa is the diffusion width and ν is viscosity, respectively. The value of viscosity depends on
<figure><img src="image_6.png" /><figcaption>Figure 4. Evolution of rms eccentricity (left panel) and number of particles in the system (right panel) in the N-body simulation after the impact.</figcaption></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Figure 5. Collision timescale as a function of particle size and velocity dispersion (Equation (13)).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 6. Spreading timescale when $Q > 1$ as a function of particle size and velocity dispersion (Equation (20)).</figcaption></figure>

Toomre's Q parameter (Toomre 1964)

$$Q = \frac{c_r \kappa}{3.36 G \Sigma} \qquad (16)$$

where $c_r$ is the velocity dispersion in the radial direction and $\kappa$ is the epicyclic frequency, respectively. Initially, $Q$ is much larger than 1 and thus gravitationally stable. Therefore, the viscosity can be expressed as

$$\nu_{Q>1} = \nu_{\text{trans}} + \nu_{\text{col}}, \qquad (17)$$

where $\nu_{\text{trans}}$ is the translational viscosity (Goldreich & Tremaine 1978)

$$\nu_{\text{trans}} = \frac{c^2}{\Omega} \frac{\tau}{1 + \tau^2} \qquad (18)$$

and $\nu_{\text{col}}$ is the collisional viscosity (Araki & Tremaine 1986)

$$\nu_{\text{col}} = \Omega r_p^2 \tau, \qquad (19)$$

respectively. Then, the spreading timescale $\tau_{\text{spr},Q>1}$ can be written as

$$T_{\text{spr},Q>1} = \frac{\Delta a^2}{\nu_{Q>1}}. \qquad (20)$$

Figure 6 shows a spreading timescale when $Q > 1$ (Equation (20)) as a function of velocity dispersion and size of particle. Figure 7 shows the ratio of the collision timescale to the spreading timescale $T_{\text{col}}/T_{\text{spr},Q>1}$. We find that the collisional damping significantly dominates over the spreading in most of the parameter space considered here ($r_p = 10^{-3}-10^5$ m and $c = 1-10^4$ m s$^{-1}$). Thus, the initial hot debris disk is expected to flatten without significant spreading and accretion.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Ratio of the collision timescale to the spreading timescale $T_{col}/T_{spr,Q>1}$, when $Q > 1$. In the case ($T_{col}/T_{spr,Q>1}$) > 1, it is plotted with white and in the case ($T_{col}/T_{spr,Q>1}$) < 1 × 10⁻¹⁰, it is plotted with black.</figcaption></figure>

## 4.2. Accretion under Gravitational Instability

### 4.2.1. Spreading Timescale with Gravitational Instability

The ratio of the size of the Hill sphere to the sum of the particle radii is written as (see also Hyodo & Ohtsuki 2014)

$$ \tilde{r}_{\mathrm{H}} = \frac{R_{\mathrm{H}}}{2r_{\mathrm{p}}} = 0.82 \left( \frac{\rho}{900 \mathrm{kg} \mathrm{m}^{-3}} \right)^{1/3} \left( \frac{a}{100,000 \mathrm{km}} \right), \quad (21) $$

where $\rho$ is the density of the particle. Using $\rho = 1200 \mathrm{kg} \mathrm{m}^{-3}$ and $a = 500,000 \mathrm{km}$, we get $\tilde{r}_{\mathrm{H}} \sim 5$. As discussed above, the initial velocity dispersion decreases due to the collisional damping. Once the velocity dispersion becomes small enough, gravitational scattering becomes effective and increases the velocity dispersion. When $\tilde{r}_{\mathrm{H}} > 0.5$, the velocity dispersion at the steady state becomes comparable to the escape velocity of particles (Salo 1995; Ohtsuki 1999) as

$$ c \sim v_{\mathrm{esc}} = \sqrt{\frac{2Gm_{\mathrm{p}}}{r_{\mathrm{p}}}} = 0.07 \mathrm{m} \mathrm{s}^{-1} \left( \frac{\rho}{900 \mathrm{kg} \mathrm{m}^{-3}} \right)^{1/2} \left( \frac{r_{\mathrm{p}}}{100 \mathrm{m}} \right), \quad (22) $$

where $m_{\mathrm{p}}$ is the particle mass. In this second stage, the $Q$ parameter can become small. Figure 8 shows the value of Toomre's $Q$ parameter, assuming $c_r = v_{\mathrm{esc}}$ and $\kappa = \Omega$ with $a = 500,000 \mathrm{km}$. We find that $Q$ becomes smaller than 1 when particle radius $\lesssim 100 \mathrm{m}$. (Note that $r_{\mathrm{p}} = 100 \mathrm{m}$ corresponds to $\tau = 0.01$ in our work). Therefore, when $r_{\mathrm{p}} \lesssim 100 \mathrm{m}$, gravitational instability occurs ($Q < 1$). In this case, the gravitational viscosity dominates over that of collision, and the viscosity can be expressed as (Daisaka et al. 2001)

$$ v_{Q<1} = 2v_{\mathrm{grav}} = 52\tilde{r}_{\mathrm{H}}^5 \frac{G^2\Sigma^2}{\Omega^3}. \quad (23) $$

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Toomre's $Q$ parameter as a function of the size of the particle, assuming $c_r = v_{\mathrm{esc}}$ in Equation (16).</figcaption></figure>

Thus, the spreading timescale can be written as

$$ T_{\mathrm{spr},Q<1} = \frac{\Delta a^2}{v_{Q<1}} \sim 1 \times 10^4 \mathrm{year} \left( \frac{\Delta a}{3.5 \times 10^5 \mathrm{km}} \right)^2 \times \left( \frac{\Sigma}{2000 \mathrm{kg} \mathrm{m}^{-2}} \right)^{-2} \left( \frac{\Omega}{2 \times 10^{-5} \mathrm{s}} \right)^3 \left( \frac{\tilde{r}_{\mathrm{H}}}{5} \right)^{-5}. \quad (24) $$

Using $\tilde{r}_{\mathrm{H}} = 5$, $\Sigma = M_{\mathrm{tot}}/(2\pi a \Delta a_{\mathrm{ini}}) \sim 2000 \mathrm{kg} \mathrm{m}^{-3}$, $\Omega \sim 2 \times 10^{-5} \mathrm{s}$, and $\Delta a = 3.5 \times 10^8 \mathrm{m}$, we get $T_{\mathrm{spr},Q<1} \sim 1 \times 10^4 \mathrm{year}$. Compared to the case of $Q > 1$ (see Figure 6 and Equation (20)), the spreading timescale is significantly shorter for this small velocity dispersion (comparable to escape velocity). Thus, spreading may occur in this second stage. Next, we will compare this timescale to the accretion timescale in the next subsection.

### 4.2.2. Accretion Timescale

Since the velocity dispersion is now small, particles can accrete and grow. As discussed above, once the particle becomes larger than a 100 m sized body, the system becomes gravitationally stable ($Q > 1$ and see Figure 8). When the accretion timescale grows up to the size of $R$ and mass $M$, it is written as $T_{\mathrm{grow}} = M/\dot{M}$ and the growth rate $\dot{M}$ can be expressed as mass that is swept up per unit time as

$$ \dot{M} = (\Sigma/H) v_{\mathrm{rel}} \pi R^2 (1 + F_{\mathrm{grav}}), \quad (25) $$

where $H$ is the scale height and written as $H = v_{\mathrm{rel}}/\Omega$ and $F_{\mathrm{grav}} = v_{\mathrm{esc}}^2/v_{\mathrm{rel}}^2$ is the gravitational focusing factor, respectively. Thus, the growth timescale becomes

$$ T_{\mathrm{grow}} \sim \frac{\rho R}{\Omega \Sigma (1 + F_{\mathrm{grav}})} = 0.07 \mathrm{year} \frac{1}{1 + F_{\mathrm{grav}}} \left( \frac{\rho}{900 \mathrm{kg} \mathrm{m}^{-3}} \right) \times \left( \frac{R}{100 \mathrm{m}} \right) \left( \frac{\Omega}{2 \times 10^{-5} \mathrm{s}} \right)^{-1} \left( \frac{\Sigma}{2000 \mathrm{kg} \mathrm{m}^{-2}} \right)^{-1}. \quad (26) $$
The timescale to form a 100 m sized body is independent of the initial size of the particle, as seen in Equation (26), and, considering $\rho = 1200 \text{ kg m}^{-3}$, $R = 100 \text{ m}$, $\Sigma = 2000 \text{ kg m}^{-2}$, $\Omega = 2 \times 10^{-5} \text{ s}$, and $F_{\text{grav}} = 1$, we get $T_{\text{grow}} \sim 0.05$ year, which is much shorter than what we obtained for spreading with $Q < 1$ (Equation (24)). Therefore, accretion takes place quickly without significant spreading even under the gravitational instability and forms particles larger than 100 m. Thus, the system again becomes gravitationally stable ($Q > 1$).

### 4.3. Accretion under Gravitational Stability

As discussed above, once the typical size of the particle becomes larger than 100 m at $a = 500,000 \text{ km}$, the system becomes gravitationally stable. Thus, the spreading timescale is regulated by $T_{\text{spr},Q>1}$ (Equation (20)), which is much longer than the accretion timescale of $T_{\text{grow}} \sim 500$ years even for 1000 km body using Equation (26).

Cuk et al. (2016) assume that the system is always gravitationally instable ($Q < 1$) and estimate the spreading timescale by using Equation (24) as about 2000 years, which is comparable to the timescale to form a 1000 km sized object (Equation (26)). Then, they proposed that the debris may spread all the way inside the Roche limit and form Saturn’s rings (Cuk et al. 2016). However, as we have shown above, the system is expected to accrete into several large objects without significant spreading. This is also confirmed by N-body simulations (Section 3) in the case where we start with large particles ($Q > 1$).

## 5. Conclusion and Discussion

Several scenarios exist for the origin of Saturn’s rings. Rings may form during the gas accretion phase (~4.5 Gyr ago) by tidal disruption of a gas-driven inward-migrating primordial satellite (Canup 2010) or it may have formed during LHB (~3.8 Gyr ago) by tidal disruption of passing large KBOs (Hyodo et al. 2017). In contrast, rings could be much younger than the solar system (Cuzzi & Estrada 1998). Recently, Cuk et al. (2016) proposed that Saturn’s moon system has experienced a catastrophic impact between Rhea-sized objects about 100 Myr ago around its current location and that the disk of debris may spread all the way inward to form rings. They also proposed that the current eccentricity of Titan could be induced by the orbital resonance with small moons that formed at the edge of the disk and migrate outward due to the interaction with the spreading disk.

In this paper, using both direct numerical simulations and analytical arguments, we investigated the hypothesis that is proposed in Cuk et al. (2016). First, we performed SPH simulations of the giant impact between Rhea-sized objects with an impact velocity of 3 km s⁻¹. We found that the outcome of the collision, if catastrophic (for impact angle 45°), in general, forms only two large remnants containing about 40% of the initial total moons’ mass. These fragments are embedded in a debris disk (Section 2). Then, we performed N-body simulations using the data obtained from SPH simulations to investigate the longer-term evolution of the debris disk (Section 3). N-body simulations suggest that the system quickly reaccretes into a single object without significant spreading of the debris.

However, in the N-body simulations, the effect of fragmentation is not included. After the giant impact, the debris

particles have large eccentricities and thus successive collisional grinding may occur. In addition, the size of fragments depends on the impact angle even though the impact velocity is the same (see Figure 1). Thus, using analytical arguments, we investigate the fate of the debris in the case in which they consist of only small particles (Section 4). We find that the system follows three different stages of dynamical evolution. Just after the impact, the system is significantly excited. At this time, Toomre’s Q parameter is larger than 1 and thus the viscosity of the debris is written as $\nu_{Q>1} = \nu_{\text{trans}} + \nu_{\text{col}}$ (Equation (17)). At this first stage, collision damping dominates over viscous spreading. Therefore, the system flattens until the velocity dispersion becomes comparable to the particle’s escape velocity (Section 4.1). Second, when the velocity dispersion becomes comparable to the escape velocity, the Q parameter can become smaller than 1 as long as the radius of particles is smaller than 100 m. Under this condition, the viscosity is regulated by gravitational interaction as $\nu_{Q<1} = 2\nu_{\text{grav}}$ (Equation (23)). Then, we calculated the accretion timescale up to a 100 m sized body and we found that the accretion timescale is much shorter than that of the spreading timescale. Therefore, at this second stage, the accretion dominates over the spreading (Section 4.2). After particles grow to sizes larger than 100 m, the system becomes $Q > 1$ again. Thus, the viscous spreading is regulated by $\nu_{Q>1}$. Comparing the timescale of viscous spreading to the accretion timescale to a 1000 km sized body, the accretion timescale is again much shorter than the spreading timescale as long as the velocity dispersion is comparable to or smaller than the escape velocity of particles. Thus, at this third stage, the accretion further takes place without significant spreading of the system (Section 4.3).

We find that the impact between the two moons is indeed catastrophic as suggested by Cuk et al. (2016). However, we do not find significant spreading, but rather rapid reaccretion of the system. Difference from Cuk et al. (2016) comes from the viscosity formula that is used. Cuk et al. (2016) assumes that the system is always gravitationally unstable ($Q < 1$) and applied the formula $\nu_{Q<1}$ to estimate the spreading timescale to compare the accretion timescale up to a 1000 km body. However, as we have shown above, the system is mostly gravitationally stable ($Q > 1$) and $\nu_{Q>1}$ should be considered.

In conclusion, this study shows that the debris is expected to reaccrete very quickly to form a new-Rhea or/and new-Dione and that spreading is very inefficient after the impact and before complete reaccretion. Therefore, as discussed above, the disk hardly spreads to form Saturn’s rings. Thus, the origin of Titan’s current eccentricity by disk-driven migration of small moons into orbital resonance with Titan as suggested by Cuk et al. (2016) is also less likely to occur.

We thank H. Genda for kindly providing us with an SPH code. R.H. thanks Shugo Michikoshi for discussion. We also thank L.W. Esposito for useful comments on the manuscript. This work was supported by JSPS Grants-in-Aid for JSPS Fellows (17J01269). Part of the numerical simulations were performed using the GRAPE system at the Center for Computational Astrophysics of the National Astronomical Observatory of Japan. Also, numerical computations were partly performed on the S-CAPAD platform, IPGP, France. We acknowledge the financial support of the UnivEarthS Labex
program at Sorbonne Paris Cité (ANR-10-LABX-0023 and ANR-11-IDEX-0005-02).

References

Araki, S., & Tremaine, S. 1986, Icar, 65, 83

Canup, R. 2010, Natur, 468, 943

Canup, R. M., & Esposito, L. W. 1995, Icar, 113, 331

Charnoz, S., Salmon, J., & Crida, A. 2010, Natur, 465, 752

Charnoz, S., Crida, A., Castillo-Rogez, J. C., et al. 2011, Icar, 216, 535

Crida, A., & Charnoz, S. 2012, Sci, 338, 1196

Cuk, M., Dones, L., & Nesvorny, D. 2016, ApJ, 820, 97

Cuzzi, J. N., & Estrada, P. R. 1998, Icar, 132, 1

Daisaka, H., Tanaka, H., & Ida, S. 2001, Icar, 154, 296

Elliott, J. P., & Esposito, L. W. 2011, Icar, 212, 268

Esposito, L. W., Albers, N., Meinke, B. K., et al. 2012, Icar, 217, 103

Genda, H., Kokubo, E., & Ida, S. 2012, ApJ, 744, 137

Goldreich, P., & Tremaine, S. 1978, Icar, 34, 227

Hyodo, R., Charnoz, S., Genda, H., & Ohtsuki, K. 2016, ApJL, 828, L8

Hyodo, R., Charnoz, S., Ohtsuki, K., & Genda, H. 2017, Icar, 282, 195

Hyodo, R., & Ohtsuki, K. 2014, ApJ, 787, 56

Hyodo, R., & Ohtsuki, K. 2015, NatGe, 8, 686

Hyodo, R., Ohtsuki, K., & Takeda, T. 2015, ApJ, 799, 40

Ida, S. 1990, Icar, 88, 129

Jackson, A. P., & Wyatt, M. C. 2012, MNRAS, 425, 657

Kokubo, E., & Ida, S. 1998, Icar, 131, 171

Kokubo, E., Ida, S., & Makino, J. 2000, Icar, 148, 419

Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141

Michikoshi, S., & Kokubo, E. 2016, ApJL, 825, L28

Morishima, R., & Salo, H. 2006, Icar, 181, 272

Ohtsuki, K. 1999, Icar, 137, 152

Rein, H., & Liu, S.-F. 2012, A&A, 537, A128

Salo, H. 1995, Icar, 117, 287

Toomre, A. 1964, ApJ, 139, 1217