A recent impact origin of Saturn's rings and mid-sized moons

We simulate the collision of precursor icy moons analogous to Dione and Rhea as a possible origin for Saturn's remarkably young rings. Such an event could have been triggered a few hundred million years ago by resonant instabilities in a previous satellite system. Using high-resolution smoothed particle hydrodynamics simulations, we find that this kind of impact can produce a wide distribution of massive objects and scatter material throughout the system. This includes the direct placement of pure-ice ejecta onto orbits that enter Saturn's Roche limit, which could form or rejuvenate rings. In addition, fragments and debris of rock and ice totalling more than the mass of Enceladus can be placed onto highly eccentric orbits that would intersect with any precursor moons orbiting in the vicinity of Mimas, Enceladus, or Tethys. This could prompt further disruption and facilitate a collisional cascade to distribute more debris for potential ring formation, the re-formation of the present-day moons, and evolution into an eventual cratering population of planeto-centric impactors.


INTRODUCTION
Whatever their origin, Saturn's rings appear to be young. Key observations made during the Cassini mission provided new measurements of the ring mass (Iess et al. 2019), the fraction of non-icy material in the rings (Zhang et al. 2017b,a), and the extrinsic micrometeoroid flux at Saturn (Kempf et al. 2023). Together, these three factors constrain the ring age to be less than a few 100 Myr (Estrada et al. 2018;Kempf et al. 2023;. This is because the primarily icy rings (>95% by mass, Doyle et al. 1989;Cuzzi & Estrada 1998;Zhang et al. 2017b,a) are continuously subjected to micrometeoroid bombardment and are darkened over time, so an upper bound on the rings' exposure age can be determined by assuming that they began as pure ice (Cuzzi & Estrada 1998;Estrada et al. 2015). The present mass of the main rings is ∼0.4 Mimas masses (∼1.5 × 10 19 kg).
Moreover, the rings appear to be losing mass at a remarkable rate (Waite et al. 2018;Hsu et al. 2018;O'Donoghue et al. 2019), suggesting they are not only young, but ephemeral as well.  showed that the process of direct deposition of micrometeoroids and the ballistic transport of their impact ejecta can account for these mass-loss rates, and  used these results to demonstrate that mass loading and ballistic transport due to micrometeoroid bombardment continue to drive the rings' dynamical evolution once viscosity becomes too weak. These studies concluded that if the rings have been losing mass at more or less their current rate, then they may have begun with a mass of one to a few times the mass of Mimas a few 100 Myr ago. They also find that, when ballistic transport is included, the rings' dynamical age and exposure age are similar -contrary to the view that these timescales can be vastly different, which would otherwise allow for the rings to look young but be dynamically old (Crida et al. 2019).
The rings' geological youth and apparent transience are at odds with most current ideas about their origin, in which they are either primordial (Pollack et al. 1976;Canup 2010), or derived from the remnants of either a collisionally disrupted, Mimas-mass moon (Harris 1984;Charnoz et al. 2009;Dubinski 2019) or a tidally split comet (Dones 1991;Dones et al. 2007;) -catastrophes that are highly unlikely in the current solar system and hence presumed ancient. In the latter tidal disruption case, a comet or Centaur a few hundreds of kilometres in size is required to pass close enough to Saturn to be tidally disrupted, with the rings forming from some of the debris. The collisional disruption of a Mimas-mass moon requires a body tens of kilometres in size. Furthermore, the precursor moon must be kept close enough to Saturn long enough for the impact event to occur, although this could be facilitated by mean-motion resonances with Enceladus and Dione (Dubinski 2019). However, the flux of such bodies decreases sharply after the Late Heavy Bombardment (Zahnle et al. 2003), giving only a probability of occurrence of ∼10 −4 in the past few 100 Myr (Ip 1988;Lissauer et al. 1988;Dones 1991;Dones et al. 2007;Charnoz et al. 2009). Wisdom et al. (2022) recently explored a novel scenario for recent ring formation in which Saturn once had an additional moon between Titan and Iapetus. This Iapetus-mass satellite helped to maintain a secular resonance with Neptune that drove up Saturn's obliquity. Subsequent destabilisation of the moon's orbit by resonances with Titan could have led to a grazing encounter with Saturn and a recent ring origin ∼100 Myr ago, in addition to kicking Saturn out of the resonance. This would imply that Saturn is no longer in the resonance with Neptune, but confirming this may require more accurate measurements of Saturn's pole precession rate (Jacobson 2022). Detailed predictions are also yet to be made for the fate of the destabilised moon's rock and ice material, if it were sent on an elliptical orbit that crosses into the Roche limit (Dones 1991), or for its potentially disruptive effects on the rest of the satellite systemperhaps along similar lines to the collisional scenario we consider below.
In addition to the rings, Saturn's current system of icy moons (interior to Titan) also presents some evidence of recent formation. Ćuk et al. (2016, hereafterĆ16) demonstrated that if the present-day mid-sized satellite system were ancient, then recent tidal evolution through dynamical resonances would have excited the moons' orbits beyond what is seen today. In particular, the "smoking gun" for the moons not being ancient is Rhea's present inclination, which is an order of magnitude lower today than if it had survived crossing an evection resonance 1 at around 8.3 R , where R = 60, 268 km is the equatorial radius of Saturn (at 1 bar, Lindal et al. 1985), and such inclinations are not readily damped (Murray & Dermott 1999;Chen et al. 2014). Rhea today is located outside evection, at 8.7 R , with an inclination of only 0.35 • , so it cannot have formed long ago and migrated outwards. A very slow tidal evolution rate might also have been able to explain the present-day moons' unexcited orbits, but this would conflict with the apparent fast orbital evolution of Rhea (Lainey et al. 2020;Jacobson 2022) and the intense tidal heating of Enceladus (Meyer & Wisdom 2007). The alternative is that Rhea formed (or re-formed) outside the resonance (Ć16). This implies a similarly recent, and thus perhaps shared, formation to the rings for at least some of these icy satellites.
More recent work has shown that at least some moons of Saturn evolve much faster than predicted by the equilibrium tides assumed byĆ16 (Lainey et al. 2020;Jacobson 2022), possibly indicating resonance locking of the moons with gravity modes or inertial waves within Saturn (Fuller et al. 2016). The apparent presence of resonant tides weakens the constraints on the age of the system from the lack of a past Tethys-Dione 3:2 mean-motion resonance found byĆ16, as it is no longer clear that these moons are currently on converging orbits. However, the observed fast migration of Rhea, if driven by resonance locking (Lainey et al. 2020), still implies that an ancient Rhea would have crossed and had its orbit excited by the solar evection resonance only a few hundred Myr ago.
Here we examine a new scenario based on the work ofĆ16, which hypothesizes a recent dynamical instability in a precursor satellite system ∼100 Myr ago. Any pair of moons with orbits and masses loosely comparable with present-day Dione and Rhea will be suddenly excited when the Rhea-analogue outer moon encounters the evection resonance during migration. Both moons have their orbits excited to significant eccentricities (∼0.1-0.2) and inclinations of several degrees. This paired excitation can naturally arise from the coupled migration of moons in mean-motion resonance; from the inner moon catching up owing to the stronger tidal effects nearer Saturn; and from an initially separate outer moon being trapped just interior to evection for a long time, by a cycle of temporary inwards migration caused by the damping of evection-driven eccentricity (Ć16).
Regardless of the precise evolution to this point, the common outcome of large eccentricities for both moons leads to orbit crossing and high-velocity collisions. Fragments and debris ejected by this impact could then reaccrete into the present-day moons, significantly erode other precursor moons leading to a collisional cascade, and perhaps even deliver material directly inside Saturn's Roche limit to form or rejuvenate rings.
This scenario does not rule out previous hypotheses for ancient satellite formation around Saturn, such as by accretion from primordial rings (Crida & Charnoz 2012. Such mechanisms would now simply apply to the formation of the precursor satellite system that was then recently destabilised and disrupted to (re-)form the moons we see today.
Some aspects of theĆ16 scenario were investigated using smoothed particle hydrodynamics (SPH) and Nbody simulations of two colliding Rhea-mass moons , with results that suggested that the debris would re-accrete into a satellite without forming rings. However, this conclusion could be premature for a few reasons, including: only a single impact angle of 45 • was considered as an example input for the longerterm models; the SPH simulations used a low resolution of 2 × 10 5 particles; the subsequent evolution assumed a simplified hard-sphere and dissipative coefficient-ofrestitution model for the re-accretion of fragments,without the possibility of further fragmentation; and the presence of other precursor moons in the system and potential cascade interactions with them were neglected.
Regarding the last issue of other satellites, theĆ16 scenario is based on a multiple-moon primordial system comparable to the current one. Such a precursor satellite system would likely have a broadly similar mass and architecture to the present system -perhaps with a somewhat larger total mass. The current systems of both Saturn and Uranus (and probably at one time, Neptune; Goldreich et al. 1989) have families of mid-sized moons that extend inward close to their parent planets, so this expectation appears well justified. Furthermore, the initial distribution of moons cannot be random, since their radial locations must satisfy orbital stability criteria (e.g. Wisdom 1980;Lissauer 1995) that restrict the regions in which moons can orbit stably (e.g. Fig. 1b, Mosqueira & Estrada 2003a). We also can expect size sorting (increasing radius with distance), as it is a natural outcome of models that form mid-sized moons in a gaseous circumplanetary disk (Mosqueira & Estrada 2003b) or if the moons were "spun out" from viscously evolving, massive, primordial rings after planet formation (Crida & Charnoz 2012, while a collisional origin might allow some variation in size and composi-tion (Asphaug & Reufer 2013). We therefore take the current Saturn system as a guide for the context and implications of impacts in a precursor one.
Here we use SPH simulations with over two orders of magnitude higher resolution than  to model collisions based onĆ16 integrations. We find that such an impact can send significant debris throughout the system, including interior to Saturn's Roche limit, and can produce an extended distribution of massive fragments. These objects are likely to continue to fragment or erode and circularise with one another on future orbits, and they also have semi-major axes and eccentricities that will lead to crossing orbits, close encounters, and possible collisions with other precursor moons. This opens the possibility of a collisional cascade that leads to additional rubble piles and debris crossing within the Roche limit, which may be tidally disrupted or collisionally captured (if some material is already present) into rings.

Smoothed particle hydrodynamics simulations
Smoothed particle hydrodynamics (SPH, Lucy 1977;Gingold & Monaghan 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 modelling 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 5 -10 6 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. 2019Kegerreis et al. , 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 2 (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 Sandnes et al. 2023). As in previous works, we also neglect material strength, as gravitational stresses are expected to domi-nate (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).

Initial conditions
To distinguish the Dione-/Rhea-like moons we collide from their present-day analogues, 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 and 2.31 × 10 21 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 modelled with Tillotson (1962) basalt and water EoS (Melosh 2007), or with ANEOS forsterite and AQUA water (Stewart et al. 2020;Haldemann et al. 2020) for a small subset of comparison simulations.
We generate the moons' internal profiles by integrating inwards while maintaining hydrostatic equilibrium 3 , then place the roughly equal-mass SPH particles to match the resulting density profiles using the stretched equal-area method 4 . Before simulating the impact, we first run a brief 2 h 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 h before contact to 9 h after. As this is only a modest fraction of the orbital period around Saturn (∼82 h), we neglect the effects of Saturn's gravity during the SPH simulations, as discussed in §2.4.

Impact scenarios
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 −1 . The orbital parameters are listed in §2.4. For comparison, the mutual escape speed is v esc ≈ 580 m s −1 , so the impacts are at ∼3.4 and 5.1 v 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 −1 ) 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 semi-major 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 Fig. 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 centre 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 • , 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 21 kg is 10 5 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 , although they used non-equal 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, and 8×10 21 kg and corresponding p-Rhea masses of double these values. These moons are given a core mass fraction of 1 3 . We simulate these collisions for β = 15 • , 25 • , 35 • for the four mass pairs. Finally, to further examine the numerical effects of our finite-particle models, each scenario at our secondhighest resolution of 10 7 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 non-negligible effects even at the high resolutions used here (Kegerreis et al. 2020(Kegerreis et al. , 2022 -an uncertainty that these tests will help to constrain.

Transforming SPH outputs into Saturn's frame
The SPH simulations are run in the centre-of-mass (CoM) frame of the colliding moons, with the relative SPH simulation Figure 1. The pre-impact orbits of p-Rhea, p-Dione, and their centre of mass (CoM) around Saturn in dark red, blue, and green, respectively, in x-y and x-z projections (where z = 0 is Saturn's equatorial plane), for the 3 km s −1 scenario detailed in §2.4. The grey circle shows Saturn, to scale. The inset zoom-in panels show, as labelled, the moons and their velocities shortly before impact in Saturn's frame and in the centre of mass frame. The bottom-right panel shows the orientation for the SPH simulation, with the velocity at contact in the x direction and an impact angle β (ignoring any tidal deformation), plus a possible extra rotation about the x axis in this frame by an angle ϕ, to change the orientation of the collision once placed in Saturn's frame. The initial separation for the simulation is set such that the time to impact is 1 hour.
velocity at contact in the simulation-box x direction, as illustrated in Fig. 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 Appx. B.
The pre-impact orbital elements of p-Rhea and p-Dione are, for the 2 km s −1 collision: a pR = 8.3 R , e pR = 0.12, i pR = 3 • , a pD = 6.09 R , e pD = 0.2, i pD = 5.25 • , where a, e, and i are the semi-major 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: The centre of mass 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 39.5 • and 58.6 • for the two scenarios.
3. Compute the position and velocity of the centre of mass at the appropriate time from the input orbital elements, then add them to the particles' positions and velocities.
As a sanity check, we confirm that these transformations are correct by analysing the SPH initial snapshots to estimate the SPH moons' orbits, using the equations in Appx. 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 Illustrative final snapshots in the SPH box frame from simulations with 10 7.5 particles for impact angles and speeds at infinity of (left) β = 15 • , v∞ = 2.0 km s −1 and (right) β = 35 • , v∞ = 3.0 km s −1 . The lower and upper panels show the x-y and x-z planes, respectively. Orange and blue show p-Rhea's core-rock and mantle-ice material, respectively, and yellow and purple the same for p-Dione. Animations are available at icc.dur.ac.uk/giant impacts/icy moons 2kms py.mp4 and ... 3kms py.mp4, and with the same data rendered in 3D at icc.dur.ac.uk/giant impacts/icy moons 2kms.mp4 and ... 3kms.mp4.
from the centre of mass. Using the final snapshots from our highly grazing impacts as examples, the semi-major axes remain close to the initial values, but by 9 h 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 ( §3.2), and regardless would not affect the overall conclusions of this proof-ofconcept study, since significant masses of debris can still be delivered throughout the system. Nevertheless, this simplification should be redressed in future work.

Identifying debris objects
The surviving remnants and accumulated fragments produced by the SPH simulations are identified using a standard friends-of-friends (FoF) 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 Fig. 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.

RESULTS AND DISCUSSION
We first examine the general outcomes of collisions between precursor moons in the center-of-mass frame in which the SPH simulations were run in §3.1, before then considering the context of the full Saturn system in §3.2 (see Fig. 1). We examine numerical and other-scenario comparisons for the primary-simulation results in §3.3,  Figure 3. The cumulative mass distributions of objects produced by impacts with v∞ = 3 km s −1 , after 9 h, simulated with 10 7.5 SPH particles, coloured by the impact angle. Additional material is ejected as diffuse debris or smaller, unresolved objects that bring the cumulative total mass in all cases to the same value beyond the low end of the horizontal axis shown. The faint lines above and below each primary result show the distributions for 20% smaller and larger linking lengths.
with a discussion of the limitations and speculative implications of the results in §3.4.

Impact outcomes
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 Fig. 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.
Much of the ice is dispersed into diffuse debris, while the rock remains or more rapidly re-accretes 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 Fig. 3. Smaller and diffuse debris are not reliably resolved into discrete objects by the SPH simulations, so are not included in the figure, and the detailed results below ∼10 16 kg should be interpreted with caution. As such, although the total mass remains the same in all simulations, the cumulative distributions in Fig. 3 show only the mass of resolved fragments. The corresponding results for 2 km s −1 impacts are qualitatively similar, but with a more rapid trend to fewer large fragments at grazing angles (see Appx. A).
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 (Fig. 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 5 , 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 18 and even above 10 19 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 §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 re-accrete 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 §2.5, and any minor differences are negligible compared with those between impact angles and speeds, as illustrated in Fig. 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.

Trajectories in the Saturn system
We now place the SPH results in the context of the Saturn system, following the procedure described in §2.4. Material is sent far throughout the system by head-on to mid-angle impacts, as shown in Figs. 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 Fig. 4 (see Appx. A). . Each point is one SPH particle with a mass of ∼10 14 kg. The dotted lines and corresponding particle colours indicate where an orbit's periapsis crosses the Roche limit or the orbits of other moons at their locations in the present-day system, or the apoapsis in the case of Titan. Other particles with semi-major axes inside and outside the impact point are coloured in grey and light grey, respectively. The dark-red and dark-blue pairs of axis ticks indicate the orbital parameters of the pre-impact moons.
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). As discussed in §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 modelling 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 17 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 19 -10 20 kg) reaches the orbits of the three inner moons.  Table 1. The plotted values are exclusive, so e.g. the mass shown for Enceladus's orbit does not include the material that also reaches Mimas or the Roche limit.
The mass that crosses each location depends on the impact angle and speed, as shown in Fig. 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 (ϕ, §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 towards forming rings (Dones et al. 2007), as discussed further in §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.
The Roche-entering material contains negligible rock, matching the present-day rings' nearly pure-ice composition (Zhang et al. 2017a,b), while the material that could encounter Tethys and Enceladus can be several tens of percent rock by mass, depending on the scenario (Fig. 6). The results for 2 km s −1 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 destabilisation like that demonstrated byĆ16 to any pair of loosely Dione-and Rhea-like analogues 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 §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 analogueor even from an outer moon destabilised 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 19 -10 20 kg mass of debris is also sent outwards 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 modelled 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 mo- mentum. However, disruptive collisions produce debris with a range of angular momenta and orbital energies, and detailed modelling 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 upperand lower-bound estimates for the mass likely to evolve there.
With these caveats in mind, Fig. 7 shows the mass distribution of material with equivalent circular radius, a eq = a(1 − e 2 ) cos(i) 2 . A significant mass of pure-ice material can circularise 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 19 kg of ice with a eq < R Roche , compared with ∼5 × 10 19 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 evo-lution 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 §3.3.
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 destabilised moon (Dones 1991;Wisdom et al. 2022), with significantly greater initial semi-major 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. Fig. 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 icerich progenitor to match the observations (Dones 1991;Wisdom et al. 2022;. 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.

Convergence and comparisons
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 7 particles for disruptive scenarios, up to a few tens of percent at other impact angles, as shown in Fig. 8 for 3 km s −1 impacts, with complete results in Table 1. Low resolutions below ∼10 6 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 −1 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 7 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 7 particles (Fig. 8). Lower resolutions systematically underestimate the mass of rock delivered.
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 7 particles. The distribution of smaller objects with masses below ∼10 18 kg is still broadly similar for ≥10 7 particles, but matches less reliably. The highest resolution used here enables the probing of objects down to masses of about 10 16 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 ( §2.3), which with infinite resolution should give identical outcomes, but can produce non-negligibly different results even at high resolutions (Kegerreis et al. 2022). Fortunately, in these highspeed scenarios the differences are minimal for the orbitcrossing 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 Fig. 3.
A potentially more significant limitation is the isolation of the SPH simulations from Saturn, as discussed in §2.4. Tidal effects would likely affect the separation and merging behaviour 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 h (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 Fig. 3 (and Fig. 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 presentday 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 centres 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 analogues 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 Rheamass bodies and found less spread of debris through the Saturn system, with a 45 • impact delivering no mate-rial 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 re-accrete into a single body on a timescale of a few kyr.
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 centre-ofmass's orbit, HC17 set a circular, equatorial orbit with a 3 km s −1 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 §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 longterm distribution of material, even when fragmentation is accounted for.
HC17 also estimated analytically that the ∼kyr 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 Rochecircularising 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 re-accrete 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 destabilisation and collisions could further distribute debris.
Both our primary simulations and those of HC17 used Tillotson (1962) equations of state (and neglected material strength), which do not include phase boundaries and can prompt unrealistic behaviour 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 (Stewart et al. 2020;Haldemann et al. 2020). The overall behaviour 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 17 -10 18 kg, yielding a somewhat shallower size distribution and typical ice fractions in fragments ranging from ∼20-100%, compared with a previous ∼5-40% for these p-Rheap-Rhea impacts. Most of the total masses of orbitcrossing debris are within a few tens of percent of the primary results, with some greater differences for the slower 2 km s −1 collisions ( Table 1). The variations due to the impact angle and speed continue to be more important than the choices of equation of state, moon masses, compositions, and numerical resolution, though the differences may become more significant for future work on the detailed subsequent evolution of the debris.

Implications and future work
The outputs of simulations like these provide the inputs for the next key stage of modelling 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 circularised rings within Saturn's Roche limit.
Returning first to a brief discussion of the initial conditions for this scenario, the original hypothesis of recent re-accretion of the inner moons byĆ16 was based on the assumption of equilibrium tides. However, as discussed in §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 re-accretion (cf.Ć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 destabilisation 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 equations of state 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 towards its current state. However, once the present-day moons are mostly in place, any debris  Figure 9. The cumulative number distributions of debris objects produced by head-on to mid-angle 3 km s −1 impacts, as also characterised by Fig. 3. The lines show fitted power laws N ∝ M b , with exponents b given in the legend, for the mass range 10 16 -2 × 10 18 kg. The fits do not include the many smaller objects that are less robustly resolved, nor the lower numbers of even larger objects that are less well represented by a single power law in some cases.
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 helio-centric 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 craterforming population, since they can be expected to evolve significantly over time, but they can represent an initial approach towards future comparisons and constraints.
With this caveat in mind, Fig. 9 shows the cumulative number distributions of debris objects from low-to-mid impact angles at an impact velocity of 3 km s −1 (see also Fig. 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 how this population of impactors could affect the re-accreting 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, planeto-centric population may provide a sufficient "flu-ence" of objects to saturate satellite surfaces in a relatively short period of time, perhaps ∼kyr, so that they may appear as cratered as ancient, Gyr-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 planeto-centric 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.
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 Myr 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, C16 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 destabilisation and disruption of a satellite system as a potential evolutionary pathway for other planets' moons and rings. However, amongst the other giant planet satellite systems in our own solar system, a similar evection-resonance destabilisation would likely not have occurred recently, if at all. In the case of Jupiter, where the evection resonance lies between Io and Europa 6 , 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 are likely to have crossed owing to the planet's large tidal Q (Ćuk et al. 2020). Lastly, whether such a destabilisation could have happened at Neptune is unknown, since any pre-existing satellite system can be expected to have been disrupted owing to to the capture of Triton (Goldreich et al. 1989;Agnor & Hamilton 2006).

CONCLUSIONS
Dynamical instabilities in a precursor system of Saturnian mid-sized icy moons can lead to high-velocity collisions that scatter fragments and debris throughout the system (Ćuk et al. 2016). Using SPH simulations with over two orders of magnitude higher resolution than previous studies, we find that a range of plausible impacts can deliver significant mass directly inside the Roche limit, with a ring-like composition of pure ice. Furthermore, more than a Mimas mass of material -and even more than an Enceladus mass in some cases -is placed onto crossing orbits with present-day Mimas, Enceladus, and Tethys (and Titan), facilitating the possibility of a collisional cascade to further distribute material across the system.
Head-on to mid-angle collisions produce large numbers of cohesive fragments and re-accreted objects with masses above 10 18 -10 19 kg (∼0.1-1 Mimas masses). Most of these are -at 9 hours after impact -primarily composed of rock in a more diffuse field of ejected ice. Grazing impacts produce little debris but the colliding moons' orbits are also negligibly changed, so the scenario retains the opportunity of a more disruptive collision on a future orbit as the satellites continue to cross paths. Slower collisions are less effective at distributing material, but they can still send small amounts of debris across the Roche limit, and send significantly more that would intersect with other precursor moons to potentially cause further disruption. Impacts between two equal-sized moons yield similar results and trends.
Even before considering subsequent collisions and future modelling of the further creation and redistribution of debris, we find that the example impacts considered here can place over 5×10 19 kg of ice onto Roche-crossing orbits (and no rock), and ∼10 18 -10 19 kg of pure ice with equivalent circular orbits inside the Roche limit. However, that these masses are factors of a few above and below the mass of the present-day rings is less important for this proof-of-concept study than the qualitative placement of significant masses onto Roche-and moon-crossing orbits across a range of impact scenarios -to which the production of immediately Rochecircularising material is an encouraging bonus. Furthermore, the mass on Roche-crossing orbits scales close to proportionally with the mass of the colliding moons, so the larger moons that might be expected in a precursor system could more readily produce massive rings.
The mass of debris entering the Roche limit converges for moderately high numerical resolutions in most cases, while the mass distribution requires at least ∼10 7 SPH particles to converge reliably on the population of midto-large fragments, let alone smaller objects.
We conclude that the impact of two destabilised icy moons is a promising scenario for the recent formation or rejuvenation of Saturn's rings and re-accretion of mid-sized moons. Future work on the long-term evolution of the orbit-crossing debris, combined with further and more detailed modelling of collisions between both icy moons and smaller fragments, will help to constrain the implications of this scenario for Saturn's rings, its moons, their craters, and other surface environments. ACKNOWLEDGMENTS J.A.K. acknowledges support from a NASA Postdoctoral Program Fellowship, administered by Oak Ridge Associated Universities. J.N.C. was supported for part of this work by the Cassini project, and partly by a NASA Cassini Data Analysis Program grant to P.R.E. We thank K. Zahnle for valuable input and discussion. We also thank A. Rhoden for an early look at a pending chapter on cratering evidence for icy moon evolution. L.F.A.T. would like to dedicate this article to the memory of Michael S. Warren. We thank the anonymous reviewer for their helpful comments. The research in this paper made use of the SWIFT open-source simulation code (Schaller et al. 2018(Schaller et al. , 2023 Software: SWIFT (www.swiftsim.com, Schaller et al. 2018.0); WoMa (pypi.org/project/woma/, Ruiz-Bonilla et al. 2021). For comparison with the head-on to mid-angle impacts in Fig. 2, example outcomes for more grazing collisions are illustrated in Fig. 10. Less debris is ejected, and the mostly intact remnants of the two moons continue on little-changed trajectories so are likely to re-collide on a future orbit, as discussed in the main text.
The mass distributions from slower impacts are shown in Fig. 11, for comparison with the faster ones in Fig. 3. The general outcomes and trends are similar, but the debris become dominated by the two remnants of the initial moons at a less grazing impact angle.
The resulting Saturn-frame orbits of the debris from the slower, more head-on scenario illustrated in Fig. 2(left) are shown in Figs. 12 and 13. The overall results are similar to the faster, higher-angle example of Figs. 4 and 5. The debris is more spread out around the impact point, which is representative for more head-on collisions in general, with somewhat less material reaching the Roche limit, which is representative for lower speeds in general. The inclinations are also significantly less high, because the initial moons' material has been more thoroughly mixed, and the debris therefore lies on orbits closer to the near-equatorial one of the centre of mass.
B. ORBITAL ELEMENTS After the SPH particle data are transformed into Saturn's frame, as described in §2.4, we estimate the orbit of each particle or fragment from its position, ⃗ r, and velocity, ⃗ v. The vis-viva equation yields the semi-major axis, a: where r = |⃗ r|, v = |⃗ v|, and µ = G(M + m), with m the particle or fragment's mass and G the gravitational constant.
The eccentricity, e, and inclination, i, are then derived from the specific angular momentum, ⃗ h = ⃗ r × ⃗ v: The longitude of ascending node, Ω, is then where the sign of Ω is flipped if h x < 0.

C. RESULTS TABLES
The primary results from all simulations are listed in Table 1, arranged in groups by subset, impact angle, and speed. The SPH simulation data will be made available on reasonable request. Table 1 is available in the online version in a machine-readable format. Table 1. Outcomes of icy-moon impact simulations. Columns list the impact angle, β; speed at infinity, v∞; the mass of debris on orbits crossing the Roche limit and present-day moons; and the fraction of which is ice. For the base suite of simulations, each entry is followed by three more lines with the SPH frame rotated about the x axis by ϕ = 90 • , 180 • , 270 • (see Fig.1). The 'type' notes other parameters that are changed from the base scenario for the following subsets of simulations, as described in §2.3, within which the angle and speed are then varied, namely: the number of particles, N (base 10 7.5 ); the body colliding with p-Rhea (base p-Dione); the EoS (base Tillotson), and the mass of p-Dione for the separate tests of moon masses, MpD in units of 10 21 kg, with corresponding p-Rhea mass MpR = 2 MpD. Note: the orbit-crossing masses include only the mass that reaches that orbit, and exclude material that also reaches further ones (e.g. the Enceladus values do not include Mimas-or Roche-crossing material).