Accretion of Saturn's Inner Mid-sized Moons from a Massive Primordial Ice Ring

The core numerical model used here is based on one developed to study the accretion of the Earth's Moon from a protolunar disk (Salmon & Canup 2012, 2014); additional details are contained in Salmon & Canup (2012), and the appendices therein. The code couples an analytical model of a viscous interior ring to the N-body code SyMBA (Duncan et al. 1998), which is used to simulate the accretion of moons exterior to the ring. The inner ring extends from the planet's surface at radius R_p to an outer edge r_out, which is initially set equal to the Roche limit, ${a}_{R}=1.524{({M}_{P}/\rho )}^{1/3}$, where M_P is the planet's mass and ρ is the density of ring material. We consider ice ring particles with $\rho =0.9\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ so that ${a}_{R}=2.24{R}_{\saturn}$. The ring's surface density, σ, is assumed to be constant across the radial extent of the ring, with

$$\begin{eqnarray}&&\sigma =\displaystyle \frac{{M}_{r}}{\pi ({r}_{\mathrm{out}}^{2}-{R}_{p}^{2})},\end{eqnarray} \tag{ 1 }$$

where M_r is the ring's total mass. We emphasize the distinction between Saturn's early radius (R_p) and its current mean radius (R_♄), where in our simulations R_p is larger than R_♄ because we consider a primordial ring and a young Saturn. An initial ice ring with ${M}_{r}={10}^{22}\,\mathrm{kg}$, ${r}_{\mathrm{out}}={a}_{R}$, and ${R}_{p}=1.4{R}_{\saturn}$ has $\sigma \approx 3\times {10}^{4}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. With time, M_r, σ, and r_out vary due to the ring's viscosity and interactions with outer moons.

The ring spreads with a viscosity ν that includes the effects of self-gravity (Ward & Cameron 1978; Salo 1995; Daisaka & Ida 1999; Daisaka et al. 2001),

$$\begin{eqnarray}&&\nu \approx \displaystyle \frac{{\pi }^{2}{G}^{2}{\sigma }^{2}}{{{\rm{\Omega }}}^{3}},\end{eqnarray} \tag{ 2 }$$

where G = 6.67 × 10⁻¹¹ m³ kg⁻¹ s⁻² is the gravitational constant and ${\rm{\Omega }}=\sqrt{({{GM}}_{P}/{r}^{3})}$ is the orbital frequency at distance r from the planet's center. In our calculation of Ω for the viscosity, we set $r={r}_{\mathrm{out}}$. Near the ring's inner edge, the viscosity would be lower for a fixed surface density because r is smaller. More detailed models of the rings' viscous evolution predict formation of a inner density peak (Salmon et al. 2010), which would tend to increase the viscosity through a higher value of σ. Our model does not resolve the radial structure of the rings, and applies the same viscosity across the entirety of the ring. Viscous spreading causes the ring to lose mass through its inner edge as mass is accreted by the planet and causes the outer edge of the ring to expand (see Appendix A in Salmon & Canup 2012 for details).

The N-body portion of the code tracks the orbital and collisional evolution of discrete objects beyond the Roche limit. Each outer object interacts with the inner ring at its strongest Lindblad resonances, resulting in a positive torque on the object and a negative torque on the ring, which causes r_out to contract. The total torque T_res exerted by the rings on an exterior satellite per unit satellite mass is found by summing the torques due to all the 0th order resonances that fall in the disk (Salmon & Canup 2012),

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{\mathrm{res}}}{m}=\left(\displaystyle \frac{{\pi }^{2}}{3}\mu G\sigma a\right)C(p),\end{eqnarray} \tag{ 3 }$$

where m is the satellite's mass, $\mu =m/{M}_{\saturn}$, $C(p)\,={\sum }_{p=2}^{{p}_{* }}2.55{p}^{2}(1-1/p)$, and p_* is the highest p for which resonance (p:p − 1) falls in the disk. Once an object is far enough from the ring that its strongest resonances no longer fall within the ring, which occurs for orbital radii $\geqslant 1.6{r}_{\mathrm{out}}$, it no longer interacts directly with the ring. The net change in the position of the ring's outer edge at each time step is found by considering the combined effect of resonant torques due to all of the exterior objects and the ring's viscosity. If the former dominates, the ring edge contracts, while if the latter dominates, r_out expands.

Ring material that spreads beyond the Roche limit can clump into tidally stable fragments due to local gravitational instabilities, which then mutually collide and accrete into still larger objects. The mass m_f of a fragment formed via local instability is

$$\begin{eqnarray}&&{m}_{f}\approx \displaystyle \frac{16{\pi }^{4}{\xi }^{2}{\sigma }^{3}{r}_{\mathrm{out}}^{6}}{{M}_{p}^{2}}\\ &&\,\approx \,5.8\times {10}^{13}{\rm{g}}{\left(\displaystyle \frac{{r}_{\mathrm{out}}}{{a}_{R}}\right)}^{6}\times {\left(\displaystyle \frac{\sigma }{3\times {10}^{4}{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{3},\end{eqnarray} \tag{ 4 }$$

where ξ is a factor of order but less than unity (Goldreich & Ward 1973). This mass is of order 10⁻¹³ times Saturn's mass, which is too small to be feasibly treated in our N-body simulations. Small fragments would likely collide and merge rapidly into bigger objects (Charnoz et al. 2010). In our simulations, we set the mass of objects spawned at the Roche limit to m_f = 10⁻⁸ M_♄, which is about one-tenth Mimas' mass, so that the growth of Mimas-sized moons is (marginally) resolved. When ${r}_{\mathrm{out}}\gt {a}_{R}$, we remove a mass m_f from the ring, and add a new discrete object having this mass to the N-body code at $r\approx {r}_{\mathrm{out}}$. The position of the ring's outer edge is then decreased to conserve angular momentum, such that ${L}_{d}+{L}_{f}={L}_{d,0}$, where ${L}_{d,0}$ and L_d are the angular momentum of the ring before and after the formation of the fragment, respectively, and L_f is the orbital angular momentum of the newly formed object.

We use tidal accretion criteria (Ohtsuki 1993; Canup & Esposito 1995) to determine if collisions between objects in the N-body code will result in a merger or intact rebound, assuming completely inelastic collisions. The outcome of a collision then depends on the impact energy, the mass ratio of the colliding bodies, and the orbital distance of the impact with respect to the Roche limit. We use a "total accretion" criterion, in which we assume that collisions occur in the radial direction along the widest axis of the Hill sphere of the colliding bodies, which is the most favorable case for accretion.

It is possible that interactions among orbiting bodies can cause an object to be scattered onto an orbit whose pericenter is close to the planet. In the limiting case of an inviscid fluid object on a parabolic orbit, tidal disruption will occur in a single pass once its pericenter r_p satisfies (Sridhar & Tremaine 1992)

$$\begin{eqnarray}&&{r}_{p}\lt 1.05{\left(\displaystyle \frac{{M}_{p}}{\rho }\right)}^{1/3}\approx 0.7{a}_{R}.\end{eqnarray} \tag{ 5 }$$

When an object satisfies this criterion, we remove it from the N-body code and add its mass and angular momentum to the ring. In practice, such events occur rarely in our simulations.

Because the evolution of the ring and the associated growth of spawned satellites occurs over $\geqslant {10}^{8}\,\mathrm{years}$, evolution of the satellite orbits due to tidal interaction with Saturn must be considered. Tides raised on Saturn by a satellite orbiting exterior (interior) to synchronous orbit produce a positive (negative) torque on the satellite's orbit, causing its orbit to expand (contract). The current synchronous orbit—where the orbital period equals Saturn's rotational day—lies within the Roche limit, with ${a}_{\mathrm{sync}}=1.9{R}_{\saturn}$. However early Saturn's radius was larger, and by conservation of angular momentum, it would have been rotating more slowly, with a_sync outside the Roche limit (Canup 2010). Tides raised on a satellite by the planet also modify the satellite's orbit, predominantly acting to decrease its orbital eccentricity.

We include the modification of satellite orbits due to tidal evolution in the N-body portion of our code by applying an additional accelerating "kick" to each orbiting object at every time step (Canup et al. 1999). We utilize the constant time delay tidal model of Mignard, which has a relatively straightforward analytic form and is valid for orbits near or that cross a_sync, and for high orbital eccentricities.

2.2.1. Planetary Tides

The acceleration of a satellite of mass m due to the second-order distortion it raises on the planet is given by (Mignard 1980; Touma & Wisdom 1994):

$$\begin{eqnarray}{\displaystyle \frac{{d}^{2}{\boldsymbol{r}}}{{{dt}}^{2}}|}_{p} & = & -\displaystyle \frac{3{k}_{2}{{GmR}}_{p}^{5}}{{r}^{10}}\left(1+\displaystyle \frac{m}{{M}_{p}}\right){\rm{\Delta }}t\\ & & \times [2({\boldsymbol{r}}\cdot {\boldsymbol{v}}){\boldsymbol{r}}+{r}^{2}({\boldsymbol{r}}\times {\boldsymbol{\omega }}+{\boldsymbol{v}})],\end{eqnarray} \tag{ 6 }$$

where ${\boldsymbol{r}}=(x,y,z)$ and ${\boldsymbol{v}}=({v}_{x},{v}_{y},{v}_{z})$ are the planetocentric satellite's position and velocity, k₂, is the planet's second-order Love number, and ${\boldsymbol{\omega }}=\omega {{\boldsymbol{u}}}_{z}$ is the planet's spin vector that we assume lies along the z-axis. The early value of k₂ is unknown; we adopt its current value for Saturn, ${k}_{2}=0.32$. The time lag ${\rm{\Delta }}t$ is defined as the time between the tide raising potential and when the equilibrium figure is achieved in response to this potential. The relation between the tidal time lag and the tidal dissipation factor Q is $Q\sim {(\psi {\rm{\Delta }}t)}^{-1}$ for a system oscillating at frequency ψ. For the planet, the dominant frequency is $\psi =2| \omega -n|$, where n is the satellite's mean motion, with ${\rm{\Delta }}t\sim 1/(2| \omega -n| Q)$.

2.2.2. Satellite Tides

The acceleration on a satellite due to tides raised by the planet on the satellite is (Mignard 1980)

$$\begin{eqnarray}{\displaystyle \frac{{d}^{2}{\boldsymbol{r}}}{{{dt}}^{2}}|}_{s} & = & -\displaystyle \frac{3{k}_{2}{{GmR}}_{p}^{5}}{{r}^{10}}\left(1+\displaystyle \frac{m}{{M}_{p}}\right){ \mathcal A }{\rm{\Delta }}t\\ & & \times [2({\boldsymbol{r}}\cdot {\boldsymbol{v}}){\boldsymbol{r}}+{r}^{2}({\boldsymbol{r}}\times {{\boldsymbol{\omega }}}_{{\boldsymbol{s}}}+{\boldsymbol{v}})],\end{eqnarray} \tag{ 7 }$$

where ${{\boldsymbol{\omega }}}_{{\boldsymbol{s}}}$ is the satellite's spin vector. The factor ${ \mathcal A }$ reflects the strength of satellite versus planetary tides, with

$$\begin{eqnarray}&&{ \mathcal A }={\left(\displaystyle \frac{m}{{M}_{p}}\right)}^{-2}{\left(\displaystyle \frac{{R}_{s}}{{R}_{p}}\right)}^{5}\left(\displaystyle \frac{{k}_{2s}}{{k}_{2}}\right)\left(\displaystyle \frac{{\rm{\Delta }}{t}_{s}}{{\rm{\Delta }}t}\right),\end{eqnarray} \tag{ 8 }$$

where k_2s, ${\rm{\Delta }}{t}_{s}$, and R_s are the satellite's Love number, tidal time lag, and physical radius.

The appropriate value for ${ \mathcal A }$ is very uncertain. In our simulations, $10\leqslant {({M}_{p}/m)}^{2}{({R}_{s}/{R}_{p})}^{5}\leqslant {10}^{2}$. Estimates suggest ${10}^{-3}\leqslant {k}_{2s}\leqslant {10}^{-1}$ for icy satellites (Murray & Dermott 1999, Table 4.1). For satellite tides, $\psi \approx n$ and ${\rm{\Delta }}{t}_{s}\sim 1/({Q}_{s}n)$, where we consider a satellite tidal dissipation factor ${Q}_{s}\sim {10}^{2}$. For $| \omega /n-1| \sim {10}^{-1}$, the final term in the expression above for ${ \mathcal A }$ is of order $10\leqslant ({\rm{\Delta }}{t}_{s}/{\rm{\Delta }}t)\leqslant {10}^{2}$ for ${10}^{4}\leqslant Q\leqslant {10}^{5}$. Thus the plausible range for ${ \mathcal A }$ is of order ${10}^{-1}\leqslant { \mathcal A }\leqslant {10}^{3}$. As such, we perform two sets of simulations that consider limiting cases: one without satellites tides (${ \mathcal A }=0$), and one with strong satellite tides (${ \mathcal A }=1000$).

When computing ${d}^{2}{\boldsymbol{r}}/{{dt}}^{2}{| }_{s}$, we make the simplifying assumption that satellites are rotating synchronously, so that ${\omega }_{s}\approx n=\sqrt{{{GM}}_{p}/{a}^{3}}$, where a is a semimajor axis. A non-synchronously rotating, uniform density satellite on a circular, non-inclined orbit will experience a torque $N=C{\dot{\omega }}_{s}$, where $C=(2/5){{mR}}_{s}^{2}$ is the satellite's moment of inertia and ${\dot{\omega }}_{s}$ is the time rate of change of its rotation rate, given by (Peale & Canup 2015)

$$\begin{eqnarray}\displaystyle \frac{d{\omega }_{s}}{{dt}} & = & -\displaystyle \frac{3{k}_{2s}{{GM}}_{p}^{2}{R}_{s}^{5}}{{{Ca}}^{6}}{\rm{\Delta }}{t}_{s}({\omega }_{s}-n)\\ & = & -\displaystyle \frac{15}{2}{k}_{2s}\displaystyle \frac{{M}_{p}}{m}\displaystyle \frac{{R}_{s}^{3}}{{a}^{3}}\displaystyle \frac{n}{{Q}_{s}}({\omega }_{s}-n).\end{eqnarray} \tag{ 9 }$$

If n is nearly constant, the quantity $({\omega }_{s}-n)$ decays exponentially with a time constant

$$\begin{eqnarray}{\tau }_{\mathrm{despin}} & = & \displaystyle \frac{2}{15}\displaystyle \frac{1}{{k}_{2s}}\displaystyle \frac{m}{{M}_{p}}\displaystyle \frac{{a}^{3}}{{R}_{s}^{3}}\displaystyle \frac{{Q}_{s}}{n}\approx 2\left(\displaystyle \frac{{Q}_{s}/{k}_{2s}}{{10}^{3}}\right)\left(\displaystyle \frac{m/{M}_{p}}{{10}^{-7}}\right)\\ & & \times {\left(\displaystyle \frac{a}{3{R}_{\saturn}}\right)}^{9/2}{\left(\displaystyle \frac{{250}^{}\mathrm{km}}{{R}_{s}}\right)}^{3}\,\mathrm{years}.\end{eqnarray} \tag{ 10 }$$

Thus ${\omega }_{s}$ will likely approach a synchronous value on a timescale short compared to orbital migration timescales.

We consider the evolution of a ring formed soon after the end of Saturn's gas accretion, at which time the planet will still be substantially larger than its current size due to the energy of its formation. To estimate the physical radius of Saturn, we use results from Fortney et al. (2007). The green dotted–dashed line in their Figure 5(B) represents the evolution of a 0.3 Jupiter mass planet, which is about the mass of Saturn, assuming a $25\,{M}_{\oplus }$ core. A fit to this data is shown as the green line in Figure 1. A core mass of $\sim 20\,{M}_{\oplus }$ (Hubbard et al. 2009) results in the red curve (data provided by W. Fortney for Canup 2010), an approximate fit to which is $R{(t)={A}_{0}+{A}_{1}\mathrm{log}{(t)+{A}_{2}\mathrm{log}(t)}^{2}+{A}_{3}\mathrm{log}(t)}^{3}$, where R is in units of R_♄, t is in years, ${A}_{0}\approx 9.576$, ${A}_{1}\approx -2.418$, ${A}_{2}\approx 0.231$, and ${A}_{3}\approx -0.0075$.

Zoom In Zoom Out Reset image size

From conservation of its spin angular momentum, one can estimate Saturn's early spin rate as a function of its physical radius and moment of inertia. We assume an early moment of inertia constant comparable to that of current Saturn, ${K}_{\saturn}\approx 0.23$ (Helled 2011; Nettelmann et al. 2013; see the Appendix). The resulting predicted evolution of the synchronous orbit with time is shown in Figure 2. While currently a_sync lies well inside the Roche limit, for the first  $\sim {10}^{9}\,\mathrm{years}$ of Saturn's history, synchronous orbit is shifted outward due to the slower rotation of the planet. For moons near the Roche limit, there will thus be a competition between the negative torque due to tides (causing orbital contraction) and the positive torque due to resonant interactions with the rings (causing orbital expansion). For $Q\geqslant {10}^{4}$ and $\sigma \geqslant {10}^{3}$ g cm⁻³, the latter are much stronger, allowing spawned satellites to evolve away from the rings.

Zoom In Zoom Out Reset image size

Table 2 lists parameters for our 12 baseline cases. We consider R_p = 1.3, 1.4, or $1.5\,{{\rm{R}}}_{\saturn}$, with corresponding planetary rotational periods of 17.9, 20.7, and $23.8\,\mathrm{hr}$, respectively. For orbits far beyond synchronous, $| \omega | \gg n$ and Saturn's tidal time lag is approximately ${\rm{\Delta }}t\approx 1/(2| \omega | Q)$. We set ${\rm{\Delta }}t$ using this expression so that $Q={10}^{4}$. We consider initial ring masses between $3\times {10}^{21}\,\mathrm{kg}$ and $1.1\times {10}^{22}\,\mathrm{kg}$, motivated by models of tidal stripping from a Titan-sized satellite (Canup 2010). We consider ice ring particles with density ${\rho }_{r}=0.9\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, which sets the Roche limit at a_R ≈ 2.24R_♄. Our initial ring extends from the planet's physical radius R_p to the Roche limit. It is possible that the ring may have been more concentrated initially, but it would rapidly viscously spread (Salmon et al. 2010).

Table 2.  Simulation Parameters

Run	Rp	T	async	aR	${\rm{\Delta }}t$	Md
	$({R}_{\saturn})$	(hr)	$({R}_{p})$	$({R}_{p})$	(s)	$({10}^{-5}\,{M}_{\saturn})$
1	1.5	23.8	2.19	1.50	0.68	0.5
2	1.5	23.8	2.19	1.50	0.68	1
3	1.5	23.8	2.19	1.50	0.68	1.5
4	1.5	23.8	2.19	1.50	0.68	2
5	1.4	20.7	2.14	1.60	0.59	0.5
6	1.4	20.7	2.14	1.60	0.59	1
7	1.4	20.7	2.14	1.60	0.59	1.5
8	1.4	20.7	2.14	1.60	0.59	2
9	1.3	17.9	2.09	1.73	0.51	0.5
10	1.3	17.9	2.09	1.73	0.51	1
11	1.3	17.9	2.09	1.73	0.51	1.5
12	1.3	17.9	2.09	1.73	0.51	2

Note. R_p is the planet's mean physical radius in units of Saturn current mean radius ${R}_{\saturn}=\mathrm{58,232}\,\mathrm{km}$. T is the spin period of the planet. a_sync and a_R are the position of the synchronous orbit and of the Roche limit, respectively, in units of the planet's physical radius. ${\rm{\Delta }}t$ is the tidal lag. M_d is the disk's initial mass.

Download table as:  ASCIITypeset image

We complete four sets of 12 baseline simulationsapjaa5ac2t2. The first and second sets assume no pre-existing exterior satellites, and consider either no satellite tides (${ \mathcal A }=0;$ "Set A") or strong satellite tides (${ \mathcal A }={10}^{3};$ "Set B"). The third and fourth sets include Dione and Rhea at their current locations, intended to represent the earlier formation of these outer moons as direct accretional products from the Saturnian subnebula, both with no satellite tides ("Set C") and with strong satellite tides ("Set D").

Figure 3 shows the system at different evolution times for the first ${10}^{8}\,\mathrm{years}$ for Run 6A. As the ring spreads, it starts producing new moonlets that, through resonant interaction, confine the ring inside the Roche limit. In turn, they recoil and the ring is progressively freed to viscously spread again. When a moonlet reaches ≥3.56R_♄, its 2:1 Lindblad resonance lies beyond a_R, and it thus stops directly interacting with the ring.

Zoom In Zoom Out Reset image size

When a moonlet is spawned at the Roche limit in the presence of exterior satellites, it will encounter their MMRs as it recoils outward due to ring torques. Initially, the ring is massive enough that its torques typically cause moonlets to recoil too rapidly for capture into resonance. As a result, as inner moonlets expand outward, they can have close encounters with outer satellites that can result in a merger and the growth of increasing massive moons. Figure 4 shows the masses of various objects in Run 6A, as a function of time. Colors correspond to indexes in our output mass array: black for moonlet #1 (which formed first and is the oldest), red for moonlet #2, green for moonlet #3, and purple for moonlet #4. Other bodies may be present at times, but have not been plotted for readability. Color changes occur when two objects merge. For example, at $t\approx {10}^{6}\,\mathrm{years}$, moonlet #3 (in green) merges with #2 (in red), and then moonlet #4 becomes the new moonlet #3, changing color from purple to green.

Zoom In Zoom Out Reset image size

Our simulations display the general behavior predicted in Crida & Charnoz (2012). Initially a moon spawned near the Roche limit directly accretes small ring material as it spreads across the Roche limit; this is defined as the "continuous regime" of growth (Crida & Charnoz 2012). As the moon rapidly recoils outward due to ring torques, its separation from the ring edge becomes large enough that a second inner satellite can begin to grow near the Roche limit. This second satellite also recoils outward and is eventually accreted by the outer satellite. This process repeats so long as the first satellite is relatively close to the ring's edge, with the first satellite growing at the same average rate as in the continuous regime, only through larger discrete steps; accordingly this mode of growth is called the "discrete regime" (Crida & Charnoz 2012). Finally, as the first satellite continues to evolve outward, it can become distant enough that it can no longer directly accrete a moon spawned from the ring, and a system of three or more moons results. Mergers are then characterized by collisions between similar-mass bodies in the so-called pyramidal regime (Crida & Charnoz 2012). The transition between the discrete and pyramidal regimes is predicted to occur when a moonlet of mass m reaches a distance r such that $r-2{r}_{H}\gt {r}_{c}$, where ${r}_{H}=r{(m/(3{M}_{\saturn}))}^{1/3}$ is the moonlet's Hill radius and ${r}_{c}={a}_{R}(8.4{M}_{\mathrm{disk}}/{M}_{\saturn}+1)$ (Crida & Charnoz 2012). In our simulations, ${M}_{\mathrm{disk}}/{M}_{\saturn}\sim {10}^{-5}$ such that moonlets should transition to the pyramidal regime after only little outward migration. This is indeed observed, but at somewhat larger distances than predicted by the previous expression. We find that several moonlets can be in the discrete regime at the same time (e.g., black and red curves before 10⁵ years on Figure 4). Also, we find that a younger moonlet can grow larger than an older one (e.g., black and red curves around 10⁵ and 10⁶ years in Figure 4). Configurations with an inner moon that is larger than an outer one are, however, generally only transient, as merging events between large objects eventually produce a system with larger satellites at larger distances, consistent with the Crida & Charnoz (2012) expectations.

As the ring mass decreases due to mass loss on the planet and by formation of moonlets at the Roche limit, two things can be noted. First, the ring's viscosity decreases and the time needed for the ring to spread back out to the Roche limit increases, such that the time between the spawning of new moonlets lengthens. Second, the Lindblad resonant torque becomes weaker and the orbital expansion of inner bodies slows, such that they can be captured into MMRs with outer bodies. When this happens, the inner object continues to recoil outward as it is torqued by the ring, and in turn it drives the outer object outward, as well due to the resonant configuration. This process allows for a transfer of angular momentum from the ring to outer objects that do not themselves have direct resonant interactions with the ring. This is a key process not included in the Charnoz et al. (2011) and Crida & Charnoz (2012) models.

The black line in Figure 5 shows the evolution of an object's semimajor axis in Run 6A. The object is initially spawned at the Roche limit at $t\sim {10}^{4}\,\mathrm{years}$ and moves outward due to resonant interactions with the rings. This object is not massive enough to confine the rings, such that secondary objects are subsequently spawned and also recoil (red and green curves in Figure 5). As they catch up with the outer object, mergers can occur and cause the outer object's semimajor axis to decrease somewhat due to the accreted object having a lower specific angular momentum (Salmon & Canup 2012). At 10⁸ years, the two most massive satellites in Run 6A have masses of $2.18\times {10}^{-6}\,{M}_{\saturn}$ and $8\times {10}^{-8}\,{M}_{\saturn}$, and semimajor axes of $3.94\,{R}_{\saturn}$ and $2.78\,{R}_{\saturn}$, respectively. They are similar to Tethys and Enceladus in mass, but their semimajor axes are smaller. Further expansion will be achieved over longer timescales due to tides and/or MMR interactions. This run did not produce a Mimas-equivalent satellite within ${10}^{8}\,\mathrm{years}$.

Zoom In Zoom Out Reset image size

Close encounters between satellites do not always result in a merger, as assumed in Crida & Charnoz (2012), but can instead lead to scattering, inward or outward. Such an event can be seen at ${10}^{4}\,\mathrm{years}$ in Run 6A (Figures 4 and 5, black and red lines). This leads to an orbital architecture in which the outermost body is less massive than the one immediately inside. This situation is transient, as the two moons re-exchange orbits (at $\sim 3\times {10}^{5}\,\mathrm{years}$ in this case).

Figure 6 shows the evolution of the rings in Run 6A: position of the outer edge (solid line), mass (dashed line), and mass fallen on the planet (dotted line). When a satellite is spawned at the ring's outer edge, resonant interactions cause the latter to slightly contract inside the Roche limit. Early on, the ring is still massive enough that its viscous torque is greater than the resonant torque from outer satellites, and the ring viscously spreads outward. As satellites grow larger through mutual collisions, and as the ring's mass decreases, the resonant torque can at times overwhelm the viscous torque, resulting in a prolonged contraction of the ring's outer edge (Figure 6, solid line at, e.g., ${10}^{6}\,\mathrm{years}$). As the confining satellite's orbit expands due to the resonant torque (Figure 5), the torque decreases both because the distance between the satellite and disk edge increases and because resonances with the satellite move outward with it and migrate out of the ring. Eventually, the ring viscously spreads outward again.

Zoom In Zoom Out Reset image size

Whenever the ring is confined inside the Roche limit, it continues to lose mass onto the planet through its inner edge, while not providing any additional mass to outer satellites. As a result, in this simulation more than $70 \%$ of the ring's mass is lost onto the planet. Mimas, Enceladus, and Tethys's masses total $1.35\times {10}^{-6}\,{M}_{\saturn}$, so that a ring initially at least three times as massive with $\gt 4\times {10}^{-6}\,{M}_{\saturn}$ would be necessary to produce these objects. In this run the rings still contain $\sim 2.5\times {10}^{-7}\,{M}_{\saturn}$ at 10⁸ years, which is about four times the mass of Mimas. Formation of a Mimas-equivalent may thus occur on longer timescales (see Section 4).

Figure 7 shows the distribution of satellites obtained in our 10⁸ year simulations without satellites tides and without pre-existing Dione and Rhea at different times of evolution. Spawned satellites have masses broadly comparable to those of Mimas, Enceladus, and Tethys. After $\sim {10}^{8}\,\mathrm{years}$, some objects have nearly reached the position of Tethys's orbit, primarily due to ring torques and capture into MMRs. The evolution of the semimajor axis of a satellite of mass m due to tides can be estimated by Burns (1977),

$$\begin{eqnarray}&&\displaystyle \frac{{da}}{{dt}}=\displaystyle \frac{3{k}_{2}}{Q}\displaystyle \frac{m}{{M}_{\saturn}}\sqrt{\displaystyle \frac{{{GM}}_{\saturn}}{{R}_{p}}}{\left(\displaystyle \frac{{R}_{p}}{a}\right)}^{11/2},\end{eqnarray} \tag{ 11 }$$

which can be integrated to give

$$\begin{eqnarray}&&a(t)={R}_{p}{\left[\displaystyle \frac{13}{2}\displaystyle \frac{3{k}_{2}}{Q}\displaystyle \frac{m}{{M}_{\saturn}}\sqrt{\displaystyle \frac{{{GM}}_{\saturn}}{{R}_{p}^{3}}}t+{\left(\displaystyle \frac{{a}_{0}}{{R}_{p}}\right)}^{13/2}\right]}^{2/13},\end{eqnarray} \tag{ 12 }$$

where a₀ is the satellite's initial semimajor axis. Panel (d) in Figure 7 shows a(m) from Equation (12) for the three values assumed for R_p, with $Q={10}^{4}$ and $t={10}^{8}$ years; the satellites produced in our simulations are shown with the same color scheme. Most satellites are to the right of their corresponding curve, indicating that they have orbits larger than expected due solely to tides. Satellites beyond $\sim 3.55\,{R}_{\saturn}$ have no resonances in the rings, and their orbital expansion beyond the dashed curves has been achieved by trapping the inner satellites into MMRs (i.e., by indirect angular momentum transport from the ring).

Zoom In Zoom Out Reset image size

Table 3 shows results from the Set A simulations. On average, they yield 2.6 ± 0.7 final satellites at ${10}^{8}\,\mathrm{years};$ at that time, the rings have an average mass of $3.83\times {10}^{-7}\,{M}_{\saturn}$, $4.78\times {10}^{-7}\,{M}_{\saturn}$, and $5.96\times {10}^{-7}\,{M}_{\saturn}$, for the runs with a planet radius of 1.5, 1.4, and $1.3\,{R}_{\saturn}$, respectively. Runs with larger planets have a lower ring mass at a given time because the flux onto the planet has been larger for a given initial ring mass. For each value of R_p, the ring mass is very similar at $t={10}^{8}$, even though the initial ring masses vary by a factor of 4. We find that about 20% of the ring's initial mass is incorporated into satellites, with a slightly higher fraction for smaller values of R_p. The average angular momentum of our satellites is higher than that of current Mimas, Enceladus, and Tethys (although standard deviation is important), mostly because our Enceladus and Tethys analogs tend to be a factor of a few times more massive than the current moons.

Table 3.  Set A Data at t = 10⁸ Years

	${R}_{p}=1.5{R}_{\saturn}$	${R}_{p}=1.4{R}_{\saturn}$	${R}_{p}=1.3{R}_{\saturn}$
$\langle {M}_{\mathrm{ring}}/{M}_{\mathrm{MET}}\rangle$	0.28 ± 0.01	0.36 ± 0.01	0.44 ± 0.001
$\langle {L}_{\mathrm{ring}}/{L}_{\mathrm{MET}}\rangle$	0.18 ± 0.01	0.22 ± 0.01	0.27 ± 0.002
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{MET}}\rangle$	1.85 ± 0.96	2.03 ± 1.1	2.11 ± 1.1
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{ring},0}\rangle$	0.2 ± 0.01	0.22 ± 0.01	0.23 ± 0.01
$\langle {L}_{\mathrm{sats}}/{L}_{\mathrm{MET}}\rangle$	1.7 ± 0.89	1.88 ± 1.04	1.98 ± 1.02

Note. Average values for our Set A runs at $t={10}^{8}$ years. M_ring is the rings' mass; M_MET is the total mass of Mimas, Enceladus, and Tethys; L_ring is the rings' angular momentum; M_sats is the total mass of the satellites in a given run; ${M}_{\mathrm{ring},0}$ is the rings' initial mass; L_sats is the angular momentum of the satellites in a given run; and L_MET is the total angular momentum of Mimas, Enceladus, and Tethys.

Download table as:  ASCIITypeset image

In all cases, sufficient mass and angular momentum remains in the rings to spawn additional satellites over longer timescales. The average angular momentum left in the rings is $4.4\times {10}^{32}\,\mathrm{kg}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$, $5.49\times {10}^{32}\,\mathrm{kg}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$, and $6.74\,\times {10}^{32}\,\mathrm{kg}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$ for the runs with a planet radius of 1.5, 1.4, and $1.3\,{R}_{\saturn}$, respectively. This is a few times larger than the angular momentum necessary to bring a Tethys-mass satellite from 4 to $5\,{R}_{\saturn}$ ($2.2\times {10}^{32}\,\mathrm{kg}\,{{\rm{m}}}^{2}\,{{\rm{s}}}^{-1}$). It appears then likely that the distribution of satellites will continue to evolve significantly on longer timescales, a point we return to in Section 4.

Eccentricities and inclinations of objects trapped in MMRs increase as the inner object is driven outward by the disk. The run shown in Figure 3 shows satellites with small eccentricities, but is a bit of an outlier in this regard compared to the satellites formed in the whole set. Damping of eccentricities occurs when objects collide, but significant values are generally reached in the absence of satellite tides. Bigger objects have on average smaller eccentricities and inclinations, since these objects have experienced more collisions, and it is more difficult to excite e and i for a bigger satellite. At $t={10}^{8}\,\mathrm{years}$, satellites with masses smaller than ${10}^{-6}\,{M}_{\saturn}$ have $\langle e\rangle \sim 0.073\pm 0.05$ (the median is $0.068$), while those with masses larger than ${10}^{-6}\,{M}_{\saturn}$ have $\langle e\rangle \sim 0.022\pm 0.02$ (median is $0.018$). Inclinations of satellites also become substantial, with an average of $\langle i\rangle \sim 3.32^\circ \pm 3.62^\circ$ (median is 1.95°) for small satellites, and $\langle i\rangle \sim 2.71^\circ \pm 5.14^\circ$ (median is 0.23°) for large satellites.

We perform a second set of simulations with the same initial parameters as in Table 2, but including satellite tides with ${ \mathcal A }=1000$. Figure 8 shows a system at different times of evolution. This case's evolution is similar to that in Figure 3, the main difference being the smaller eccentricities at all times, a direct consequence of satellite tides. Compared to the case without satellite tides, the largest satellite is smaller. This is due to a factor of 2 difference in the mass of the disk, which results in less massive satellites being formed and weaker orbital expansion rates, which overall decreases the delivery rate of material beyond the Roche limit.

Zoom In Zoom Out Reset image size

The overall results at $t={10}^{8}\,\mathrm{years}$ for the Set B runs are given in Table 4. Figure 9 shows the distribution of satellites at different times of evolution. The lower eccentricities of the satellites noted in Run 5B (Figure 8) can be observed across all runs. While without satellites tides some moons in Set A also have low eccentricities at $t={10}^{8}\,\mathrm{years}$, they went through a phase of high values before they were damped by accretional collisions (Figure 7).

Zoom In Zoom Out Reset image size

Table 4.  Set B Data at t = 10⁸ Years

	${R}_{p}=1.5{R}_{\saturn}$	${R}_{p}=1.4{R}_{\saturn}$	${R}_{p}=1.3{R}_{\saturn}$
$\langle {M}_{\mathrm{ring}}/{M}_{\mathrm{MET}}\rangle$	0.22 ± 0.07	0.33 ± 0.04	0.45 ± 0.005
$\langle {L}_{\mathrm{ring}}/{L}_{\mathrm{MET}}\rangle$	0.13 ± 0.04	0.20 ± 0.02	0.27 ± 0.004
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{MET}}\rangle$	2.19 ± 1.48	2.47 ± 1.69	2.69 ± 1.86
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{ring},0}\rangle$	0.22 ± 0.04	0.25 ± 0.05	0.27 ± 0.06
$\langle {L}_{\mathrm{sats}}/{L}_{\mathrm{MET}}\rangle$	2.01 ± 1.4	2.29 ± 1.59	2.51 ± 1.78

Note. Average values for our Set B runs at $t={10}^{8}$ years. M_ring is the rings' mass; M_MET is the total mass of Mimas, Enceladus, and Tethys; L_ring is the rings' angular momentum; M_sats is the total mass of the satellites in a given Run; ${M}_{\mathrm{ring},0}$ is the rings' initial mass; L_sats is the angular momentum of the satellites in a given Run; and L_MET is the total angular momentum of Mimas, Enceladus, and Tethys.

Download table as:  ASCIITypeset image

Figure 10 shows the mean eccentricities of the satellites as a function of time, weighted by the mass of the satellites, for the Set A versus Set B runs. The initial eccentricities of a fragment spawned at the Roche limit in our simulations is set to be approximately the ratio of the fragment's escape velocity to the local orbital velocity (Lissauer & Stewart 1993), which is $\sim {10}^{-3}$ for a ${10}^{-8}\,{M}_{\saturn}$ mass fragment. Set B runs with tidal dissipation in the satellites (dashed line); efficient damping occurs on a $\geqslant {10}^{6}$ year timescale, resulting in much smaller average eccentricities at ${10}^{8}\,\mathrm{years}$, compared to our runs without satellites tides (solid line).

Zoom In Zoom Out Reset image size

Compared with Set A, there is a somewhat higher average fraction of the ring's initial mass incorporated into satellites, which in Set B approaches 30%. On average, the Set B simulations have 3.75 ± 1.05 satellites at ${10}^{8}\,\mathrm{years}$. A larger final number of spawned moons is a consequence of the satellites' lower eccentricities: their pericenter is larger, and they do not "sweep" the region close to the rings as efficiently, allowing more small moons to form and survive. The combination of more mass put into satellites and lower eccentricities contributes to a higher total angular momentum in spawned moons compared with Set A.

For a traditional Saturn tidal parameter of $Q\sim {10}^{4}$, mid-sized moons spawned from a Roche-interior ring do not reach distances consistent with those of Dione and Rhea. We assume Dione and Rhea formed from a different process, which we argue is the most straightforward way to explain the much larger total mass of rock in these moons compared to that in Mimas, Enceladus, and Tethys. If Dione and Rhea were present as the inner mid-sized moons (or their progenitors) were spawned from the rings, the inner moons would have encountered MMRs with the outer satellites as the spawned moons recoiled outward due to ring interactions. For example, Dione's 2:1 MMR currently lies near $4.08\,{R}_{\saturn}$ and Rhea's 3:1 MMR currently lies around $4.35\,{R}_{\saturn}$, positions that a Tethys analogue would need to cross to reach Tethys's current orbit at $5.06\,{R}_{\saturn}$.

In the Set C runs, we include Dione and Rhea, with their current masses and positions with ${ \mathcal A }=0$ (no satellite tides); results are shown in Table 5 and Figure 11.

Zoom In Zoom Out Reset image size

Table 5.  Set C Data at t = 10⁸ Years

	${R}_{p}=1.5{R}_{\saturn}$	${R}_{p}=1.4{R}_{\saturn}$	${R}_{p}=1.3{R}_{\saturn}$
$\langle {M}_{\mathrm{ring}}/{M}_{\mathrm{MET}}\rangle$	0.28 ± 0.01	0.36 ± 0.01	0.45 ± 0.01
$\langle {L}_{\mathrm{ring}}/{L}_{\mathrm{MET}}\rangle$	0.17 ± 0.01	0.22 ± 0.01	0.27 ± 0.01
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{MET}}\rangle$	1.88 ± 1.0	2.02 ± 1.05	2.15 ± 1.16
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{ring},0}\rangle$	0.2 ± 0.01	0.22 ± 0.01	0.23 ± 0.01
$\langle {L}_{\mathrm{sats}}/{L}_{\mathrm{MET}}\rangle$	1.72 ± 0.94	1.87 ± 1.01	1.99 ± 1.1
$\langle {\rm{\Delta }}{L}_{\mathrm{DR}}/{L}_{\mathrm{MET}}\rangle$	0.02 ± 0.01	0.01 ± 0.01	0.01 ± 0.01

Note. Average values for our Set C runs at $t={10}^{8}$ years. M_ring is the rings' mass; M_MET is the total mass of Mimas, Enceladus, and Tethys; L_ring is the rings' angular momentum; M_sats is the total mass of the spawned satellites in a given run (excluding Dione and Rhea); ${M}_{\mathrm{ring},0}$ is the rings' initial mass; L_sats is the angular momentum of the spawned satellites in a given run (excluding Dione and Rhea); L_MET is the total angular momentum of Mimas, Enceladus, and Tethys; and $\langle {\rm{\Delta }}{L}_{\mathrm{DR}}\rangle /{L}_{\mathrm{MET}}$ is the variation of angular momentum of Dione and Rhea.

Download table as:  ASCIITypeset image

The inclusion of Dione and Rhea does not dramatically alter the distribution of spawned satellites, and most results here are similar to those in Set A. The average number of satellites per run at 10⁸ years is 2.08 ± 0.3 (not including Dione and Rhea). No spawned moon collides with either Dione or Rhea in any of the set C runs. However, some of the spawned satellites do become captured into MMR with Dione and Rhea, as expected. These resonant configurations can be transient or still present at 10⁸ years. This has two effects. First, the inner satellites transfer some angular momentum to Dione and Rhea, resulting in Tethys analogs that have somewhat smaller semimajor axes at 10⁸ years compared to Sets A and B. Second, in some cases, Dione and Rhea experience substantial eccentricity growth, a result of MMRs with inner moons in the absence of tidal dissipation in the satellites. At 10⁸ years, Dione and Rhea have mean eccentricities of 0.02 ± 0.01 and 0.01 ± 0.01, respectively.

Set D simulations begin with Dione and Rhea in their current positions, and include tidal dissipation in the satellites with ${ \mathcal A }={10}^{3}$. Table 6 and Figure 12 show results at $t={10}^{8}\,\mathrm{years}$. Set D simulations produce 3.4 ± 0.5 spawned satellites at ${10}^{8}\,\mathrm{years}$ (excluding Dione and Rhea). Compared with the case without satellites tides (Set C), we find that less mass from the ring is placed into satellites by 10⁸ years, and the spawned satellites have lower total angular momentum. In Set D, tidal damping of eccentricities keeps the eccentricities of Dione and Rhea lower (of order 10⁻³ on average), and maintains some of the spawned moons in MMRs with Dione and Rhea. We find that the largest spawned moon is captured in Dione's 2:1 MMR in 9 of the 12 runs from Set D, with the resonance configuration still present at 10⁸ years. As inner spawned moons transfer more of the ring's angular momentum to Dione and Rhea, the spawned moon orbits do not expand as far as in the case without satellite tides, resulting in greater confinement of the rings and somewhat less total mass incorporated into the spawned moons.

Zoom In Zoom Out Reset image size

Table 6.  Set D Data at t = 10⁸ Years

	${R}_{p}=1.5{R}_{\saturn}$	${R}_{p}=1.4{R}_{\saturn}$	${R}_{p}=1.3{R}_{\saturn}$
$\langle {M}_{\mathrm{ring}}/{M}_{\mathrm{MET}}\rangle$	0.26 ± 0.05	0.34 ± 0.03	0.4 ± 0.05
$\langle {L}_{\mathrm{ring}}/{L}_{\mathrm{MET}}\rangle$	0.16 ± 0.03	0.21 ± 0.02	0.24 ± 0.03
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{MET}}\rangle$	1.66 ± 0.94	1.84 ± 1.01	2.04 ± 1.12
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{ring},0}\rangle$	0.18 ± 0.01	0.19 ± 0.01	0.22 ± 0.01
$\langle {L}_{\mathrm{sats}}/{L}_{\mathrm{MET}}\rangle$	1.49 ± 0.86	1.67 ± 0.95	1.85 ± 1.05
$\langle {\rm{\Delta }}{L}_{\mathrm{DR}}/{L}_{\mathrm{MET}}\rangle$	0.03 ± 0.01	0.02 ± 0.01	0.02 ± 0.02

Note. Average values for our Set D runs at $t={10}^{8}$ years. M_ring is the rings' mass; M_MET is the total mass of Mimas, Enceladus, and Tethys; L_ring is the rings' angular momentum; M_sats is the total mass of the satellites in a given run (excluding Dione and Rhea); ${M}_{\mathrm{ring},0}$ is the rings' initial mass; L_sats is the angular momentum of the satellites in a given run (excluding Dione and Rhea); L_MET is the total angular momentum of Mimas, Enceladus, and Tethys; and $\langle {\rm{\Delta }}{L}_{\mathrm{DR}}\rangle /{L}_{\mathrm{MET}}$ is the variation of angular momentum of Dione and Rhea.

Download table as:  ASCIITypeset image

Figure 12 shows the obtained distribution of satellites, at different times of evolution. As with the case without Dione and Rhea, the inclusion of tidal dissipation in the satellites efficiently damps the eccentricities of the growing moons. As a consequence, the radial "feeding" zone of each satellite is narrower, as they remain on quasi-circular orbits, which allows on average a larger number of satellites to survive, in particular Mimas progenitors lying close to the rings.

We ran the 12 Set A simulations for ${10}^{6}\,\mathrm{years}$ with the accelerated code and compared the obtained distribution of satellites with that obtained over ${10}^{7}\,\mathrm{years}$ with the standard code. The results are shown in Figure 13. The average number of satellites per run is 3.17 ± 1.03 with the accelerated code, compared to 2.42 ± 0.79 with the standard code. This is mainly due to a greater number of small satellites in the accelerated code. If we consider only satellites with a mass $\geqslant 2\times {10}^{-8}\,{M}_{\saturn}$, thereby removing the stochasticity of newly spawned moonlets, then the average number of satellites per run is 2.25 ± 0.96 with the standard code, and 2.58 ± 0.79 with the accelerated code. The average eccentricity of satellites is 0.05 ± 0.066 for the accelerated runs, and 0.05 ± 0.058 for the normal ones. Thus the overall distribution of spawned satellites produced by our accelerated code is similar to those obtained with the normal code.

Zoom In Zoom Out Reset image size

For our initial runs, we have assumed that the planet radius remained constant over ${10}^{8}\,\mathrm{years}$, with a value of 1.3, 1.4, or 1.5R_♄. For our billion-year runs, we use as our starting condition the outputs of Runs 5 to 12 in Sets A through D. Thus the initial planetary radius is either 1.3 or $1.4\,{R}_{\saturn}$. We then include the contraction of the planet and increase its spin (which moves the synchronous orbit inward) in the 10⁹ year runs. This is done by first estimating the initial age t_ev of the planet, given its radius (i.e., inverting Figure 1), and then by computing the new planetary radius at time t given by $R({t}_{{ev}}+t)$. To save computational time, we do not update the radius at every time step. We find that updating the radius every $\sim {10}^{4}\,\mathrm{years}$ gives a good approximation of our analytical derivation from Section 2.3.

Figure 14 shows the later evolution of the system from Run 6A on the left panel and for Run 5B on the right panel. In both cases the largest satellite spawned in the first 10⁸ years does not accrete additional mass but keeps evolving away due to tides and capture of inner moons in MMRs. The second largest satellite has continued growing, and reaches its final mass after a few 10⁸ years. The formation of a Mimas-equivalent satellite is not complete in Run 6A. In Run 5B, there are two moons around the position of current Mimas, and they will likely merge over longer timescales. This is also seen in the other runs where we get good Enceladus and Tethys equivalent satellites. On the other hand, cases where we form a good Mimas analogue at 10⁹ years have either too many moons outside Mimas, or two moons that are much more massive than Enceladus and Tethys.

Zoom In Zoom Out Reset image size

In Run 6A, the rings still contain $\sim 1.4\times {10}^{-7}\,{M}_{\saturn}$ and should be able to provide enough material to complete the formation of a Mimas-type satellite over longer timescales. A similar conclusion can be made for the other runs that have good Enceladus and Tethys equivalents at 10⁹ years. However,this implies that Mimas may be at least a billion years younger than Tethys. A similar point was made by Charnoz et al. (2011), who claimed that "a Mimas-like satellite could be about 1–1.5 Gy younger than a Rhea-like satellite."

Figure 15 shows the distribution of satellites at ${10}^{9}\,\mathrm{years}$ for all cases considered previously (with/without Dione and Rhea, and with/without satellite tides). For the cases without Dione and Rhea (panels a and c), the distribution of satellites obtained in our various runs agrees reasonably well with the masses and positions of current Mimas, Enceladus, and Tethys. The runs that include tidal dissipation in the satellites (panels c and d) are better at producing Mimas equivalents. This is due to the fact that tidal dissipation keeps the eccentricity of the outer satellites small, such that their pericenter lies further from the edge of the rings, preventing them from "sweeping" this area. The average number of satellites per run at 10⁹ years is 3.25 ± 0.7 for Set A, 4.87 ± 0.64 for Set B, 2.12 ± 1.13 for Set C, and 3.75 ± 0.71 for Set D. As expected, the number of satellites per run is larger in cases that include tidal dissipation in the satellites, as the latter keeps eccentricities low and limits the feeding zone of a given satellite, allowing more satellites to coexist.

Zoom In Zoom Out Reset image size

For the set of runs that include Dione and Rhea (panels b and d), the match to Mimas, Enceladus, and Tethys is best when we include satellite tides (Set D). Through capture into MMR, the inner satellites have transferred some of their angular momentum to Dione and/or Rhea, and have thus experienced weaker orbital migration. On some cases, Dione and Rhea have migrated away significantly. As for the runs without Dione and Rhea, the inclusion of tidal dissipation into the satellites allows the survival of Mimas-like satellites.

Table 7 lists average quantities for the rings and produced satellites at 10⁹ years across our four sets of simulations. The average mass of the rings at ${10}^{9}\,\mathrm{years}$ is $1.5\times {10}^{-7}\,\pm 2\times {10}^{-8}\,{M}_{\saturn}$ across all our runs. Despite a factor of six variation in the initial mass, all rings have roughly the same mass at the end of our simulations, in good agreement with the expected asymptotic evolution of the ring mass (Salmon et al. 2010). The average ring mass at ${10}^{9}\,\mathrm{years}$ is only slightly larger than the mass of Mimas, in good agreement with current estimates of the mass of Saturn's rings (Esposito et al. 1983; Salmon et al. 2010).

Table 7.  Data for All 4 Sets at t = 10⁹ Years

	Set A	Set B	Set C	Set D
$\langle {M}_{\mathrm{ring}}/{M}_{\mathrm{MET}}\rangle$	0.12 ± 0.02	0.12 ± 0.01	0.12 ± 0.02	0.11 ± 0.02
$\langle {L}_{\mathrm{ring}}/{L}_{\mathrm{MET}}\rangle$	0.07 ± 0.01	0.07 ± 0.01	0.07 ± 0.01	0.06 ± 0.01
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{MET}}\rangle$	2.12 ± 1.02	2.62 ± 1.66	2.13 ± 1.02	1.95 ± 1.01
$\langle {M}_{\mathrm{sats}}/{M}_{\mathrm{ring},0}\rangle$	0.22 ± 0.01	0.25 ± 0.05	0.22 ± 0.01	0.2 ± 0.01
$\langle {L}_{\mathrm{sats}}/{L}_{\mathrm{MET}}\rangle$	2.18 ± 1.1	2.67 ± 1.78	2.2 ± 1.13	1.87 ± 1.03
$\langle {\rm{\Delta }}{L}_{\mathrm{DR}}/{L}_{\mathrm{MET}}\rangle$	N/A	N/A	0.09 ± 0.1	0.11 ± 0.08

Note. Average values at 10⁹ for the four sets of simulations. M_ring is the rings' mass; M_MET is the total mass of Mimas, Enceladus, and Tethys; L_ring is the rings' angular momentum; M_sats is the total mass of the spawned satellites in a given run (excluding Dione and Rhea); ${M}_{\mathrm{ring},0}$ is the rings' initial mass; L_sats is the angular momentum of the spawned satellites in a given run (excluding Dione and Rhea); L_MET is the total angular momentum of Mimas, Enceladus, and Tethys; and $\langle {\rm{\Delta }}{L}_{\mathrm{DR}}/{L}_{\mathrm{MET}}\rangle$ is the variation of angular momentum of Dione and Rhea.

Download table as:  ASCIITypeset image

Figure 16 shows the range of satellites formed across all our simulations. We plot the average mass and semimajor axis of our Tethys, Enceladus, and Mimas analogs. For the latter, we sum the masses of the third largest bodies with any other moonlets present inside its orbit, as they will likely collide and merge on longer timescales, and we compute a mass-weighted mean for the semimajor axis. For our outermost satellite, the mass range is similar in sets A, B, and C, and a little lower for set D. For the outermost satellite, sets not including tidal dissipation in the satellites (A and C) have consistently larger average semimajor axes than the corresponding set that includes them (B and D). Our Enceladus analogs are systematically larger than current Enceladus, which is consistent with later mass loss from this body (e.g., Hansen et al. 2008).

Zoom In Zoom Out Reset image size

For the innermost satellite there is a clear separation between sets that do not include satellite tides, which produce smaller and closer moons, and sets that include strong satellite tides, which produce larger and more distant satellites. While we have in this study explored extreme cases (no satellite tides or very strong ones), this suggests that moderate satellite tides could produce an innermost satellite more similar to Mimas.

We find that many of our produced satellite systems show resonant configurations at 10⁹ years: three runs in Set A, five runs in Set B, three runs in Set C, and six runs in Set D. For Set C and D, a resonant configuration between Dione and Rhea occurs in two and five runs, respectively. In Set D, the largest moon formed from the rings (i.e., our Tethys equivalent) is captured in a Dione's 2:1 MMR in five runs. Despite efficient capture into MMRs in all our cases, survival of the system of Saturn's moons through their formation appears likely over 1 billion years. Some tidal dissipation into the growing satellites seems a prerequisite to allow the formation of the innermost satellites, in particular Mimas.

In the so-called Nice model for the origin of the dynamical structure of the outer solar system (Tsiganis et al. 2005), the giant planets are initially in a compact orbital configuration. Their orbits slowly migrate and diverge due to dynamical interactions with a planetesimal disk of initial mass M_disk. When Jupiter and Saturn cross their mutual 2:1 MMR, their orbital eccentricities increase. This leads to a period of dynamical instability, during which Uranus and Neptune (and perhaps a fifth outer planet as well; Nesvorný & Morbidelli 2012) are scattered outward to their current positions and planetesimals are scattered across the solar system. The scattered planetesimals become a population of impactors, which could, for example, account for the spike in the impact rate believed by many to be necessary to explain the formation of the large lunar impact basins at ∼3.9 Gyr during a so-called late heavy bombardment on the Moon (Tera et al. 1974; Cohen et al. 2000; Stöffler & Ryder 2001; Kring & Cohen 2002; Gomes et al. 2005; Levison et al. 2011). The total mass of scattered planetesimals scales roughly with M_disk. The disk must be massive enough to decrease the eccentricities of Uranus and Neptune to their current values through friction effects, but not so massive as to cause excessive migration. Forming a planetary system like ours appears most likely for $20\,{M}_{\oplus }\leqslant {M}_{\mathrm{disk}}\leqslant 50\,{M}_{\oplus }$ (Batygin & Brown 2010; Nesvorný & Morbidelli 2012). For ${M}_{\mathrm{disk}}=35\,{M}_{\oplus }$, $\sim {10}^{19}\,\mathrm{kg}$ of planetesimals are scattered onto the Moon, consistent with that needed to explain the lunar LHB (Gomes et al. 2005). If the instability was responsible for the lunar LHB, the timing of the event is constrained to occur some 700 Myr after the origin of the solar system.

A Nice-like instability consistent with our specific solar system structure may be a relatively improbable case among the broad range of possible outcomes (Nesvorný & Morbidelli 2012). However, the Nice model remains the most detailed dynamical history available for the early outer solar system, and it has been remarkably successful in explaining a variety of solar system features, including the giant planet eccentricities (Tsiganis et al. 2005), the structure of the Kuiper Belt (Levison et al. 2008), the capture of Jupiters Trojans (Morbidelli et al. 2005), the capture of the irregular satellites (Nesvorný et al. 2007), and the Ganymede/Callisto dichotomy (Barr & Canup 2010). While we consider predictions of the specific Nice model as a guide, the existence of an enhanced bombardment period in the outer solar system probably applies more generally. It has long been recognized that the structure of the Kuiper Belt requires that Neptune migrated outward via planetesimal scattering (Malhotra 1995), and that this implies both an initially more compact giant planet configuration and a planetesimal disk containing between 10 and 100M_⊕ (Fernandez & Ip 1984; Hahn & Malhotra 1999). The interaction of the giant planets with such a massive disk would have produced an enhanced bombardment period, even if the details of the evolution differed from that of the Nice model.

The composition of the impactors originating from the planetestimal disk is uncertain. Jupiter Trojans are thought to be captured bodies that originated in the region beyond Neptune (Nesvorný et al. 2013). These bodies show a diversity of densities: $0.8\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ for 617 Patroclus (Marchis et al. 2006), and $2.5\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ for 624 Hektor (Lacerda & Jewitt 2007). Jupiter-family comets are also thought to originate from the Kuiper Belt (Fernandez 1980; Duncan et al. 1988), and show densities lower than $0.6\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ (Lowry et al. 2008). These low densities indicate high porosities. The rock fraction of these objects would be between about 40% and 70% if they reflected a bulk solar composition of outer solar system solids (Simonelli et al. 1989). Observations from the Deep Impact mission imply that the rock-to-ice ratio in comet 9P/Tempel 1 is larger than 1, suggesting that comets are "icy dirtballs" rather than "dirty iceballs" (Küppers et al. 2005). The size distribution of impactors is also uncertain, but constraints can be derived from the cratering rate on Iapetus, and the observed size distribution of comets and KBOs (Charnoz et al. 2009).

An LHB in the outer solar system would have affected any satellite of Saturn that existed at that time. We here consider a bombardment that occurs at 700 Myr, and estimate the mass of rock that would have collided with Saturn's inner moons as a function of the assumed mass of the transneptunian disk, assuming that the LHB impactors contain between 40% and 60% rock by mass.

For an initial disk mass ${M}_{d,0}=35\,{M}_{\oplus }$, the Moon accretes an average of $(8.4\pm 0.3)\times {10}^{18}\,\mathrm{kg}$ (Gomes et al. 2005), which is $\sim 4\times {10}^{-8}\,{M}_{{\rm{d}},0}$. Levison et al. (2001) found that during and LHB-type event, Callisto and Ganymede are respectively impacted 40 and 110 times more than the Moon. Using the impact probabilities from Table 1 of Zahnle et al. (2003), this implies that a mass $\sim 3.14\times {10}^{-2}\,{M}_{{\rm{d}},0}$ collides with Jupiter, while $\sim 1.32\times {10}^{-2}\,{M}_{{\rm{d}},0}$ collides with Saturn, where these values assume the planets have their current mean radii.

Given the probability of impact with Saturn, P_♄, the impact probability onto a Saturnian satellite, P_s, can be approximated by (Zahnle et al. 2003)

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{s}}{{P}_{\saturn}}\approx \displaystyle \frac{{R}_{s}^{2}}{{R}_{\saturn}{a}_{s}},\end{eqnarray} \tag{ 13 }$$

where R_s and a_s are the satellite's physical radius and its semimajor axis.

We apply this formula to the the distribution of moons in each of our simulations at t = 700 Myr to determine the probability of impact with each moon. For the spawned moons, we calculate R_s, assuming a density appropriate for pure ice; for Dione and Rhea, we use their current physical radii. We estimate the total mass of rock delivered to the satellites during the LHB by computing the total of these probabilities times ${M}_{d,0}$. Figure 17 shows results as a function of ${M}_{d,0}$. The points indicate the average mass of rock (and its standard deviation) by satellites spawned from the rings (left panel), and by Dione and Rhea in runs that included them (right panel). The vertical dashed lines represent variations for impactors containing 40%–60% rock.

Zoom In Zoom Out Reset image size

The mass of rock that impacts the satellites spawned from the rings is consistent with the total mass of rock in Mimas, Enceladus, and Tethys (shown as the horizontal dashed lines in Figure 17) for ${M}_{d,0}\gt 20\,{M}_{\oplus }$. Thus if an LHB in the outer solar system occurred at roughly the same time as the lunar LHB, it would have delivered a rock mass comparable to the total rock in Saturn's inner moons. This suggests that these moons were initially much more rock-poor than they are today. In contrast, the LHB delivers a mass in rock to Dione and Rhea that is about an order-of-magnitude less than the current rock content in these moons (Figure 17, right panel), suggesting that they were already (relatively) rock-rich when they formed.

In the above we estimate the total mass of rock from external impactors that collides with the spawned satellites as a group, but we have not calculated the rock mass that would ultimately be accreted by each satellite. Using the mean values for each equivalent satellite (Figure 16), we can estimate the expected rock mass to collide with each object. We find that "Mimas" and "Tethys" would be impacted by their currently estimated mass of rock, provided the initial mass of the transneptunian disk was $\gt 20\,{M}_{\oplus }$, but that "Enceladus" would be impacted by about a factor of two, with too little rock mass. We envision several possibilities to explain this discrepancy. First, most of the mass was delivered by large impactors with radius $\sim 100\,\mathrm{km}$ or larger (Charnoz et al. 2009). We estimate that stochastic variations could then produce the observed rock distribution in the system with a probability up to 10%. Second, Enceladus may have initially formed with more rock due to the presence of some large chunks in the rings, as suggested by Charnoz et al. (2011). Third, the amount of rock colliding with each moon was different than the rock ultimately retained by each moon.

Per the third possibility, impact velocities onto the spawned moons will be dominated by gravitational focusing by Saturn, and so will greatly exceed each moon's escape velocity. For a typical impact velocity ∼20 km s⁻¹, the impactor will be destroyed and its rock shock heated to temperatures ∼2000–8000 K, implying a primarily melt-vapor state for all but the most highly oblique impacts (Pierazzo & Melosh 2000). While much of the ejecta in a hypervelocity cratering event (e.g., that originated from the target in the"far field") has ejection velocities comparable to the target's escape velocity (Alvarellos et al. 2005), the impactor material itself will have an ejection velocity within a factor of several of the impact velocity (Melosh 1989; Pierazzo & Melosh 2000). Thus most or all of the impacting material will initially escape the satellite, although most will still be in bound Saturn orbit (Movshovitz et al. 2015). For small ejecta sizes, as would be expected for droplets condensing from vapor, ejecta–ejecta collisions may occur on much smaller timescales than re-accretion onto the target, allowing ejecta to collisionally damp to a ring that overlaps the orbits of neighboring satellites. Thus while collisions with Saturn's inner satellites should effectively capture rock from external impactors into the inner satellite region, where the rock ultimately accretes will depend on the ejecta's size and velocity distribution and on its post-impact evolution. This will require follow-on models for assessment.

External impactors would have also encountered the rings. Large impactors would pass through the rings, while those small enough to encounter a ring mass comparable to their own during a single passage could be directly captured. For a ring at 700 Myr with a surface density ∼few ${10}^{2}\,{\rm{g}}\,{\mathrm{cm}}^{-2}$, this corresponds to impactor radii smaller than about $2\,{\rm{m}}$. The size distribution of small impactors during the LHB is uncertain, but it is inferred to be quite shallow, with a cumulative size index of −1.5, based on the cratering record on Iapetus (Charnoz et al. 2009; Movshovitz et al. 2015). With 0.002 times the total impactor population expected in objects $\leqslant 1\,\mathrm{km}$ in radius (Movshovitz et al. 2015), the fraction in objects smaller than $2\,{\rm{m}}$ would be $\sim {10}^{-7}$. This implies a relatively small total rock mass captured by the rings, comparable to or smaller than the upper limit on the rock in the rings today (Nicholson et al. 2008).