# The peculiar shapes of Saturn's small inner moons as evidence of mergers of similar-sized moonlets

A. Leleu, M. Jutzi, M. Rubin

► To cite this version:

A. Leleu, M. Jutzi, M. Rubin. The peculiar shapes of Saturn's small inner moons as evidence of mergers of similar-sized moonlets. Nature Astronomy, 2018, 2 (7), pp.555-561. 10.1038/s41550-018-0471-7. hal-02316323

HAL Id: hal-02316323

https://hal.science/hal-02316323v1

Submitted on 8 Oct 2025

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

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

<figure><img src="image_1.png" /></figure>
# The peculiar shapes of Saturn's small inner moons as evidence of mergers of similar-sized moonlets

A. Leleu<sup>1,2*</sup>, M. Jutzi<sup>1</sup> and M. Rubin<sup>1</sup>

The Cassini spacecraft revealed the spectacular, highly irregular shapes of the small inner moons of Saturn<fnref n="1" />, ranging from the unique 'ravioli-like' forms of Pan and Atlas<fnref n="2" />,³ to the highly elongated structure of Prometheus. Closest to Saturn, these bodies provide important clues regarding the formation process of small moons in close orbits around their host planet⁴, but their range of irregular shapes has not been explained yet. Here, we show that the spectrum of shapes among Saturn's small moons is a natural outcome of merging collisions among similar-sized moonlets possessing physical properties and orbits that are consistent with those of the current moons. A significant fraction of such merging collisions take place either at the first encounter or after 1-2 hit-and-run events, with impact velocities in the range of 1-5 times the mutual escape velocity. Close to head-on mergers result in flattened objects with large equatorial ridges, as observed on Atlas and Pan. With slightly more oblique impact angles, collisions lead to elongated, Prometheus-like shapes. These results suggest that the current forms of the small moons provide direct evidence of the processes at the final stages of their formation, involving pairwise encounters of moonlets of comparable size⁴⁻⁶. Finally, we show that this mechanism may also explain the formation of lapetus' equatorial ridge⁷, as well as its oblate shape⁸.

The small inner moons Atlas, Prometheus, Pandora, Janus and Epimetheus are repelled by the rings of Saturn at a rate that is proportional to their mass and decreases with their distance⁴,⁹. It has thus been proposed that these moons were formed in a pyramidal regime (that is, by a series of mergers of similar-sized bodies) as they migrated away from the rings⁴,⁶. This scenario is supported by the observations of the small inner Saturnian satellites, as bodies with similar semi-major axis have comparable masses: *m*<sub>Prometheus</sub>/*m*<sub>Pandora</sub> ~1.16 and *m*<sub>Janus</sub>/*m*<sub>Epimetheus</sub> ~3.6 (Jet Propulsion Laboratory Solar System Dynamics database, not accounting for error bars), and the mass of these bodies increase with their distance to the rings. The pyramidal regime provides an alternative to the formation by gradual accretion of small aggregates of ring material onto a protomoon. The latter scenario would result in Roche ellipsoids³,¹⁰, not consistent with the observed shapes of the small moons. For instance, while Atlas and Pan require a mechanism that makes their shapes flatter¹⁰, Prometheus is over-elongated and its long axis extends beyond the Roche lobe³.

Here, we develop a model for the late stages of the formation of the small inner moons of Saturn, from Pan to Janus and Epimetheus, assuming that they formed in the pyramidal regime⁴,⁶. We investigate if the shapes resulting from the final collisional mergers are consistent with the current shapes of these moons. To this end, we combine N-body simulations to estimate the possible range

of impact angles and velocities between the precursors of a given moon, and smooth particle hydrodynamics (SPH) simulations to obtain the outcome of the collisions.

For low-velocity impacts, which are expected given the range of eccentricities and semi-major axis of the small Saturnian satellites ($e < 0.01$, $a \sim 140,000$ km), the shape resulting from a collision depends on the impact angle $\theta$, the angle between the relative position and relative velocity of the moonlets, and the ratio $V_{\text{imp}}/V_{\text{esc}}$ between the impact velocity and the mutual escape velocity¹¹. Considering similar-sized moonlets of low density ($\rho = 700 \text{ kg m}^{-3}$) and of semi-major axis $a = a_{\text{Atlas}}$, this ratio can be estimated as $V_{\text{imp}}/V_{\text{esc}} \sim 0.58e/\mu^{1/3}$ (see Methods), where $e = \max(e_1, e_2)$ is the eccentricity of the moonlets before collision, and $\mu = (m_1 + m_2)/m_{\text{saturn}}$ is the mass ratio between the colliding moonlets and Saturn. Using the current masses and eccentricities of the small moons, we obtain $V_{\text{imp}}/V_{\text{esc}} \sim 3.0$ for $m_1 + m_2 = m_{\text{Atlas}} \sim 1.9$ for $m_{\text{Prometheus}}$ and $\sim 3.9$ for $m_{\text{Pandora}}$. There is a priori no reason to believe that the eccentricity during the mergers was comparable to the current eccentricity of the moons. We show, however, that it was probably the case: although the collisional growth of the moonlets efficiently damped their eccentricity, and the oblateness of Saturn prevented their excitation through

<figure><img src="image_2.png" /><figcaption>**Fig. 1** | Shapes of small moons of Saturn compared to model outcomes.</figcaption></figure>

Top: Cassini observations of Pan and Atlas revealed pronounced equatorial ridges. The elongated structure of Prometheus shows very different characteristics. Bottom: results of SPH calculations of low-velocity mergers of similar-sized moonlets. Left: $V_{\text{imp}} = 2V_{\text{esc}}$, $\theta = 2^\circ$, $m_1/m_2 = 1$; middle: $V_{\text{imp}} = 2.25V_{\text{esc}}$, $\theta = 0^\circ$, $m_1/m_2 = 1$; right: $V_{\text{imp}} = 2V_{\text{esc}}$, $\theta = 6^\circ$, $m_1/m_2 = 1$. White regions correspond to highly strained and/or ejected material (total integrated strain > 1), fine-grained material, whereas dark regions indicate less affected, original surface. Credit: top row, NASA/JPL-Caltech/Space Science Institute.

<footnote n="1">¹Physics Institute, Space Research and Planetary Sciences, Center for Space and Habitability - NCCR PlanetS - University of Bern, Bern, Switzerland.</footnote>

<footnote n="2">²IMCCE, Observatoire de Paris - PSL Research University, UPMC Univ. Paris 06, Univ. Lille 1, CNRS, Paris, France. *e-mail: adrien.leleu@space.unibe.ch</footnote>
<figure><img src="image_3.png" /></figure>

**Fig. 2** | Impact probability distribution and collision regimes. The solid line in panel (i) marks the transition between merging/non-merging collisions, defined as $V_{\text{imp}}/V_{\text{esc}} = 0.4 \sin(45^{\circ})/\sin(\theta^{[\circ]})$, corresponding to a constant initial angular momentum$^{1}$. Examples of shapes resulting from collisional mergers (Fig. 3; both mass ratios) are included. At very low impact speeds ($< -0.75V_{\text{esc}} $, grey area), distinct bilobate-type shapes are produced. These configurations are unstable and either separate due to tidal and Coriolis forces into two components, or might evolve into elongated bodies. The evolution of $10^{4}$ pairs of moonlets of mass $m_{1}/m_{2} = 6/5$ and $m_{1} + m_{2} = m_{\text{Atlas}}$ have been studied for three different distributions of initial eccentricities. Panels (ii) and (iii) show, respectively, the distributions in $V_{\text{imp}}/V_{\text{esc}}$ and $\theta$ for the moonlets that merged with $m_{1}/m_{2} < 7/4$, resulting in shapes comparable to those presented in panel (i). The fractions that reached $m_{1}/m_{2} > 7/4$ before the end of the simulation are not represented and we assume that the significant disruption/reaccretion they undergo results in random shapes. 30% of the systems with the $e_{\text{medium}}$ distribution (equivalent to $10^{5}$ years spent in the chaotic area, the mean value of the distribution being close to the current eccentricity of Atlas $< e_{\text{medium}} > -e_{\text{Atlas}}$) merged with $m_{1}/m_{2} < 7/4$, 50% for $e_{\text{low}}$ (equivalent to $-10^{4}$ years in the chaotic area, $< e_{\text{low}} > -e_{\text{Atlas}}/2$), and 18% for $e_{\text{high}}$ (equivalent to $-10^{6}$ years, $< e_{\text{high}} > -3e_{\text{Atlas}}/2$).

secular perturbation from the main moons (Methods), their eccentricity would increase before collision—as the semi-major axes of two moonlets get close, they enter a chaotic area due to the overlap of first-order mean-motion resonances$^{1,2,13}$. We studied the evolution of the eccentricity in this area in the case of $\mu_{\text{Atlas}}$, $\mu_{\text{Prometheus}}$ and $\mu_{\text{Janus}}$, using numerical integration of the three-body problem, taking into account the oblateness of Saturn (the effects of the rings and tides are negligible on the considered timescale$^{1,9,14}$) (Methods). We found that the excited eccentricities ($e_{\text{exc}}$) rapidly reach $3\mu^{1/3} < e_{\text{exc}} < 9\mu^{1/3}$, values that are comparable to the current eccentricities of the small moons. As the excited eccentricity is proportional to $\mu^{1/3}$, this leads to collisions in the range $0 < V_{\text{imp}}/V_{\text{esc}} \le 5$, regardless of the size of the considered moonlets. For Atlas-like objects, these calculations indicate typical impact velocities of order of a few $10 \text{ m s}^{-1}$.

We explore the outcome of such low-velocity (comparable to a few times the escape velocity) collisions among similar-sized moonlets using a three-dimensional SPH code$^{15-17}$ (Methods). Our code includes sophisticated models for material strength and porosity,

effects not considered in previous studies$^{5,10,18}$. We investigate a range of target-to-impactor mass ratios $m_1/m_2$, impact angles $\theta$ and velocities $V_{\text{imp}}$. The small moonlets are modelled as porous self-gravitating aggregates with densities of $\rho = 500-700 \text{ kg m}^{-3}$ and radii of $R_1 \sim 11 \text{ km}$. Since the small moons are located very close to Saturn, at or within the Roche distance, tidal effects are significant$^{18}$ and are included (Methods). Our three-dimensional SPH calculations show that only collisions with a close to head-on impact angle lead to merging into a stable structure (Figs. 1 and 2). Impacts with slightly larger impact angles lead to elongated, rotating structures that are split subsequently into two components. Because of the small Hill radius ($R_{\text{Hill}}/R_1 \sim 1.7$ at a distance of $\sim 1.4 \times 10^8 \text{ m}$, assuming a density of $\sim 700 \text{ kg m}^{-3}$), the two components quickly become unbound and separate. For even more oblique impact angles, collisions lead to hit-and-run events with not much interaction$^{19}$. The transition from merging to non-merging collisions can be defined as a function of the angular momentum$^{11}$ as given by the relative velocity $V_{\text{imp}}$ and the impact angle $\theta$ (Fig. 2).
<figure><img src="image_4.png" /><figcaption>**Fig. 3** | Shapes resulting from merging collisions. Shown are the results of SPH impact code calculations of low-velocity mergers among similar-sized moonlets, taking place within the Roche distance of Saturn. Collisions of bodies with an initial mass ratio of $m_1/m_2=1$ (left) and $m_1/m_2=1.5$ (right) are investigated for a range of impact angles $\theta$ and velocities $V_{imp}/V_{esc}$. The calculations start with initially spherical objects modelled as porous aggregates with initial bulk densities of $\rho=700 \text{ kg m}^{-3}$ and include self-gravitation as well as tidal forces. The resulting shapes are displayed as iso-density surfaces ($\rho_{iso}=500 \text{ kg m}^{-3}$). White regions correspond to highly strained/ejected material (total integrated strain $>1$), smoothed material; while dark regions indicate less affected, original surface. The simulations were carried out to $t \sim 17 \text{ h ($>1$ orbital period)}$. At this time, the objects produced in the collisional mergers are not yet aligned with the orbital plane, but will be reoriented later on due to the torques resulting from the tidal forces.</figcaption></figure>

For random coplanar orbits, the frequency function of the impact angle $\theta$ is given by $dP = \cos(\theta)d\theta$ (Methods). This distribution is hence maximum for head-on collisions, which favours the merging of the moonlets. However, as the merging domain is rather narrow (Fig. 2), it is probable that a moonlet pair undergoes a few hit-and-run collisions before the final merger, changing their mass ratio in the process. Directly before the first collision, and during the hit-and-run phase, the moonlets are on intersecting orbits in the inner part of the previously discussed chaotic area. Given that the timescale between collisions is significantly shorter than any perturbation (a few $10^4$ years, see Methods), we assume that the gravitational interaction between the moonlets dominates their dynamics. To estimate their mass ratio at the time of merging, along with $V_{imp}/V_{esc}$ and $\theta$, we hence designed a model that takes into account the gravitational interactions between the two moonlets and Saturn, the oblateness of Saturn, and the collisions between the moonlets. When a collision occurs, $V_{imp}/V_{esc}$ and $\theta$ are compared to the domains defined in Fig. 2(i). In the case of non-merging collisions, the mass ratio $m_1/m_2$ after the collision is derived from the SPH collision model (Methods), and the new velocities are computed considering an elastic collision with dissipation²⁰ (Methods). We point out that in the case of Pan, we would need to model the rings as well to check if the eccentricity of the moonlets is excited to the same degree before collision. We note, however, that the chaotic area and the radial excursion corresponding to the considered eccentricities is comparable to the size of the Encke division (~300 km). Pan's precursors could hence have been on slightly eccentric orbits before collision without encountering ring material.

The model was run for different initial distributions of eccentricity, with $m_1/m_2=6/5$ and $m_1+m_2=m_{\text{Atlas}}$ as initial condition.

However, it can be used to describe the formation of the other small inner moons as well, because for the size range of these small bodies, the results are scale-free for a given mean density and semi-major axis (Methods). In each case, we randomly selected $10^4$ initial conditions on intersecting orbits, such that 95% of the systems either merged or exceeded a threshold mass ratio within $10^5$ years. We set this threshold to $m_1/m_2 \sim 7/4$ as it is roughly the limit for which the formation of an equatorial ridge by merging of similar-sized bodies is still feasible (Fig. 3). Because the size of the moonlet is comparable to its Hill sphere, ejections due to scattering events are extremely rare. For the systems that exceeded the threshold mass ratio before the final merger, we assume that the disruption/re-accretion processes occurring during the multiple collisions tend to produce featureless/random shapes within the Roche lobe, as is the case for Pandora, Janus and Epimetheus.

The distribution of $V_{imp}/V_{esc}$ and $\theta$ for the mergers occurring with $m_1/m_2 < 7/4$ are shown in Fig. 2. These mergers occur typically after 1–2 merge-and-split or hit-and-run collisions. These collisions lead to a stable merged structure that result in a variety of shapes (Fig. 3). Regardless of the chosen initial eccentricity distribution (Fig. 2), close to head-on impacts lead to flattened structures with large-scale ridges, resembling the observed shapes of Atlas and Pan. Another frequent type of resulting structures is elongated shapes with characteristics similar to Prometheus, which result from impact parameters close to the boundary between merging and merge-and-split collisions. Depending on the chosen eccentricity distribution, 20 to 50% of the systems result in flat or elongated moons (these probabilities depend on the initial mass ratio between the precursors and increase as this mass ratio gets closer to one). These results are consistent with the fraction of small moons possessing peculiar shapes:
<figure><img src="image_5.png" /><figcaption>**Fig. 4** | Ridge formation on Iapetus. Top: Iapetus as observed by Cassini has an oblate spheroid shape and an equatorial ridge. Bottom: shown is the result of an SPH calculation of a head-on merger of two equal-sized bodies with half of the mass of Iapetus taking place with a velocity of $V_{imp}/V_{esc} = 1.25$, leading to a distinct equatorial ridge and an overall shape with an oblateness of the order of a few per cent, comparable with the observed one (-4%). Credit: top, NASA/JPL/Space Science Institute</figcaption></figure>

Pan, Atlas and Prometheus exhibit such features, while Pandora, Janus and Epimetheus do not. For completeness, the effect of inclination of the moonlets was also investigated. Inclination comparable to the current value for Atlas does not change significantly the ratio of successful merger, however, it shifts $V_{imp}/V_{esc}$ towards lower values as it slightly increases the probability of hit-and-run encounters (see Methods).

The extended ridges observed on Saturn's small moons represent very smooth terrains, while the main bodies not covered by the ridge show more rocky/rough terrain². In collisions, smooth terrain may result from ejecta reaccumulation as well as restructuring due to large deformations of the original surface (Methods). According to our model, for bodies resulting from ridge-forming mergers, such material (white in Fig. 3) is mainly focused at the equatorial areas. The bulk densities resulting from the mergers are around $\rho = 400-600$ kg m⁻³ (Methods), in agreement with the observations²¹. These results suggest that the spectrum of structural features observed on Saturn's small moons—from flattened objects with extensive ridges to over-elongated shapes—can be explained by the final pairwise accretion of comparable-sized moonlets, supporting the pyramidal regime formation scenario⁴,⁶, and suggest that no other significant process changed the shape of the small moons since the pyramidal

regime. Our study also implies that hit-and-run collisions between similar-sized objects were frequent throughout the formation of the small inner moons. Although such an event between Prometheus and Pandora was previously considered to explain the origin of the F ring⁷, in our model the F ring material could originate from the hit-and-run collision of the precursors of Prometheus (or Pandora), then be trapped in the stability area between Pandora and Prometheus after the merging of the precursors.

The unique shapes of the small Saturnian moons are a consequence of the near-coplanarity of the system. The current inclinations of Jupiter's small moons, for example, do not favour head-on collision and hence make the merger of similar-sized moonlets inefficient. The similarity between the shapes of Amalthea and Prometheus, however, suggests a common origin for these objects, and would imply that the inclination of Amalthea's orbit evolved after its formation, perhaps crossing resonances with Io²². The mid-sized Saturnian moons have also been proposed to have formed in the pyramidal regime⁴,⁶. Although the features we produce here are not observed on most of them, our study does not contradict this scenario: the current inclinations of these moons would favour multiple collisions before the final merger, which may produce featureless/random shapes. Moreover, the mass of these objects is large enough that their shapes have evolved closer to hydrostatic equilibrium. Finally, these moons are older and it is difficult to infer if the pyramidal regime would have been the last process to significantly affect their shapes.

Interestingly, Iapetus is the only large moon that displays an oblate shape and possesses an equatorial ridge, although this moon was not proposed to have formed in the pyramidal regime⁶. Its ridge was previously suggested to have formed either by endogenic or exogenic processes²³⁻²⁸. To explore if these features can result from the merging of similar-sized bodies, we performed additional collisional simulations with sizes and densities adjusted for the Iapetus case (Methods). We note that tidal forces are negligible at Iapetus' location. Our modelling results suggest that the oblate shape of Iapetus as well as the equatorial ridge may be a result of a merger of similar-sized moons, taking place with velocities around $V_{imp} = 1.2-1.5V_{esc}$ and a close to head-on impact angle (Fig. 4). The probability for such an impact is difficult to estimate, as it may not have occurred on the current orbit of the moon. We note, however, that if the peculiar shape of Iapetus is due to the merging of similar-sized bodies, possibly after a series of hit-and-run encounters, its precursors must have been on a similar orbital plane. In that case, Iapetus probably formed before gaining its significant inclination with respect to Saturn's equatorial plane, or out of bodies that had another common dynamical origin²⁹.

## Methods

**Orbital trajectories.** *Implementation of the forces.* The interaction of two moonlets orbiting Saturn was determined by the integration of the equation of the three-body problem in astrocentric coordinates, where the potential of the centre body was modified to take into account the first harmonic of the oblateness ($J_2$) of Saturn³⁰:

$$V_{\text{Sat}}(r, \varphi) = -\frac{Gm_{\text{Sat}}}{r} \left[ 1 - \frac{J_2}{2} \left( \frac{R_{\text{Sat}}}{r} \right)^2 (3\cos^2\varphi - 1) \right]$$

where $\varphi$ is the colatitude, $r$ is the distance from the moonlet to the centre of Saturn, $R_{\text{Sat}}$ is the equatorial radius of Saturn and $G$ is the universal constant of gravitation. The integrations were performed using the variable-step integrator DOPRI, which uses the Runge-Kutta method of order 8(7), setting an integration precision of $10^{-14}$.

**Secular perturbation from the main moons.** The secular perturbations induced by one of the main Saturnian moons on a moonlet near the rings of Saturn can be estimated by the Laplace-Lagrange solution of the secular problem generalized to take into account the oblateness of the central body³⁰. Considering the first order in $J_2$, we checked the effect of Mimas (the innermost main moon) and Titan (the largest), on both Atlas (closest to Saturn) and Epimetheus (the outermost considered small moon); in each case, the forced eccentricities due to the secular
perturbation are lower than 10⁻⁵. The effects of the main moons can hence be neglected as long as the considered semi-major axes are not too close to the 3:2 and 4:3 mean-motion resonances with Mimas (currently located at ~141,500 km and ~153,000 km, respectively).

Eccentricity excitation in the chaotic area. In the case of the coplanar three-body problem, for initially circular orbits, the moonlets will be in the chaotic area due to the overlap of first-order mean-motion resonances (MMR) if their semi-major axis $a_i$ satisfy¹²,¹³ ($\mu = (m_1+m_2)/m_{sat}$):

$$|a_1 - a_2|/a_1 \le 1.5\mu^2$$

and they are not protected by the 1:1 (co-orbital) MMR. At $a_i = a_{\text{Atlas}}$, this yields a half-width of ~150 km for $m_1 + m_2 = m_{\text{Atlas}}$, ~385 km for $m_1 + m_2 = m_{\text{Prometheus}}$ and ~860 km for $m_1 + m_2 = m_{\text{Fama}}$. The precession induced by the oblateness of Saturn displaces individual MMRs and might split them into corotation and Lindblad resonances that could enlarge the chaotic area³¹.

To estimate the orbital evolution of two moonlets in the chaotic area, we integrate the trajectories of a pair of moonlets for different initial values of $a_1 - a_2$. The two moonlets are initially placed on circular orbits such that the semi-major axis of the outer moonlet is the current semi-major axis of Atlas, and the initial mean longitudes were chosen randomly. The maximum time step is set to 1/20 of the orbital period, and the integrator precision to 10⁻¹⁴. The results are presented in Supplementary Fig. 1.

Supplementary Fig. 1a-c represent the evolution of the eccentricity (colour code in log scale) of the body $m_1$ as a function of the time and the initial value of $a_1-a_2$, in the case $m_1=m_2$. Three cases are represented: $m_1+m_2=m_{\text{Atlas}}$ (Supplementary Fig. 1a), $m_1+m_2=m_{\text{Prometheus}}$ (Supplementary Fig. 1b) and $m_1+m_2=m_{\text{Fama}}$ (Supplementary Fig. 1c). For these objects on initially circular orbit, the stable 1:1 mean-motion area extend roughly up to $a_1-a_2 \sim 0.7a_1\mu^{1/3}$, but decrease for larger eccentricities³². Supplementary Fig. 1d represents the mean value of the eccentricity of $m_1$ after a given time spent in the chaotic area for each pair of half-moons. On average, eccentricities comparable with the current eccentricities of these moonlets is reached within 10⁶ years.

Supplementary Fig. 1e represents the same quantity as in the Supplementary Fig. 1d, but with $1/t$ and $<e_1>$ now normalized by $\mu^{1/3}$. We found that the average eccentricities are excited proportional to $\mu^{1/3}$ after a time $\mu^{1/3}t$ spent in the chaotic area, a scaling that is present in the Hamiltonian of individual resonances at low eccentricities³³.

Estimation of the ratio $V_{\text{imp}}/V_{\text{esc}}$. The mutual escape velocity of two moonlets of equal mass $m_1=m_2=\mu m_{\text{sat}}/2$ is given by:

$$V_{\text{esc}} = \sqrt{\frac{2G(m_1 + m_2)}{R_1 + R_2}} = \sqrt{2G \left(\frac{4\pi\rho}{3}\right)^{1/3} m_{\text{Sat}}^{1/3} \left(\frac{\mu}{2}\right)^{1/3}}$$

where $R_1$ and $R_2$ are the radius of $m_1$ and $m_2$ in the case of spherical bodies and $\rho$ is the mean density of the moonlets. At the considered densities ($\rho \sim 500 \text{ kg m}^{-3}$), the Hill radius of the moonlets is comparable to their physical radius; the impact velocity $V_{\text{imp}}$ is hence comparable to the relative velocity. At a given eccentricity, the impact velocity can be estimated by:

$$V_{\text{imp}} \approx V_{\text{relative}} \approx nae = \sqrt{Gm_{\text{Sat}}/a^3 e}$$

where $n$ is the mean motion of the moonlet. We hence obtain:

$$\frac{V_{\text{imp}}}{V_{\text{esc}}} \approx m_{\text{Sat}}^{\frac{1}{6}} \left( \frac{4\pi\rho}{3} \right)^{\frac{1}{6}} \frac{1}{2^{\frac{1}{6}}\sqrt{a}} \frac{e}{\mu^{1/3}}$$

Taking the semi-major axis of the considered small moonlets outside the main ring $a \sim 137,500$ km and $\rho = 700 \text{ kg m}^{-3}$, we obtain $V_{\text{imp}}/V_{\text{esc}} \sim 0.58e/\mu^{1/3}$.

**Collision modelling.** To model the collisions we use a SPH impact code¹⁶,¹⁷,¹⁸. This code includes self-gravity and is especially suited to modelling collisions between rocky/granular bodies.

We assume that the moonlets are spherical, cohesionless porous low-density aggregates. In the model, a pressure-dependent shear strength (friction) is included by using a standard Drucker-Prager yield criterion¹⁷. Granular flow problems (of cohesionless material) are well reproduced using this method¹⁷,¹⁸. A coefficient of friction $\mu = 0.8$ (ref. ³⁶) is used. Porosity is modelled using a P-alpha model with crush-curve parameters corresponding to pumice³⁷. We use the Tillotson equation of state with a reduced bulk modulus of $A = 2.7 \times 10^9 \text{ Pa}$ (leading term in the equation of state) to take into account the smaller elastic wave speeds in porous materials compared with solid rock. We note that for the impact velocities considered here (up to a few $10 \text{ m s}^{-1}$), no significant heating occurs and the bulk modulus together with the crush-curve parameters solely define the material response in compression.

$5 \times 10^5-10^6 \text{ SPH particles}$ are used in the simulation of the merging collisions. To explore the parameter space of non-merging collisions, low-resolution runs with $\sim 5 \times 10^4$ particles are performed.

To include tidal effects, we use the linearized equation of motion in a rotating Cartesian coordinate system³⁸ orbiting Saturn. The colliding small moonlets are assumed to have a tidally locked rotation state, that is, they have zero spin in the rotating coordinate system.

We use an $x$ axis pointing away from Saturn, the $y$ axis tangential to the direction of motion and the $z$ axis perpendicular to the orbital plane. This yields the following forces per unit mass, which are added as external accelerations in each time-step to the SPH particles:

$$dx^2/dt^2 = 3\Omega^2 x + 2\Omega\Omega dy/dt$$

$$dy^2/dt^2 = -2\Omega\Omega dx/dt$$

$$dz^2/dt^2 = -\Omega^2 z$$

where $\Omega^2 = GM_{\text{sat}}/a^3$, $G$ is the gravitational constant and $a$ is the distance to Saturn. In our nominal case, $a$ is set to $1.33 \times 10^9 \text{ m}$ and the impact geometry is oriented so that the bodies are moving along the $x$ axis before the collision. The effects of the strength of the tidal force as well as the orientation are shown in Supplementary Fig. 2.

Previous studies of small-body collisions suggest that structures resulting from large-scale collisions may have an increased porosity¹⁹,⁴⁰. We therefore start with bodies with initial bulk densities ($\rho = 700 \text{ kg m}^{-3}$) that are slightly larger than those observed for the small moons ($\rho \sim 500 \text{ kg m}^{-3}$). It has been suggested that the small moons formed from accretion of ring debris around a denser core of $\rho \sim 900 \text{ kg m}^{-3}$, but the shapes of the moons are also consistent with a homogenous structure¹⁹. Here, we assume an intermediate density (and initially homogenous structure) for the colliding moonlets. The structures resulting from merging collisions have non-homogenous density distributions and average densities lower than the initial ones, as discussed below.

We first investigate the possible evolution of the spherical shapes of the individual bodies under the influence of self-gravitation and tidal forces. We find that these bodies, supported by the pressure-dependent shear strength (friction), keep their spherical shapes, and we therefore use initially spherical bodies in our simulations.

After non-merging hit-and-run collisions, the bodies may not be perfectly spherical any more, which could affect the outcome of the subsequent merger. Due to computational limitations, it is not possible to systematically track and compute all the shapes resulting from hit-and-run events. We compute for one characteristic case (with $m_1/m_2 = 1.125$, $V_{\text{imp}}/V_{\text{esc}} = 2$, $\theta = 45^\circ$) the post hit-and-run shape of the largest remnant. This shape is then used to study the outcome of a merger of two irregularly shaped bodies (using $m_1/m_2 = 1$, $V_{\text{imp}}/V_{\text{esc}} = 2$, $\theta = 90^\circ$). We find that the main features (large-scale ridges) are similar as in the case of a merger of spherical bodies with otherwise the same initial conditions (Supplementary Fig. 3).

**Post-impact analysis.** Properties of merged bodies. We compute the total strain experienced by the material during the final collision by integrating the second invariant of the strain-rate tensor for each SPH particle (over the entire simulation). We use this measure to distinguish between highly sheared (integrated strain > 1), possibly fined-grained material, and material that was less affected (integrated strain < 1), and kept its original properties. Material that was ejected and reaccumulated on the final merged structure experienced high strains ($\gg 1$). We note that the potential (small) deformations experienced during the previous hit-and-run events are not considered in these calculations.

As found in previous studies¹¹, the low-velocity collisions considered here do not lead to any significant compaction. However, several mechanisms, such as shear dilatation or the ejection of material followed by reaccumulation, can lead to a final bulk density that can be significantly lower than that of the colliding bodies³⁹,⁴⁰.

To estimate the final density of the merged structure, we assume that the highly strained material ($> 1$) experienced an increase of porosity either due to shear dilatation⁴¹ or the process of reaccumulation that introduces macroporosity³⁹,⁴⁰. The corresponding material is assigned an increase of distention by $1/(1-0.4) \sim 1.7$, following ref. ³⁹. The resulting average densities depend on the impact parameters (the density generally decreases with increasing impact velocity) and the resulting structure. For initial densities of $\rho = 700 \text{ kg m}^{-3}$, the structures formed with $V_{\text{imp}}/V_{\text{esc}} > 1.75$ (Fig. 3) are in the range of 400–600 kg m⁻³.

We note that this analysis is performed as a post-processing step; the bodies shown in Figs. 1–4 are based on the original densities. As a result, the material shown in white would have a larger volume and therefore the ridges would be slightly more pronounced.
**Change of mass ratio in collisions.** In the case of the non-merging collisions, the post-impact mass ratio (mr') depends on the initial mass ratio (mr) as well as the impact parameters. We perform a set of simulations using different initial mass ratios to determine the change of the mass ratio in non-merging collisions as a function of impact angle and velocity. This allows for an empirical determination of mr' = f(mr) in different areas of the parameter space with varying degree of interaction.

The post-impact mass ratio is estimated by computing the masses of the largest two remnants after separation, considering SPH particles that are part of the individual bodies. For this, all accreted SPH particles with $\rho / \rho_0 > 0.7$ (where $\rho_0$ is the initial bulk density) are considered in the calculation of the mass; the lower density ejecta are not counted as part of the body.

While the transition from merging to non-merging collisions is relatively well defined by the angular momentum given by the collision geometry and velocity, the mass ratios resulting from the non-merging collisions have a more complex dependence on the impact parameters. In addition to the angular momentum, the strong tidal and Coriolis forces lead to a highly nonlinear behaviour. Collisions close to the merging/non-merging boundary typically lead to an initial merging of the two bodies, resulting in an unstable configuration that leads to a subsequent splitting into two components. On the other hand, collisions taking place far from the transition line show the characteristics of a typical hit-and-run event, with not much interaction between the two bodies.

The relative change of the mass ratio in the non-merging collisions is shown in Supplementary Fig. 4 as a function of impact angle and normalized velocity. We can roughly define three zones to characterize the relative change of the mass ratio in a non-merging collision, corresponding to different degrees of change.

For each zone, we compute the average relative change of the mass ratio. Supplementary Fig. 5 shows the dependence of the relative change of mass ratio as a function of the initial mass ratio for the three different zones. We find that the linear function $(mr' - mr)/mr = ar(mr - 1)$ reproduces well the dependence on the initial mass ratio in the considered range. Here, $mr = m_1/m_2$ is the mass ratio before the collision and $mr' = m_1'/m_2'$ is the mass ratio after the collision. The best-fit values are $ar = 0.11$ (zone 1), $ar = 0.65$ (zone 2) and $ar = 2.81$ (zone 3). Using the linear function as defined above for the different zones, we can compute the post-collision mass ratio as follows:

$$mr' = mr(ar(mr-1) + 1)$$

**Merging probability.** Hit-and-runs. The collisions occurring in the hit-and-run regime are treated as inelastic collision with dissipation$^{30,31}$ with a damping factor of 0.8 in the normal direction and 1 in the tangential direction, typical for the collision of rubble piles of the considered density$^{30,42}$. The new mass ratio between the moonlets is computed using the relations given in the previous section.

**Initial conditions.** We explore the effect of a different distribution of eccentricity that accounts for different durations spent in the chaotic area. For the cases displayed in Fig. 2, $10^4$ initial conditions were randomly chosen as following: the initial eccentricities $e_1$ and $e_2$ follow the distributions $e_1(t) = 10^k$ yr) and $e_2(t) = 10^k$ yr) that represent the eccentricity distributions of two bodies on initially circular orbits that remained for $t = 10^k$ years in the chaotic area, where $k$ is a positive integer. It is the mean values of these distributions that are represented in Supplementary Fig. 1d. The semi-major axes are then randomly chosen in order for the moonlets to be on intersecting orbits: $a_1 = a_{\text{Atlas}}$ and $a_2$ in $[a_1 + R_{\text{tilt}} + R_{\text{impact}}]$ $a_1(1 + e_1 + e_2) + R_{\text{tilt}} + R_{\text{impact}}]$. Finally, the mean longitudes and longitude of pericentere are randomly chosen in $[0, 2\pi)$.

**Comparison of the Atlas and Prometheus cases.** We verify that the dynamical study is scale-free for a given mean density and semi-major axis by comparing the outcomes of the cases $m_1 + m_2 = m_{\text{Atlas}}$ and $m_1 + m_2 = m_{\text{Prometheus}}$. The initial distribution of eccentricity is taken at $t = 10^5$ years in the case of Atlas (the nominal case described in the paper) and the corresponding $t = (m_{\text{Atlas}} / m_{\text{Prometheus}})^{1/10^5}$ years in the case of Prometheus. $2 \times 10^3$ initial conditions are chosen as described in the previous section for the Prometheus case, to compare with the $10^4$ initial conditions for Atlas. $R_{\text{impact}}$ is set to 66 km for Prometheus. With an initial mass ratio $m_1/m_2 = 6/5$, 29% of the runs ended up with a mass ratio below 7/4, against 30% for Atlas, with roughly identical distributions in $V_{\text{imp}}/V_c$ and $\theta$ (Supplementary Fig. 6).

**Effect of inclination.** For a random impact in 3D, the probability function $dP$ for the impact angle $\theta$ is given by$^{43}$ $dP = \sin(2\theta)d\theta$, giving a maximum probability of impact at $\theta = 45^\circ$. Let us adapt the reasoning to the 2D case.

Consider a meteoroid approaching a massless sphere of radius $r$. If both objects are coplanar, within a segment of length $2r$, all points of intersection of the path of the meteoroid with a plane perpendicular to its path are equally probable, and the total probability $P$ is given by:

$$P = 2r = 1$$

The differential probability with which the meteoroid will pass through a point at a distance $x$ from the centre of this segment, where $0 \le x \le r$, is:

$$dP = 2dx$$

and we have:

$$x = r\sin(q)$$

As $dx = r \cos(\theta) d\theta$, we obtain $dP = \cos(\theta) d\theta$, yielding a maximum probability for head-on impact $\theta = 0^\circ$.

In terms of impact angle probability, the system can be considered fully 2D as long as the vertical excursion of the bodies is significantly smaller than their physical size. In the case of Atlas, the current inclination$^{41}$ $i_{\text{Atlas}} = 0.003^\circ$ leads to a vertical excursion of ~7 km from the midplane, which is comparable to the size of the body. This results in a slight shift of the maximum probability of $\theta$ from zero to a few degrees. Supplementary Fig. 7 shows the distributions of $V_{\text{imp}}/V_c$ and $\theta$ in the case where the moonlets had inclinations comparable with the current inclination of Atlas before collision. The maximum of the distribution in $\theta$ is shifted to a few degrees, while $V_{\text{imp}}/V_c$ is shifted towards smaller values.

**Iapetus modelling.** We investigate potential ridge-forming collisions at large, lapetus-sized, scales using a relatively narrow range of mass ratios from 1:1 to 5:4; impact angles ranging from 0 to $10^\circ$ and velocities $V_{\text{imp}}/V_{\text{esc}}$ from 1 to 1.5 (Supplementary Figs. 8 and 9). We use homogenous non-porous bodies with densities of 1,100 kg m$^{-3}$ and masses of $0.9 \times 10^{21}$ kg, roughly half that of Iapetus. Although the spatial resolution (~13 km) in our numerical calculation does not allow for a comparison of the small-scale features of the ridge produced with the observed ones$^{45}$, our analysis suggests that the current shape as well as the equatorial ridge may be a result of a merger of similar-sized moons, taking place with velocities around $V_{\text{imp}} = 1.2-1.5V_{\text{esc}}$ and a close to head-on impact angle.

**Data availability.** The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.

Received: 2 November 2017; Accepted: 11 April 2018; Published online: 21 May 2018

References

1. Thomas, P. C. Sizes, shapes, and derived properties of the saturnian satellites after the Cassini nominal mission. *Icarus* **208**, 395–401 (2010).

2. Charnoz, S., Brahic, A., Thomas, P. C. & Porco, C. C. The equatorial ridges of Pan and Atlas: terminal accretionary ornaments? *Science* **318**, 1622–1624 (2007).

3. Porco, C. C., Thomas, P. C., Weiss, J. W. & Richardson, D. C. Saturn's small inner satellites: clues to their origins. *Science* **318**, 1602–1607 (2007).

4. Charnoz, S., Salomon, J. & Crida, A. The recent formation of Saturn's moonlets from viscous spreading of the main rings. *Nature* **465**, 752–754 (2010).

5. Hyodo, R. & Ohtsuki, K. Saturn's F ring and shepherd satellites a natural outcome of satellite system formation. *Nat. Geosci.* **8**, 686–689 (2015).

6. Crida, A. & Charnoz, S. Formation of regular satellites from ancient massive rings in the Solar System. *Science* **338**, 1196–1199 (2012).

7. Porco, C. C. et al. Cassini Imaging Science: Initial results on Phoebe and Iapetus. *Science* **307**, 1237–1242 (2005).

8. Thomas, P. C. et al. Shapes of the saturnian icy satellites and their significance. *Icarus* **190**, 573–584 (2007).

9. Goldreich, P. & Tremaine, S. Disk-satellite interactions. *Astrophys. J.* **241**, 425–441 (1980).

10. Karjalainen, R. Aggregate impacts in Saturn's rings. *Icarus* **189**, 523–537 (2007).

11. Jutzi, M. & Asphaug, E. The shape and structure of cometary nuclei as a result of low velocity collisions. *Science* **348**, 1355–1358 (2015).

12. Wisdom, J. The resonance overlap criterion and the onset of stochastic behavior in the restricted three-body problem. *Astron. J.* **85**, 1122–1133 (1980).

13. Deck, K. M., Payne, M. & Holman, M. J. First-order resonance overlap and the stability of close two-planet systems. *Astrophys. J.* **774**, 129 (2013).

14. Poulet, F. & Sicardy, B. Dynamical evolution of the Prometheus-Pandora system. *Mon. Not. R. Astron.Soc.* **322**, 343–355 (2001).

15. Benz, W. & Asphaug, E. Simulations of brittle solids using smooth particle hydrodynamics. *Comput. Phys. Commun.* **87**, 253–265 (1995).

16. Jutzi, M., Benz, W. & Michel, P. Numerical simulations of impacts involving porous bodies. I. Implementing sub-resolution porosity in a 3D SPH hydrocode. *Icarus* **198**, 242–255 (2008).

17. Jutzi, M. SPH calculations of asteroid disruptions: the role of pressure dependent failure models. *Planet. Space Sci.* **107**, 3 (2015).
18. Hyodo, R. & Ohtsuki, K. Collisional disruption of gravitational aggregates in the tidal environment. *Astrophys. J.* **787**, 56 (2014).

19. Asphaug, E., Agnor, C. B. & Williams, Q. Hit-and-run planetary collisions. *Nature* **439**, 155–160 (2006).

20. Richardson, D. C. Tree code simulations of planetary rings. *Mon. Not. R. Astron. Soc.* **269**, 493–511 (1994).

21. Thomas, P. C. Sizes, shapes, and derived properties of the saturnian satellites after the Cassini nominal mission. *Icarus* **208**, 395–401 (2010).

22. Hamilton, D. P., Proctor, A. L. & Rauch, K. P. An explanation for the high inclinations of Thebe and Amalthea. *Bull. Am. Astron. Soc.* **33**, 1085 (2001).

23. Castillo-Rogez, J. C. et al. Iapetus' geophysics: rotation rate, shape, and equatorial ridge. *Icarus* **190**, 179–202 (2007).

24. Singer, K. N. & McKinnon, W. B. Tectonics on Iapetus: despinning, respinning, or something completely different? *Icarus* **216**, 198–211 (2011).

25. Ip, W.-H. On a ring origin of the equatorial ridge of Iapetus. *Geophys. Res. Lett.* **33**, L16203 (2006).

26. Levison, H. F., Walsh, K. J., Barr, A. C. & Dones, L. Ridge formation and de-spinning of Iapetus via an impact-generated satellite. *Icarus* **214**, 773–778 (2011).

27. Dombard, A. J., Cheng, A. F., McKinnon, W. B. & Kay, J. P. Delayed formation of the equatorial ridge on Iapetus from a subsatellite created in a giant impact. *J. Geophys. Res.* **117**, E03002 (2012).

28. Stickle, A. M. & Roberts, J. H. Modeling an exogenic origin for the equatorial ridge on Iapetus. *Icarus* **307**, 197–206 (2018).

29. Asphaug, E. & Reufer, A. Late origin of the Saturn system. *Icarus* **223**, 544–565 (2013).

30. Murray, C. D. & Dermott, S. F. *Solar System Dynamics* (Cambridge Univ. Press, Cambridge, UK, 1999).

31. El Moutamid, M., Sicardy, B. & Renner, S. Coupling between corotation and Lindblad resonances in the presence of secular precession rates. *Celest. Mech. Dyn. Astron.* **118**, 235–252 (2014).

32. Leleu, A., Robutel, P. & Correia, A. C. M. On the coplanar eccentric non-restricted co-orbital dynamics. *Celest. Mech. Dyn. Astron.* **130**, 3 (2018).

33. Petit, A. C., Laskar, J. & Boué, G. AMD-stability in presence of first order mean motion resonances. *Astron. Astro.* **607**, A35 (2017).

34. Benz, W. & Asphaug, E. Simulations of brittle solids using smooth particle hydrodynamics. *Comput. Phys. Commun.* **87**, 253–265 (1995).

35. Jutzi, M., Holsapple, K. A., Wuenneman, K. & Michel, P. in *Asteroids IV* (eds Michel, P., DeMeo, F. & Bottke, W. F.) 679–699 (Univ. Arizona, Tucson, USA, 2015).

36. Collins, G. S., Melosh, H. J. & Ivanov, B. A. Modeling damage and deformation in impact simulations. *Meteorit. Planet. Sci.* **39**, 217–231 (2004).

37. Jutzi, M., Michel, P., Hiraoka, K., Nakamura, A. M. & Benz, W. Numerical simulations of impacts involving porous bodies: II. Comparison with laboratory experiments. *Icarus* **201**, 802–813 (2009).

38. Wisdom, J. & Tremaine, S. Local simulations of planetary rings. *Astron. J.* **95**, 925–940 (1988).

39. Jutzi, M., Benz, W., Toliou, A., Morbidelli, A. & Brasser, R. How primordial is the structure of comet 67 P/C-G? *Astron. Astro.* **597**, A61 (2017).

40. Jutzi, M. & Benz, W. Formation of bi-lobed shapes by sub-catastrophic collisions. *Astron. Astro.* **597**, A62 (2017).

41. Collins, G. S. Numerical simulations of impact crater formation with dilatancy. *J. Geophys. Res. Planets* **119**, 2600–2619 (2014).

42. Leinhardt, Z. M. & Richardson, D. C. & Quinn, T. Direct N-body simulations of rubble pile collisions. *Icarus* **164**, 133–151 (2000).

43. Shoemaker, E. M. in *Physics and Astronomy of the Moon* Ch. 8 (Academic Press, New York, 1962).

44. Spitale, J. N., Jacobson, R. A., Porco, C. C. & Owen, W. M. Jr The orbits of Saturn's small satellites derived from combined historic and Cassini imaging observations. *Astron. J.* **132**, 692–710 (2006).

45. Dombard, A. J., Cheng, A. F., McKinnon, W. B. & Kay, J. P. Delayed formation of the equatorial ridge on Iapetus from a subsatellite created in a giant impact. *J. Geophys. Res.* **117**, E03002 (2012).

## Acknowledgements

The authors acknowledge support from the Swiss NCCR PlanetS and the Swiss National Science Foundation.

## Author contributions

A.L. performed the dynamical modelling and analysed the results. M.J. performed the collision modelling and analysed the results. M.R. contributed initial ideas

for the study. All authors contributed to interpretation of the results and the preparation of the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

*Supplementary information* is available for this paper at https://doi.org/10.1038/s41550-018-0471-7.

*Reprints and permissions information* is available at www.nature.com/reprints.

*Correspondence and requests for materials* should be addressed to A.L.

*Publisher's note:* Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.