A Recent Impact Origin of Saturn's Rings and Mid-sized Moons

SPH (Gingold & Monaghan 1977; Lucy 1977) is a Lagrangian method widely used to model planetary impacts and diverse astrophysical and other systems. Materials are represented by many interpolation points or "particles" that are evolved under hydrodynamical forces, with pressures computed via an equation of state (EOS), and gravity. SPH is particularly well suited for modeling asymmetric geometries and/or large dynamic ranges in density and distribution of material, such as the messy impacts considered here. However, standard resolutions in planetary simulations, using 10⁵–10⁶ SPH particles, can fail to converge on even large-scale outcomes of giant impacts (Genda et al. 2015; Hosono et al. 2017; Kegerreis et al. 2019, 2022).

Here we run simulations with up to ∼10^7.5 SPH particles, to resolve the detailed outcomes of collisions and to test for numerical convergence, using the open-source hydrodynamics and gravity code SWIFT ⁷ (Kegerreis et al. 2019). For this study, we use a "vanilla" form of SPH plus the Balsara (1995) switch for the artificial viscosity. Future work will include more sophisticated modifications recently added to SWIFT that mitigate issues that can arise from density discontinuities in SPH (Ruiz-Bonilla et al. 2022; T. D. Sandnes et al. 2023, in preparation). As in previous works, we also neglect material strength, as gravitational stresses are expected to dominate (Asphaug & Reufer 2013), though, as smaller bodies are considered, this assumption will be valuable to test in future work. Nevertheless, we tentatively expect these standard simplifications to be comparatively innocuous for the high-speed collisions (relative to the mutual escape speed) considered here (Asphaug & Reufer 2013; Kegerreis et al. 2020).

To distinguish the Dione-/Rhea-like moons we collide from their present-day analogs, we refer to the precursor satellites as "p-Dione" and "p-Rhea." While a precursor system may have been larger in mass, we start by considering moons like those today, with masses of 1.10 × 10²¹ kg and 2.31 × 10²¹ kg, respectively. Their detailed internal state cannot readily be predicted, so, for simplicity and in line with previous works, we model them both as differentiated into silicate cores and ice mantles, at an isothermal 100 K, with core mass fractions of 60% and 40%, respectively. These materials are modeled with Tillotson (1962) basalt and water EOS (Melosh 2007), or with ANEOS forsterite and AQUA water (Haldemann et al. 2020; Stewart et al. 2020) for a small subset of comparison simulations.

We generate the moons' internal profiles by integrating inward while maintaining hydrostatic equilibrium, ⁸ and then we place the roughly equal-mass SPH particles to match the resulting density profiles using the stretched equal-area method. ⁹ Before simulating the impact, we first run a brief 2 hr settling simulation for each body in isolation, to allow any final settling to occur. The specific entropies of the particles are kept fixed, enforcing that the particles relax themselves adiabatically. The impact simulations are then evolved from 1 hr before contact to 9 hr after. As this is only a modest fraction of the orbital period around Saturn (∼82 hr), we neglect the effects of Saturn's gravity during the SPH simulations, as discussed in Section 2.4.

As detailed in Ć16, a pair of precursor moons can be suddenly excited by p-Rhea encountering the evection resonance, leading to high eccentricities, inclinations, and crossing orbits. We take two specific example outcomes of Ć16 models, which result in collisions with relative velocities of 2 and 3 km s⁻¹. The orbital parameters are listed in Section 2.4. For comparison, the mutual escape speed is v_esc ≈ 580 m s⁻¹, so the impacts are at ∼3.4v_esc and 5.1v_esc. These two orbit scenarios set the velocity at infinity, v_∞, for the impact initial conditions. Given the moons' small masses and high relative velocities, little additional speed (∼20 m s⁻¹) is gained from self-gravity as the moons approach.

The moons' orbits are rapidly and chaotically varying at this time, and their radii are a negligible fraction of their orbital semimajor axes, so we freely set a variety of impact angles to examine the range of possible outcomes at these speeds. For each of the two speed scenarios, we simulate collisions at eight impact angles (see Figure 1) from nearly head-on to highly grazing, with β = 5°–75° in steps of 10°. This angle positions the point of contact along a line from the center to the edge of the target moon. The impact scenario can also be rotated by an angle ϕ (default 0°) about the x-axis in the SPH frame to position the point of contact anywhere on the target. While a ϕ ≠ 0 simulation is identical in the SPH frame, the orientation in Saturn's frame and thus the resulting distribution of orbiting debris are affected. Therefore, we repeat the analysis of the debris with ϕ = 90°, 180°, and 270° for the primary suite of simulations.

For each angle and speed scenario, the SPH particle mass, which controls the resolution, is set such that the number of particles per M_Rh + M_Di = 3.4 × 10²¹ kg is from 10⁵ up to 10^7.5 in steps of 0.5 dex. In addition to the primary suite of p-Dione–p-Rhea impacts, we repeat these same simulations for p-Rhea–p-Rhea as considered by Hyodo & Charnoz (2017), although they used nonequal ratios between rock and ice.

For a cleaner side test of how the results might scale with the precursor moons' masses, we also run a subset of impacts with p-Dione masses of (1, 2, 4, 8) × 10²¹ kg and corresponding p-Rhea masses of double these values. These moons are given a core mass fraction of $\tfrac{1}{3}$. We simulate these collisions for β = 15°, 25°, and 35° for the four mass pairs.

Finally, to further examine the numerical effects of our finite-particle models, each scenario at our second-highest resolution of 10⁷ particles is repeated four more times using a rotated orientation for the settled p-Rhea and p-Dione. Particularly for low-speed collisions, this can yield nonnegligible effects even at the high resolutions used here (Kegerreis et al. 2020, 2022)—an uncertainty that these tests will help to constrain.

The SPH simulations are run in the center-of-mass (COM) frame of the colliding moons, with the relative velocity at contact in the simulation box x-direction, as illustrated in Figure 1. The simulation results must therefore be rotated and translated into Saturn's frame to provide physically realistic pre-impact orbit conditions. The orbits of the resulting debris can then be computed, as detailed in Appendix B.

The pre-impact orbital elements of p-Rhea and p-Dione are, for the 2 km s⁻¹ collision,

$$\begin{eqnarray*}\begin{array}{rcl}{a}_{\mathrm{pR}} & = & 8.3\,{R}_{\unicode{x02644}},\,\,{e}_{\mathrm{pR}}=0.12,\,\,{i}_{\mathrm{pR}}=3^\circ ,\\ {a}_{\mathrm{pD}} & = & 6.09\,{R}_{\unicode{x02644}},\,\,{e}_{\mathrm{pD}}=0.2,\,\,{i}_{\mathrm{pD}}=5\buildrel{\circ}\over{.} 25,\end{array}\end{eqnarray*}$$

where a, e, and i are the semimajor axis, eccentricity, and inclination, respectively. The orbits intersect at the periapsis and ascending node of p-Rhea and the apoapsis and descending node of p-Dione. The relative velocity at this point of impact is

$$\begin{eqnarray*}&&{{\boldsymbol{v}}}_{\mathrm{rel}}=[1.544,0,1.273]\,\mathrm{km}\,{{\rm{s}}}^{-1},\end{eqnarray*}$$

and the COM has orbital elements

$$\begin{eqnarray*}&&{a}_{\mathrm{com}}=7.356\,{R}_{\unicode{x02644}},\,\,{e}_{\mathrm{com}}=0.0068,\,\,{i}_{\mathrm{com}}=0\buildrel{\circ}\over{.} 64.\end{eqnarray*}$$

For the 3 km s⁻¹ collision, the corresponding values are

$$\begin{eqnarray*}\begin{array}{rcl}{a}_{\mathrm{pR}} & = & 8.3\,{R}_{\unicode{x02644}},\,\,{e}_{\mathrm{pR}}=0.12,\,\,{i}_{\mathrm{pR}}=6^\circ ,\\ {a}_{\mathrm{pD}} & = & 6.09\,{R}_{\unicode{x02644}},\,\,{e}_{\mathrm{pD}}=0.2,\,\,{i}_{\mathrm{pD}}=10\buildrel{\circ}\over{.} 5,\\ {{\boldsymbol{v}}}_{\mathrm{rel}} & = & [1.608,0,2.538]\,\mathrm{km}\,{{\rm{s}}}^{-1},\\ {a}_{\mathrm{com}} & = & 7.246\,{R}_{\unicode{x02644}},\,\,{e}_{\mathrm{com}}=0.0059,\,\,{i}_{\mathrm{com}}=1\buildrel{\circ}\over{.} 30.\end{array}\end{eqnarray*}$$

The COM of the particles in the SPH simulation is at the origin at rest, and the relative velocity at impact is in the x-direction. The particles' positions and velocities are transformed into Saturn's rest frame by the following steps:

• 1.  
  Rotate about the z-axis such that v _∞ rather than v _t=0 is in the x-direction. This rotation is by <1° and has a negligible effect, but we include it for completeness.
• 2.  
  Rotate about the y-axis such that the direction of the relative velocity changes from x to that of v _rel. This corresponds to 395 and 586 for the two scenarios.
• 3.  
  Compute the position and velocity of the COM at the appropriate time from the input orbital elements, and then add them to the particles' positions and velocities.

As a sanity check, we confirm that these transformations are correct by analyzing the SPH initial snapshots to estimate the SPH moons' orbits, using the equations in Appendix B. These indeed reproduce the expected orbital elements for each moon, with slight variations caused by the shifted positions and velocities we imposed to yield the different impact angles. When the SPH simulation evolves further, the orbital estimates become less accurate, depending on how far material travels away from the COM. Using the final snapshots from our highly grazing impacts as examples, the semimajor axes remain close to the initial values, but by 9 hr after impact the extracted eccentricities can vary by ∼10%. This, alongside the less predictable effects on the debris evolution from tidal forces from Saturn, results in up to tens of percent uncertainty in the mass of debris that crosses different orbits. However, this remains a smaller potential variation than that between different impact angles and moon masses (Section 3.2); regardless, it would not affect the overall conclusions of this proof-of-concept study, since significant masses of debris can still be delivered throughout the system. Nevertheless, this simplification should be redressed in future work.

The surviving remnants and accumulated fragments produced by the SPH simulations are identified using a standard friends-of-friends algorithm, where particles within a certain distance of each other—the "linking length," l_link—are grouped together. At the relatively short time after impact of our analysis, most of these objects are still embedded in a debris field. This makes the precise choice of linking length somewhat arbitrary, because small changes in the linking length yield slightly different groups, in contrast with bodies in a vacuum with clear surfaces, which would be identified much the same for a wider range of linking lengths. For our primary resolution of 10^7.5 particles, we set l_link = 15 km, which scales for other resolutions with the cube root of the particle mass. We test the sensitivity of our results to this choice by repeating the analysis for significantly lower and higher l_link of ± 20% (12 and 18 km for 10^7.5 particles; see Figure 3). Even smaller or even larger linking lengths start to no longer include the ice mantles and rocky cores of bodies in the same group, or connect together many incohesive fragments into the same group, respectively.

For head-on to mid-angle impacts, a significant amount of ice and rock is ejected from both bodies, as illustrated for two example scenarios in Figure 2. The faster or closer to head-on the collision, the greater the disruption, the smaller the surviving remnants of the two original bodies (if any), and the further the debris is spread away from the impact point—both in and out of the orbital plane.

Zoom In Zoom Out Reset image size

Much of the ice is dispersed into diffuse debris, while the rock remains or more rapidly reaccretes into a large number of cohesive fragments. These disruptive impacts produce a relatively smooth mass distribution of many low-mass and fewer high-mass objects, as shown by Figure 3. Smaller and diffuse debris are not reliably resolved into discrete objects by the SPH simulations, so they are not included in the figure, and the detailed results below ∼10¹⁶ kg should be interpreted with caution. As such, although the total mass remains the same in all simulations, the cumulative distributions in Figure 3 show only the mass of resolved fragments. The corresponding results for 2 km s⁻¹ impacts are qualitatively similar, but with a more rapid trend to fewer large fragments at grazing angles (see Appendix A).

Zoom In Zoom Out Reset image size

For more grazing collisions, the pre-impact moons are incrementally less disrupted. Less ice and much less rock is ejected, particularly for impact angles at which the moons' cores do not intersect. The debris becomes dominated by the remnants of the two pre-impact bodies, with a smaller number of other fragments and diffuse ejecta (Figure 3). Although grazing collisions may not immediately deliver ring-forming or moon-disrupting material across the system, the orbits of the remnant moons are negligibly changed, ¹⁰ and they will likely collide again in a future orbit (Ć16), at perhaps more head-on angles.

Many fragments can be produced with masses above 10¹⁸ and even above 10¹⁹ kg (i.e., order 0.1–1 Mimas masses), which could readily disrupt other precursor moons analogous to Tethys, Enceladus, and Mimas on crossing orbits (Leinhardt & Stewart 2012), as discussed further in Section 3.2. The ice fractions of these larger fragments are typically a few tens of percent by mass, as more massive objects are better able to hold on to and reaccrete ice in addition to their rocky cores. Less ice is also ejected into distant debris by slower and/or more grazing impacts. For close to head-on collisions, fragments can be composed of up to equal mixes of p-Rhea and p-Dione material, while for mid-to-grazing angles, most large fragments are built from almost purely p-Rhea or p-Dione.

These results are similar for the different linking lengths discussed in Section 2.5, and any minor differences are negligible compared with those between impact angles and speeds, as illustrated in Figure 3. The ice fractions are similarly insensitive to the linking length, with the larger and smaller values typically yielding a few percent more and less ice by mass, respectively.

We now place the SPH results in the context of the Saturn system, following the procedure described in Section 2.4. Material is sent far throughout the system by head-on to mid-angle impacts, as shown in Figures 4 and 5. Large amounts of debris and fragments are placed onto orbits that would intersect with other satellites around the locations of those in the present-day system, including an extended population of debris that is directly heading deep into the Roche limit. This qualitative outcome is consistent across scenarios, with the angle and speed affecting, for example, the thickness and extent of the branches in Figure 4 (see Appendix A). We consider Saturn's Roche limit for water ice to be at 2.27 R_♄ = 136,800 km, which corresponds to the outer edge of the A ring. This makes it a somewhat conservative value, as, depending on the porosity and density of the debris, the effective limit could be higher (Tiscareno et al. 2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As discussed in Section 1, a precursor system can reasonably be expected to have multiple mid-sized satellites in a broadly similar architecture to the system today. The eccentric, high-velocity debris and fragments produced here could significantly erode other moons in a collisional cascade to distribute even more material throughout the system and into the Roche limit (Ć16), although this will have to be explored in detail with future modeling beyond the scope of this initial study. Head-on to mid-angle impacts send a few to a few tens of objects with masses above 10¹⁷ kg onto Enceladus-crossing orbits, for example, and an order of magnitude more for Tethys. Furthermore, a total mass of debris ranging from around that of Mimas to over that of Enceladus (order 10¹⁹–10²⁰ kg) reaches the orbits of the three inner moons.

The mass that crosses each location depends on the impact angle and speed, as shown in Figure 6. Table 1 contains the results for all simulations. The distributions of debris trajectories are highly anisotropic, so the orientation of the impact point in Saturn's frame (ϕ, Section 2.3) is also important and can significantly increase or decrease the crossing masses, depending on the impact angle. For these simulations, more than twice the present-day ring mass can be immediately sent inside the Roche limit, where it could begin to tidally and collisionally evolve toward forming rings (Dones et al. 2007), as discussed further in Section 3.4—even before considering further cascades and other material redistribution, or the wider parameter space of plausible impact scenarios beyond the few examples considered here.

Zoom In Zoom Out Reset image size

The Roche-entering material contains negligible rock, matching the present-day rings' nearly pure-ice composition (Zhang et al. 2017a, 2017b), while the material that could encounter Tethys and Enceladus can be several tens of percent rock by mass, depending on the scenario (Figure 6). The results for 2 km s⁻¹ scenarios are qualitatively similar, with the crossing masses reduced by factors of a few, approaching an order of magnitude at greater distances from the impact point, as detailed in Table 1.

In addition to the range of possible impact speeds, the applicability of a resonant destablization like that demonstrated by Ć16 to any pair of loosely Dione- and Rhea-like analogs reminds us that more mass could readily be delivered to the Roche limit by, for example, a more massive pair of colliding satellites (as examined in Section 3.3), or an impact at a distance closer to Saturn. The latter could arise from a different orbital evolution of the outer two satellites, or perhaps from a collision between p-Dione and a p-Tethys analog—or even from an outer moon destabilized as in Wisdom et al. (2022). Correspondingly, less material would likely be delivered by a more distant collision or a lower-mass pair, although this is less likely than higher precursor masses given that the present-day moons must accrete from the remnants.

A comparable ∼10¹⁹–10²⁰ kg mass of debris is also sent outward by disruptive impacts to encounter Titan. Like the material reaching the Roche limit, this is typically also pure ice, although some rock can be delivered to Titan by close to head-on collisions.

These initial trajectories for the debris are only the starting place for significant further evolution that must be examined and modeled in detail in future work. A common approach to estimating the evolution of collisional material is to assume that each piece of debris will evolve onto a circular orbit with the same angular momentum. However, disruptive collisions produce debris with a range of angular momenta and orbital energies, and detailed modeling beyond the scope of this paper will be required to constrain how much additional fragmentation and spreading of material occurs. If the system were isolated—though that is not the expectation here in a precursor satellite system—then the mass of material with periapses within the Roche limit and the mass of material with equivalent circular orbits within the Roche limit could be interpreted as crude upper- and lower-bound estimates for the mass likely to evolve there.

With these caveats in mind, Figure 7 shows the mass distribution of material with equivalent circular radius, ${a}_{\mathrm{eq}}=a{(1-{e}^{2})\cos (i)}^{2}$. A significant mass of pure-ice material can circularize inside the Roche limit even without subsequent disruption, though not immediately as much mass as that of the present-day rings for this example collision. The most dispersive scenarios tested here can yield up to ∼10¹⁹ kg of ice with a_eq < R_Roche, compared with ∼5 × 10¹⁹ kg with periapses <R_Roche. Therefore, under the conservative assumption of perfectly angular-momentum-retaining collisions and for the specific scenarios simulated here, some further redistribution of material via cascade collisions or other processes could be required to form the rings. Such evolution could be supported or even made unnecessary by a more dispersive base impact scenario, such as with larger or faster colliding moons, as discussed in Section 3.3.

Zoom In Zoom Out Reset image size

We also note that a comparable but perhaps greater challenge could be faced by hypotheses of ring formation from the tidal disruption of an external body or a destabilized moon (Dones 1991; Hyodo et al. 2017; Wisdom et al. 2022), with significantly greater initial semimajor axes and eccentricities. Such scenarios could still lead to the creation of rings but lack the range of highly reduced orbital energies and angular momenta of the initial debris that is produced by the dissipative collisions considered here.

Figure 7 also further demonstrates the lack of non-ice pollution in the initial ring-forming material, which for the other hypotheses could require a fortuitously ice-rich progenitor to match the observations (Dones 1991; Wisdom et al. 2022; Estrada & Durisen 2023). Here, no rocky material has an equivalent circular orbit within the Roche limit or even inside of Enceladus's orbit in this example, and none inside of Mimas's orbit for any of our impact simulations.

The total mass of debris that is placed onto orbits reaching the Roche limit and crossing the orbits of other moons numerically converges to within ∼1% by 10⁷ particles for disruptive scenarios, up to a few tens of percent at other impact angles, as shown in Figure 8 for 3 km s⁻¹ impacts, with complete results in Table 1. Low resolutions below ∼10⁶ particles yield unconverged masses that differ by many tens of percent for various impact angles. The masses for low-angle impacts converge at lower resolutions, apart from very close to head-on collisions, for which the mass crossing the Roche limit is lower and thus the resolution requirements for reliability are higher. This aligns with previous studies where the primary outcomes of high-speed, disruptive impacts converge relatively quickly, while impacts closer to the mutual escape speed can require much higher resolutions, as do extracted results of lower-mass features (Kegerreis et al. 2020). For 2 km s⁻¹ impacts, the mass crossing the Roche limit is lower, and thus convergence requires higher resolution, with most angles yielding differences of several percent between 10⁷ and 10^7.5 particles. In all cases, the ice fraction of the orbit-crossing debris converges to within a few percent by 10^6.5–10⁷ particles (Figure 8). Lower resolutions systematically underestimate the mass of rock delivered.

Zoom In Zoom Out Reset image size

For the mass function of debris, higher resolutions are important both for converging on the distribution of large objects and for enabling the inspection of smaller ones. We find that the overall shape of the high-mass end of the mass function is typically well converged by 10⁷ particles. The distribution of smaller objects with masses below ∼10¹⁸ kg is still broadly similar for ≥10⁷ particles but matches less reliably. The highest resolution used here enables the probing of objects down to masses of about 10¹⁶ kg, corresponding to ∼100 SPH particles, with potential numerical concerns for the analysis of smaller objects.

We also tested repeat simulations with reoriented initial conditions (Section 2.3), which with infinite resolution should give identical outcomes but can produce nonnegligibly different results even at high resolutions (Kegerreis et al. 2022). Fortunately, in these high-speed scenarios the differences are minimal for the orbit-crossing masses (standard deviations <1%), and the variation in the debris mass functions is of a similar magnitude to that between different linking lengths shown in Figure 3.

A potentially more significant limitation is the isolation of the SPH simulations from Saturn, as discussed in Section 2.4. Tidal effects would likely affect the separation and merging behavior of debris fragments, and the orbital estimates become less accurate. However, for the first several hours after the impact, the mass function is still rapidly evolving. By ∼8–9 hr (or sooner in highly grazing scenarios), the results remain broadly consistent with time. Some fragments continue to merge, collide, and separate, but this yields far less variation than found between different impact scenarios.

Replacing the impacting p-Dione with a larger p-Rhea produces similar results to the scenarios examined thus far. The mass functions of fragments follow the same trends shown in Figure 3 (and Figure 11), though as a symmetric collision there are two largest remnants rather than one, and at middle angles not quite as many large, secondary fragments are produced. The total mass that crosses the Roche limit is within a factor of two of the p-Dione–p-Rhea results in most cases but is over an order of magnitude higher for close to head-on impacts (see Table 1). The mass crossing the orbits of other present-day moons is similarly comparable to the p-Dione–p-Rhea cases but is typically somewhat higher for head-on to mid-angle impacts and lower at more grazing angles.

The above differences for p-Rhea–p-Rhea collisions are likely driven largely by the impact geometry, with some effects also expected from the slightly greater self-gravity of the moons and their post-impact remnants. At the same input impact angle, the centers of the two moons are more distant at the time of contact, so a smaller portion of the p-Rhea impactor intersects with the p-Rhea target than in the case of the smaller p-Dione impactor. This is compounded by Rhea's smaller rock fraction and lower density, such that grazing collisions can eject significant outer ice material but the rocky cores are less readily disrupted. This is reflected in the ice fractions of the orbit-crossing debris, which contain rock out to less grazing angles than in the p-Dione–p-Rhea case, as detailed in Table 1.

We also tested a few scenarios with systematically varying masses of the precursor moons, with the p-Rhea mass simplified to double that of the p-Dione. The overall outcomes are highly similar to the fiducial simulations at the same impact angle. The mass crossing the Roche limit increases with the colliding mass almost proportionally and roughly follows a power law with exponent ∼0.85–0.95 depending on the impact angle (Table 1). This substantiates the expectation that significantly more mass could be distributed throughout the system if the precursor moons were more massive than the present-day analogs considered for the primary simulations here.

In contrast to this general consistency in outcomes that we find across scenarios, our results show some differences from those of Hyodo & Charnoz (2017, hereafter HC17), who ran SPH impact simulations of two Rhea-mass bodies and found less spread of debris through the Saturn system, with a 45° impact delivering no material directly to the Roche limit, for example. They then took 10% of the particles from the 45° SPH simulation as inputs for five N-body integrations, with a different random subset of particles in each run. Further fragmentation was not included, but particles could merge, following a hard-sphere model, which they found leads most of the debris to reaccrete into a single body on a timescale of a few kiloyears.

Our SPH simulations differ in a few respects. One numerical difference is the SPH resolution used, as discussed above, which is over two orders of magnitude higher here. However, even with our test simulations at a comparable resolution to theirs, we still find more material on orbit-crossing trajectories for the same impact angle. Minor variations in the initial conditions preclude a direct comparison of the similar numerical techniques: First, for the impact point and the COM's orbit, HC17 set a circular, equatorial orbit with a 3 km s⁻¹ relative impact velocity. This places the impact point ∼13% farther away from Saturn than here. The predicted orbits of the debris are highly sensitive to the transformation of the SPH results into the Saturn system (see Section 2.4, underspecified in HC17 with respect to the orientation of the moons), and a more distant impact point reduces the inward reach of debris.

Second, HC17 collided Rhea-mass bodies with core mass fractions of 40% and 60%, compared with only 40% in both p-Rhea bodies here. This greater core fraction might also have reduced the amount of widespread ejecta. The interior structures (and masses) of the moons in a precursor system are unknown, so future work to systematically examine the outcomes from collisions of different initial bodies would help to constrain this uncertainty. Similarly, simulations that include Saturn and model impacts at various distances from the planet would reveal the importance both of the planet's gravity and of the location of the point of impact.

The conclusions that we draw from the SPH results in the context of the wider Saturn system also differ from HC17, as we speculate that collisions with other inner precursor moons could further support the delivery of material into the Roche limit, even with significantly less initial spread of debris than we find here. However, detailed N-body simulations that include both fragmentation and a full precursor satellite system will be required to make more robust predictions. It also may be that a more disruptive, less grazing impact than the 45° one examined by HC17 will be required for successful long-term distribution of material, even when fragmentation is accounted for.

HC17 also estimated analytically that the ∼kiloyear accretion timescale of their disk, neglecting fragmentation, would be shorter than the timescale of viscous spreading, which Ć16 had speculated could promote the delivery of material to inside the Roche limit. However, the immediate production of Roche-bound and Roche-circularizing debris that we find here and the possibilities of collisional cascades with other moons could avoid the requirement for a ring-forming disk to spread in this way. Furthermore, as also noted by Ć16, the first set of satellites that reaccrete from such a disk would likely be numerous and unstable. A longer timescale could be required for these satellites to merge into the relatively stable system we see today, during which the likely chaotic destablization and collisions could further distribute debris.

Both our primary simulations and those of HC17 used Tillotson (1962) EOSs (and neglected material strength), which do not include phase boundaries and can prompt unrealistic behavior for highly shocked ejecta (Stewart et al. 2020). A full exploration of the importance and effects of these systematic uncertainties is beyond the scope of this project, but we make an initial test of the sensitivity of our results using a repeated subset of the p-Rhea–p-Rhea impacts simulated with the more sophisticated ANEOS and AQUA EOS (Haldemann et al. 2020; Stewart et al. 2020). The overall behavior is similar to the Tillotson results, although in general the ice ejecta appears to be less diffusely dispersed. A greater proportion of ice also reaccretes into objects with masses around 10¹⁷–10¹⁸ kg, yielding a somewhat shallower size distribution and typical ice fractions in fragments ranging from ∼20% to 100%, compared with a previous ∼5%–40% for these p-Rhea–p-Rhea impacts. Most of the total masses of orbit-crossing debris are within a few tens of percent of the primary results, with some greater differences for the slower 2 km s⁻¹ collisions (Table 1). The variations due to the impact angle and speed continue to be more important than the choices of EOS, moon masses, compositions, and numerical resolution, though the differences may become more significant for future work on the detailed subsequent evolution of the debris.

The outputs of simulations like these provide the inputs for the next key stage of modeling the evolution of the initially eccentric and inclined fragments and debris. Combinations of N-body integrations and further impact simulations can then begin to determine the consequences for the other mid-sized moons and the possible growth or rejuvenation of circularized rings within Saturn's Roche limit.

Returning first to a brief discussion of the initial conditions for this scenario, the original hypothesis of recent reaccretion of the inner moons by Ć16 was based on the assumption of equilibrium tides. However, as discussed in Section 1, more sophisticated tidal models and the possibility of resonance locking (Lainey et al. 2020) do not avoid the recent encounter that an ancient Rhea would have had with the evection resonance. It is unlikely that this crossing could be reconciled with the relatively low-e, low-i orbit of Rhea without significant collisional damping, possibly even reaccretion (see Ć16), so an ancient age for Rhea remains highly unlikely. Additionally, a resonant tidal response offers many more possibilities for the system to be disrupted, as the energy that can be put into exciting satellite orbits by resonant tides is orders of magnitude larger than that provided by equilibrium tides. In future work, we hope to revisit the potential destablization of the system and expand on the input models that could lead to this general class of recent-collision scenarios, whether from evection excitation or, for example, the ejected outer satellite proposed by Wisdom et al. (2022).

There are also some limitations of this work that merit future study. While we explored a range of impact angles for two plausible speeds and brief tests of other moon masses, this scenario allows for a diverse variety of impact conditions, including different sizes and compositions of the colliding moons in addition to the angle, speed, and distance from Saturn of the impact. We also neglected material strength, which is unlikely to dominate in such high-speed impacts but could have uncertain effects on the details and thermodynamics of the ejecta—as would the use of more sophisticated EOSs than in the primary suite of simulations. Furthermore, the SPH simulations were performed in an isolated box before being placed around Saturn. The presence of Saturn throughout the hours of the simulation could promote the separation of fragments via tidal forces, for example, and its absence adds uncertainty to the values of the distributed mass. Therefore, while our primary conclusion that a significant amount of debris can be sent throughout the Saturn system to cross other moons and the Roche limit should not be sensitive to these limitations, the precise details of the debris and of the parameter space of successful specific scenarios may be expected to change as these simplifications are addressed.

The ejected fragments and debris will continue to interact with each other, with the other moons and any existing rings, and with Saturn's tidal field, before the system evolves toward its current state. However, once the present-day moons are mostly in place, any debris that is still in orbit—including secondary and further cascade ejecta—could imprint itself in the crater populations seen today. The craters on Saturn's moons indicate a population of at least some planeto- rather than heliocentric impactors (Ferguson et al. 2022) and do not match the distributions on the other outer planets' satellites. It is likely premature to consider the immediate fragment distributions from this study as a crater-forming population, since they can be expected to evolve significantly over time, but they can represent an initial approach toward future comparisons and constraints.

With this caveat in mind, Figure 9 shows the cumulative number distributions of debris objects from low to intermediate impact angles at an impact velocity of 3 km s⁻¹ (see also Figure 3). The distributions can be fitted crudely with a ∼1/M power law, with a slight suggestion of a trend to steeper slopes for higher impact angles, as comparatively fewer medium to large objects are produced by more grazing collisions. One can then begin to speculate on how this population of impactors could affect the reaccreting moons and what the impact rate on these bodies might be—again neglecting the fact that on these elliptical orbits significant erosion and fragmentation could first be expected. For example, a dense, planetocentric population may provide a sufficient "fluence" of objects to saturate satellite surfaces in a relatively short period of time, perhaps ∼kiloyears, so that they may appear as cratered as ancient, gigayear-old terrain might be from a heliocentric population (Lissauer et al. 1988; Zahnle et al. 2003). Furthermore, the range of impact speeds will also be lower for planetocentric orbits, so different scaling laws for impactor size to crater diameter may be required. The thermal state of the evolving moons should also be considered, in terms of the cooling time that might be required before craters can be preserved without significant relaxation (White et al. 2017). However, the low gravitational binding energy of the small moons means that large impactors may cause more erosion than heating.

Zoom In Zoom Out Reset image size

In addition to the disruption, accretion, and eventual surfacing implications for Saturn's inner mid-sized moons, approximately an Enceladus mass of debris can be placed onto orbits that intersect with Titan. Titan's unique atmosphere and surface present many unanswered questions regarding their evolution and history (Nixon et al. 2018). This includes a possible young age for the current methane system of a few tens to hundreds of megayears associated with the rapid conversion of methane into heavier hydrocarbons. The effects of this high-energy delivery of a large mass of icy debris are beyond the focus of this study but could open up new options for models of Titan's evolution. Furthermore, Ć16 also proposed that resonant interactions of Titan with potentially long-lived outer parts of the debris disk and a transient satellite that accretes there could excite Titan's eccentricity enough to explain its significant value (∼0.03) today.

Finally, we could consider this general scenario for the destablization and disruption of a satellite system as a potential evolutionary pathway for other planets' moons and rings. However, among the other giant planet satellite systems in our own solar system, a similar evection resonance destablization would likely not have occurred recently, if at all. In the case of Jupiter, where the evection resonance lies between Io and Europa, ¹¹ the resonant configuration of the inner three Galilean moons has long been known to be stable (Laplace et al. 1829) and is likely ancient, perhaps as old as the solar system (Goldreich & Soter 1966). At Uranus, evection lies between Ariel and Umbriel, neither of which is likely to have crossed owing to the planet's large tidal Q (Ćuk et al. 2020). Lastly, whether such a destablization could have happened at Neptune is unknown, since any preexisting satellite system can be expected to have been disrupted owing to the capture of Triton (Goldreich et al. 1989; Agnor & Hamilton 2006).