© The Authors 2025

1 Introduction

The existence of binary asteroid systems with small primaries is largely explained by rapidly rotating rubble-pile asteroids undergoing structural failure and shedding mass via ejection of boulders from the equator (Walsh et al. 2008, 2012), contact binary fission (Jacobson & Scheeres 2011), or ejection of a large chunk of mass (Tardivel et al. 2018), creating a satellite. The rapid rotation is likely caused by gradual spin-up via the Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect, where anisotropic re-emission of absorbed sunlight generates a weak torque that can either spin the asteroid up or down over tens of thousands to millions of years, depending on its size, shape, and distance from the Sun (Rubincam 2000). However, while these formation models explain the predominantly prolate (elongated) population of asteroid satellites (Pravec et al. 2016), the recent observations of two atypically shaped moons in binary asteroid systems have proved difficult to explain. First, NASA’s DART mission directly imaged the spherically oblate moon Dimorphos orbiting the asteroid Didymos (Rivkin et al. 2021; Daly et al. 2023; Chabot et al. 2024). Soon thereafter, the Lucy spacecraft observed what appeared to be a bilobate contact binary satellite now called Selam during a fly-by of the Dinkinesh system (Levison et al. 2021, 2024). Their shapes have proven difficult to reproduce with a single-moonlet formation model, as tidal forces render objects formed at typical distances from their primary prolate (Porco et al. 2007; Tiscareno et al. 2013). Models involving circumasteroidal debris discs caused by chaotic mass-shedding events provide an alternative path to explain atypically shaped asteroids (Madeira & Charnoz 2024; Agrusa et al. 2024; Wimarsson et al. 2024), as such a dynamically rich environment can lead to low-velocity mergers, tidal disruptions, and subsequent mass-accretion that deform satellites from their initially prolate shapes. For a more in-depth introduction to binary asteroid formation, we refer to Agrusa et al. (2024) and Wimarsson et al. (2024). These authors explore the formation of atypically shaped asteroid satellites using discrete element methods (DEMs).

More specifically, performing simulations with a polyhedral N-body code, Wimarsson et al. (2024) find that discs with a width of more than two primary radii form multiple moonlets that dynamically interact to form oblate spheroids and bilobate objects. Analysing mergers occurring in these debris disc simulations, they found that the initial non-spherical shapes, orientation, and mass difference between the moonlets highly affects the final shape. Crucially, they further identified that moonlet-moonlet mergers occurring in this regime have impact velocities, υ_impact, below the mutual escape velocity, υ_esc, leading to little initial deformation of the separate moonlets. In contrast, other studies investigating low-speed mergers between similar-sized rubble-pile structures, forming, for example, moons around Saturn (Leleu et al. 2018), bilobate cometary nuclei (Jutzi & Asphaug 2015; Jutzi & Benz 2017), contact binary asteroids (Marohnic et al. 2021), or asteroid satellites (Raducan et al. 2025), mainly consider velocities of the order of one or a few escape velocities. The very low velocities observed by Wimarsson et al. (2024), with the highest being 0.86υ_esc, can be explained by moonlets being on similar orbits, with eccentricities kept low by gravitational interactions and dissipative contacts with surrounding debris in the disc, leading to small impact angles ≤21.5°. As the impact velocity for moonlets A and B is proportional to µ^−1/3, where µ = (m_A + m_B)/M_primary is the mutual mass ratio relative to the primary (Leleu et al. 2018), this leads to the observed soft collisions. Such low impact velocities, of the order of or below υ_esc, have also been observed between clumps formed in gravitationally unstable circumstellar discs during planet formation simulations (Schib et al. 2025a,b). To our best knowledge, strictly sub-escape-velocity mergers between rubble piles remain an unexplored regime that could explain the presence of atypically shaped satellites in binary asteroid systems.

To fully understand the dynamical history of observed asteroid satellites, it is also crucial to take reshaping via tidal forces into account, given how susceptible rubble piles are to elongation and mass shedding. When assessing the stability of an object subject to tidal forces, a common approximation is to treat the body as a fluid, according to the analysis by Roche (1847). The fluid Roche limit is the closest approach where such a body can stay structurally intact. However, due to non-negligible friction levels in rubble piles, they exhibit a higher survivability and are less prone to reshaping. There have been a significant number of studies investigating tidal disruption and the distortion of rubble piles using numerical (Asphaug & Benz 1994, 1996; Richardson et al. 1998; Bottke et al. 1999; Walsh & Richardson 2006) as well as analytical methods (Holsapple & Michel 2006, 2008), and they have laid the groundwork for the current understanding. A general consensus of this research is that outcomes from tidal encounters depend on the distance to the primary object, the material’s angle of friction, particle size frequency distribution (SFD), cohesion, and the object shape. More recent works have further identified that particle packing, shape, and simulation resolution also have a non-negligible effect on the outcome of tidal disruption events of rubble piles (Movshovitz et al. 2012; Zhang & Lin 2020; Zhang & Michel 2020; Marohnic et al. 2023).

With these results in mind, it is clear that tidal evolution of rubble piles is a dynamically rich subject, given how sensitive the outcome is to the specific properties of a model. Moreover, a new dimension of complexity is added to the problem when involving non-homogeneous shapes such as bilobate structures. In this work, we make a first attempt to narrow down the formation of atypically shaped asteroid moons via sub-escape-velocity mergers and their subsequent structural evolution under the influence of tides with the use of numerical N-body simulations. Our formation scenarios are based on the results from the debris disc model of Wimarsson et al. (2024), henceforth referred to as W24. In Sect. 2, we introduce the numerical model and simulation setup, while Sect. 3 contains the results from our investigation. Section 4 provides analysis and discussion of how our results relate to the different aspects of our model, as well as their implications for potential formation scenarios of Dimorphos and Selam. Finally, Sect. 5 gives our conclusions and outlook for future studies.

2 Method

To track deformation of the moons resulting from the mergers while capturing key granular dynamics and long-lasting contacts between boulders, we chose to base our simulations on the gravitational N-body code GRAINS (Ferrari et al. 2017, 2020), which uses irregularly shaped particles. We go through the numerical details in the following section, and this is followed by the implementation of the mergers in simulations.

2.1 Numerical setup

The GRAINS code is based on the libraries of Chrono::Engine, an open source C++ physics engine (Chrono, Tasora et al. 2016) and can use non-spherical particles to study rubble-pile dynamics in space. The angular nature of the particles allowed us to replicate the complex physical effects of granular media such as particle–particle interlocking, particle spin orientation, off-centre collisions, and polyhedral contacts. Several studies have explored the benefit of using non-spherical particles, concluding that it enhances the ability to capture key dynamical behaviours of rubble-pile aggregates (Korycansky & Asphaug 2006, 2009; Ferrari & Tanga 2020; Marohnic et al. 2023). GRAINS is a DEM code that can handle tens of thousands of particles due to GPU acceleration and has been used to study rubble-pile aggregation (Ferrari et al. 2017; Ferrari & Tanga 2020) and rotational equilibrium (Ferrari & Tanga 2022), the dynamical state of the Didymos binary system (Agrusa et al. 2022; Ferrari et al. 2022), and the formation of binary asteroid systems post rotational failure (W24). Furthermore, Burnett et al. (2024) verify that it can be used to study tidal dissipation in binary asteroid systems. To resolve physical contacts between particles, Chrono offers two different options: non-smooth that is impulse-based and smooth that is force-based. The former performs well for short-lasting, high strain-rate events, while the smooth-contact method (SMC) is better suited for longer-duration contacts and low-velocity encounters, making it optimal for the purposes of this study. Moreover, using GRAINS with SMC makes it behave comparably to other DEM N-body codes that resolve contacts with soft-sphere methods (Sánchez & Scheeres 2011; Schwartz et al. 2012).

The numerical configuration for our simulations is practically identical to the one in W24. We set a static time step of 0.2 s, such that it was smaller than the characteristic time of the fundamental frequency of the spring-dashpot contact model (Ferrari et al. 2020) and used the default parallel integrator of Chrono: a Nesterov accelerated projected gradient descent method (Tasora 2017). To resolve contacts with SMC, we employed the non-linear Hertzian version of the two-way normal-tangent spring–dashpot system. The resulting forces due to contacts between two bodies depend on a set of constants that define physical properties of the material used. These include static and sliding friction, stiffness, damping, adhesion, and restitution, which have been kept at their default values determined from comparisons with laboratory studies of granular media (Tasora et al. 2016). Aggregates in our simulations are strictly kept together due to gravity and granular dynamical mechanisms such as interlocking, meaning the adhesion is kept at zero. Note that contacts with the primary body in the system have been omitted to speed up our simulations.

2.2 Particle generation

In GRAINS, each simulation is populated with randomly generated non-spherical particles created by filling a box with a pre-selected number of points, N_vertices, that act as vertices for a convex hull. The size of each box is determined by drawing a random diameter from the following exponential SFD: $$f_{D_{\text{p}}}\left(D_{\text{p}},D_{\text{min}},D_{\text{mean}}\right)=e^{-\left(D_{\text{p}}-D_{\text{min}}\right){\lambda }_D},$$(1)

where λ_D = 1/(D_mean − D_min). Due to the random nature of the vertex generation, some convex hulls will end up with fewer than N_vertices. For the simulations from W24, they employed an SFD with values of D_mean = 30 m and D_min = 5 m to keep the particle number of the generated debris disc of the order of 10³−10⁴. The mass of each particle is calculated based on its generated volume and a previously defined particle density, ρ_p, which was set to 1190 kg m⁻³ in W24. Throughout the remainder of this paper, we also make use of the property, ρ_bulk, which is the bulk density of a given aggregate.

Fig. 1 Visualisation of how new initial conditions are generated by rotating the most massive moonlet counterclockwise around its barycentre by an angle, α, here set to 45°. On the left is the original configuration from the W24 merger, with a vector â0 parallel to the aggregate’s semi-major axis. We performed the rotation for each particle, ending up with a new aggregate with its semi-major axis parallel to the vector â1.

2.3 Simulation setup

2.3.1 General configuration

The core part of our study revolved around altering the initial conditions of a specific merger and its subsequent dynamical evolution to better understand how pre-impact asymmetric, non-spherical shapes of the participating moonlets affect the shape of the final moon. The merger in question was observed during the later stages of a binary asteroid formation simulation based on the debris disc model in W24, namely Ryugu_run4, and resulted in an elongated, bilobate satellite (see their Sects. 3.3.1 and 5.2). For the remainder of this paper, we refer to this event as the W24 merger and the dynamical details of the event are provided in Table 1. Isolating the two moonlets involved in the merger from the surrounding debris disc 0.5 h prior to their impact, using the friends-of-friends clustering method introduced in W24, we obtained a more general scenario where we could disregard the dynamical influence of the disc and any smaller aggregates present in the system. The configuration is shown in the left plot of Fig. 1 while the moonlet data can be found in Table 2. The most massive moonlet, situated closer to the primary, is henceforth known as moonlet A, while the less massive, more spherically oblate body is moonlet B. The principal axes of the moonlets, along with any reference to principal axes in the remainder of the paper, have been computed with the dynamically equivalent equal volume ellipsoid (DEEVE) method used in, for example, Agrusa et al. (2024) and W24. The axes are denoted a, b and c and follow a > b, with c pointing out of the equatorial plane of the primary (for a full definition see Sect. 3.2 of W24). An oblate body has a principal axis ratio a/b close to unity, with Dimorphos having the value a/b = 1.06 ± 0.03 (Daly et al. 2024). To avoid confusion, we refer to the orbital semi-major axis as a_orb throughout the rest of this paper.

To generate a set of initial configurations for the two moon-lets that could help emphasise the effect of impact geometry for the post-merger shape of the moon, we imposed different rotations, α, onto the bodies in the equatorial plane of the primary. Information on how rotations are handled in GRAINS can be found in Appendix A. For a visual depiction of the rotation, see Fig. 1. Noting that rotating the second aggregate had little to no effect on the final shape of the merged object, we focused only on altering the pre-impact state of moonlet A. Each component of the aggregate was rotated around its own z-axis, while also transforming its position, velocity and angular velocity relative to the barycentre of the moonlet. The result was a shift in orientation of the rubble pile while keeping the magnitude and direction of its angular and linear velocity vectors intact. All-in-all, we generated 24 different setups with counter-clockwise rotations α ∈ [0, 360)° using an increment of 15°, with α = 0° being the fiducial configuration taken from W24 seen on the left side of Fig. 1.

2.3.2 Higher-resolution simulations

As mentioned, the dynamics and structural stability of rubble piles are highly dependent on particle shapes and polyhedral contacts (Korycansky & Asphaug 2006, 2009; Movshovitz et al. 2012; Ferrari & Tanga 2020; Marohnic et al. 2023), as well as SFD and boulder packing (Zhang et al. 2017, 2021; Zhang & Michel 2020; Raducan et al. 2024). Hence, it is likely that the impact and tidal forces affecting the final moon will lead to different degrees of distortion when employing finer or coarser SFDs in our system. To elaborate, when considering non-spherical particles, if there is a uniform distribution of large and small boulders, the aggregate will be rigid as compared to a case where the distribution is dominated by smaller particles. This comes from the fact that the small particles can fill out voids surrounding larger boulders and smooth out contacts between them, making the aggregate less rigid. This mechanism facilitates propagation of seismic waves in granular material (Marti et al. 2024). In turn, we set up simulations, using similar initial conditions but with aggregates of a steeper SFD with D_mean = 20 m and D_mean = 3 m, significantly increasing the particle number for each moonlet. As the new SFD is more similar to the observed SFDs on the surface of Dimorphos and Didymos (Barnouin et al. 2024; Robin et al. 2024; Pajola et al. 2024) than the old distribution, we chose not to treat our particles as porous collections of smaller sized boulders as in W24 and opted for a higher material density of 3000 kg m⁻³ (instead of 1190 kg m⁻³) more similar to that of LL chondrites (Flynn et al. 2018).

This was done by filling a box with a large cloud of 80 000 boulders with random positions and diameters, using the new SFD and material density. Letting the cloud gravitationally collapse into a clump, we could impose a three-dimensional graphical object model of each moonlet to create a higher-resolution moonlet with the correct shape. Each object model was created using the α-wrap library of CGAL (the computational geometry algorithms library, Alliez et al. 2022), essentially wrapping a thin ‘sheet’ around the moonlet, recursively subdividing its edges around the boulders until the algorithm has reached a target minimum length. Using built-in ray tracing capabilities of Chrono, we iterated through each body in the clump, sending out straight rays in the negative and positive directions of each axis from the centre-of-mass of the body. If all rays hit the alpha wrap mesh, the boulder was considered part of the new moon-let. As populating the shape with boulders from a finer SFD will lead to a different packing, we had to ensure the final configuration would have similar properties to the target moonlet in terms of bulk density and mass. This was achieved by exploiting the non-dimensionality of GRAINS. The resulting properties of each new aggregate, referred to as moonlets A2 and B2, can be found in Table 3. We point out that both aggregates had to be given the same bulk and material densities due to limitations of the implementation in GRAINS. For each simulation, we integrated the initial 0.05 h with a small time step of 0.01 s to let the body properly settle, avoiding any divergent solutions for contact computation. To account for the change in particle SFD, the rest of the simulation was performed with a time step of 0.05 s.

3 Results

3.1 Sub-escape-velocity mergers

Letting each initial configuration dynamically evolve for several orbits, it became clear that the removal of the disc significantly altered the resulting post-merger orbits for the newly formed moons, destabilising them. In Fig. 2, we have plotted the post-merger evolution of the main orbital elements for the original W24 case with a disc, as well as the scenario where the disc has been removed. The small initial offset in semi-major axis, eccentricity, and inclination can be explained by the simulations being started 0.5 h before the merger, while only the evolution of the largest remnant after the merger is shown. As the two moons evolved, the dynamical effect of the remaining disc, which at this point had a mass of 0.031 M_primary (including 0.023 M_primary in other moonlets), became substantial. The difference is clearly seen when comparing the change in eccentricity over time. The two small declines seen for the disc case at the 3.5 h and near the 5 h mark were caused by two minor tidal disruption events, where the body only lost a few per cent of its mass in total. Any further excitation of the orbit was then dampened by the disc. In contrast, the no-disc moon suffered a catastrophic tidal disruption between 4.5 and 8 h, losing more than half of its mass as the bilobate structure got pulled apart.

This accounts for the corresponding sharp increase in eccentricity and semi-major axis as only moonlet B remained, although on an eccentric orbit (we recommend watching the two complementary movies of this process¹). The final properties of the moon were m_moon = 0.35M₀, where M₀ is the post-merger mass of the progenitor, a/b = 1.48, b/c = 0.78 and a = 0.23R_primary. The merged moons for the remaining α values also suffer significant mass loss in a majority of the cases, leaving oblate or prolate remnants behind (see Fig. C.1 and Table D.1). All-in-all, this result indicates that bilobate moons formed in a disc have a higher survivability and that they can form closer to the primary than in the no-disc scenario.

To address the discrepancy in orbital behaviour, while still being able to investigate the stability and shapes of moons formed via sub-escape-velocity mergers, we decided to divide our study into two stages. First, we carried out each merger according to the 24 different initial conditions obtained via rotating moonlet A, starting each simulation 0.5 h before impact. Letting each formed moon settle after the merger, we paused the simulation after a total time of 3 h. Second, we pushed each moon further out from the primary, giving it a velocity corresponding to a circular orbit, leading to a new position of r = kr₀, where k is a factor larger than 1 and r₀ is the position relative the primary barycentre after 3 h. We cannot conclude that the mergers would lead to the same configurations when taking place further out from the primary without running new, self-consistent simulations. However, due to the chaotic nature of the debris disc formation scenario, along with the observed range of impact velocities recorded in W24, it is plausible that the shapes formed in our mergers could also form at these new distances kr₀. Even so, the results from the second stage of the study have to be interpreted with the removal of the disc in mind. To clarify, it only teaches us about the structural and orbital stability of a given object in isolation, as the presence of the disc would alter the dynamics of the problem and lead to further deformation via subsequent accretion, potentially homogenising the surface and shape as seen in Agrusa et al. (2024) and W24. That being said, in W24, the final shapes of the moons were largely determined by mergers, while subsequent accretion of single particles and small clusters had little effect on the final a/b ratio in comparison.

The initial shapes post-merger at the end of the first stage are shown in Fig. 3. For each case, we have coloured the particles depending on which moonlet they belonged to before impact. All figures depicting moon shapes in this work have been rendered in the inertial equatorial plane of the primary. For impact velocities higher than υ_esc, we would expect the initial shapes of the bodies to have little effect on the final configuration of particles as mergers in this regime lead to significant deformation (Leleu et al. 2018; Raducan et al. 2025). However, for the case of subescape-velocity, the deformation due to the impact is small, and we observe drastic differences in post-merger shape with just 15° rotations of the most elongated body, pre-merger, emphasising the influence of merger geometry. Among the resulting shapes, we can identify different types of bilobate structures, and we classify these using three different categories. The categories are quantified by γ ∈ (−90, 90]° that represents the angle between the longest principal axis of the most massive lobe DEEVE fit (see Appendix B) and the vector between the two barycentres of the lobes, r_AB, as shown in Fig. 4. The three categories are defined as follows:

• Long-axis bilobate (LAB), where |γ| < 10°, α = 0, 165, 180, 345°;

• Off-axis bilobate (OAB), where 10° ≤ |γ| < 45°, α = 15, 30, 45, 150, 195, 210, 225, 330°;

• Short-axis bilobate (SAB), where |γ| ≥ 45°, α = 60, 75, 90, 105, 120, 135, 240, 255, 270, 285, 300, 315°.

By this classification, Selam is potentially a LAB type (Levison et al. 2024), while Itokawa (Fujiwara et al. 2006) is an OAB binary. The seemingly bilobate nature of Apophis (Brozovic et al. 2018, 2022) would make it an SAB. We acknowledge that this methodology can only serve as a first attempt to classify these objects, as it fails to quantify how distinctly bilobate a merged object is. In turn, we aim to improve upon this classification system in further work.

Fig. 2 Evolution of semi-major axis (top), eccentricity (middle), and inclination (bottom) for the moons formed in the W24 merger with (solid blue) and without a disc present (dashed orange).

Fig. 3 Shapes and orientations of each formed moon after 3 h of simulation time after a sub-escape-velocity merger between the two moonlets in Table 2. The value α, represents the counter-clockwise rotation relative to the initial position of moonlet A in Fig. 1 around its barycentre. The particles originally belonging to moonlet A have been coloured in blue, while the corresponding particles from moonlet B are orange. The black arrow in each plot is pointing in the direction of the primary.

Fig. 4 Angle, γ, between the longest principal axis of the DEEVE of moonlet A (blue); a; and the vector rAB, which is the positional vector between the barycentres of moonlet A and moonlet B (orange).

Fig. 5 Same as Fig. 3 but after putting each body on a circular orbit with a semi-major axis of 1.1 r0 and evolving the system for an additional 45 h (48 h in total), where r0 is the distance from the primary barycentre at 3 h.

3.2 Survivability of formed objects

As at this point we could say little about the stability of each object, we proceeded to evolve the system in the second stage. Running test simulations without the disc, placing the α = 0° moon at different distances from the primary (further explored in Section 4.1), we found that the orbital and structural evolution of the object was highly similar to the disc-case when r ~ 1.1 r₀. Hence, for the first set of simulations, we let r₁ = 1.1r₀ = 2.673 R_primary for each α case and propagated the system for an additional 45 h. The resulting shapes after 48 h of total simulation are shown in Fig. 5. The physical and orbital properties of each final moon can be found in Table D.2, along with reference principal axis values for the shapes at the onset of this second stage of simulation from Fig. 3, subscripted by 0. We have attempted to classify the final shapes using our bilobate categories and γ prescription, along with the two additional types: oblate spheroid (OS), i.e. Dimorphos-like, and prolate spheroid (PS). The latter signifies a more compact shape where the bilobate structure of an elongated object can no longer be detected. Comparing the figures and properties of each moon, we observe how, throughout three orbits, the tidal influence from the primary has led to mass shedding and disruption in some cases, as well as homogenisation of some bilobate shapes, smoothing out their surfaces. For the cases of α = 15, 195, 210°, the disruption was catastrophic, resulting in mass loss. Studying their shapes in Fig. 5, we observe that most of moonlet A has been stripped off, along with some particles belonging to moonlet B, leaving behind an OS Dimorphos-like structure. While its DEEVE principal axes indicate a/b = 1.33, we recall that this method for evaluating physical shapes can sometimes overestimate the a/b value (Agrusa et al. 2024) and only serves as an estimate for the true physical extents. Formation of OS moons due to catastrophic tidal disruption was theorised in W24, but not shown until now. Even though this particular scenario leads to the object being put on an eccentric orbit, given the results from Fig. 2, it is plausible that the presence of a disc would be enough to dampen the orbit and circularise it. Furthermore, it also becomes evident that LABs and OABs are more susceptible to disruption and mass loss, given that all SABs stay intact. Instead, the tidal forces from the primary act to ‘straighten’ and smooth out the shapes of the short-axis merger objects (e.g. α = 120°), which is transformed from an SAB to a more homogeneous PS without any real indications of being produced by a merger. We can get a sense of how the shape changes over time by comparing the initial and final values for a/b and macroscopic porosity, estimated by evaluating the total volume of each aggregate with an α-shape scheme² based on the CGAL::Alpha_shape_3 library (Da et al. 2023). This module possesses the ability to compute an optimal α_3D value, which is the smallest value for which the algorithm will produce a ‘watertight’ mesh (for a more in-depth explanation, see Sect. 2.3.1 of W24). In turn, the porosity calculated using this method is the minimum macroscopic porosity possible for our aggregates. While the a/b values stay similar, except for extreme cases of either catastrophic disruption (α = 15, 195, 210°) or major reshaping (α = 225°), the porosity is generally significantly lower for each case, emphasising how the tidal forces serve to, with a few exceptions, homogenise the aggregates and make them less elongated at these distances.

To further improve our understanding of how the distance to the primary affected the final shapes of the merged moons, we also performed the same set of rotationally altered simulations for the case of r₂ = 1.2r₀. The corresponding shapes after 48 h of total simulation are shown in Fig. C.2 and their physical, geometrical, and orbital properties can be found in Table D.3. Notably, the moons undergo no mass shedding at these distances and mainly reshape due to homogenisation and spin, as a majority of them end up in super-synchronous rotational states (i.e. their rotational periods are shorter than their orbital periods). Studying the change in a/b values between the 3 h and 48 h mark, we find that the a/b-values consistently decrease for each α angle while the changes in porosity remain similar to that of the 1.1r₀ cases. When further shifting the post-merger distances to 1.3r₀, tides become even less influential, leading to additional reductions in a/b.

3.3 Higher-resolution shapes

Using the GPU-accelerated Barnes-Hut scheme of GRAINS, we carried out eight simulations for the moonlets in Table 3 with α increments of 45° up to the 3 h mark, as these simulations are computationally heavy with N_p = 35 729, excluding the primary (for visualisations of a few different mergers, see supplementary movies³). The initial bulk density of the cut-out aggregates did not match the target density of 626.8 kg m⁻³. Therefore, we had to scale the aggregate properties, which yielded a corresponding material density of 1134 kg m⁻³. The aggregates formed from each merger are shown in Fig. 6 and Table 4 shows the principal axes data, as well as porosities, for both the low-and high-resolution scenarios. Both from a visual comparison with Fig. 3 and the principal axes values, we observe a significant change in aggregate structure. The steeper SFD indeed serves to smooth out contacts between larger boulders, generally leading to less elongated shapes and more homogeneous surfaces. With an increase in resolution by a factor of almost 20, this also greatly contributes to the homogeneity of the surfaces. Furthermore, the DEEVE value estimates also become more accurate with larger values of N_p. Together with the smoother, less elongated shapes, this helps explain the discrepancy in a/b ratios for the objects, while the principal axes a remain largely similar. Yet, despite the reduced rigidity of the aggregate, the sub-escape-velocity nature of the mergers still leads to clear bilobate shapes. We can distinguish necks for the LAB cases α = 135 and 180° and the OAB shape α = 315°, as well as SAB structures for α = 45, 90, 225 and 270°. Moreover, after dynamically propagating the α = 0^° system in time for an additional 22 h at r₀ (25 h in total), we saw a similar catastrophic tidal disruption to the low-resolution case (shown in movie alpha0deg_post_merger_highres_627kgm3.avi), leading to a largest remnant on an eccentric orbit with m_moon = 0.40M₀, a/b = 1.30, b/c = 1.15 and a = 0.27 R_primary, as seen in the left image of Fig. 7. We also ran a comparison test for the same orbital properties and shape but higher bulk and particle densities of ρ_bulk = 1190 kg m⁻³ and ρ_p = 2154 kg m⁻³, giving the moon the same bulk density as the primary, which is predicted to be the case for the Didymos system (Daly et al. 2023). In this scenario, the moon did not undergo any tidal disruption event and retained all its mass (see movie titled alpha0deg_post_merger_highres_1190kgm3.avi). Deformation primarily occurred in the form of additional displacement of boulders. The final, more homogeneous shape can be seen in the right plot of Fig. 7 and has the properties: a/b = 1.60, b/c = 1.05, a = 0.41 R_primary and ϕ = 0.29. Tidal homogenisation and straightening made the body more elongated, similar to a prolate spheroid with an even less distinguishable boundary between the two moonlets, only evident from a shallow neck that can be detected when studying its α-wrap model.

Fig. 6 Higher-resolution versions of a few of the cases given in Fig. 3, with the moonlets consisting of 35 729 particles in total. The moonlet properties are given in Table 4.

Fig. 7 Final shapes for the case of α = 0° in Fig. 6 after an additional 21 h of simulation (25 h in total) for a bulk density of 626.8 kg m−3 (left) and 1190 kg m−3 (right).

4 Discussion

The variety in shapes obtained from small adjustments of our initial conditions reinforces our understanding that forming binary asteroid systems through the debris disc model is a highly chaotic and unpredictable process (Agrusa et al. 2024; W24). Notably, most of the pre-impact properties of the W24 merger, such as the mass ratio, impact velocity, distance, and angle, along with the moonlet shapes, were kept intact, with only the orientation of the most massive moonlet and particle SFD being altered. In turn, we have only explored a small portion of the extensive parameter space for mergers in the sub-escape-velocity regime. Looking at the distribution of properties for the major mergers observed in the W24 simulations, the most massive moonlet is always more or less prolate in nature and dominant in mass compared to the second (see their Sect. 5.2), with the chosen merger case for this study being close to the median values in both elongation and mass ratio. As for impact velocity, it is one of the most extreme cases with a value 0.84υ_esc, since it occurs closer to the primary than a majority of the other mergers. Despite this fact, we still see clear bilobate structures, even for the high-resolution simulations. Hence, a more typical impact velocity, near the median value of around 0.7υ_esc (based on the sample size of major mergers in W24) would result in even less deformation and more distinguishable bilobate features. In Fig. 8, we have plotted the scaling of the impact velocity with distance from the primary for the W24 merger case in Table 1. For the case of simplicity, we have assumed the largest moonlet to be on a circular orbit and scaled the velocity of the second moonlet according to their relative magnitudes and impact angle, i.e. rotating the velocity vector of the second, impacting moonlet in the direction of the primary relative to the circular orbit velocity. Given the low bulk densities of our moonlets, here set to a static 1190 kg m⁻³, we assume that the impact velocity is comparable to the relative velocity (Leleu et al. 2018). The three different graphs each represent a rubble-pile primary, all assumed to be spherical in nature, with the fiducial case of Ryugu in solid having ρ_mean = 1190 kg m⁻³ and R_primary = 500 m (Barnouin et al. 2019). The dashed line is representing Didymos with ρ_mean = 2400 kg m⁻³ and R_primary = 382.5 m (Daly et al. 2023), while the dotted line shows the velocities for Dinkinesh with ρ_mean = 2700 kg m⁻³ and R_primary = 395 m (Levison et al. 2024). It becomes clear that even for this extreme case with an impact angle of 21.5° occurring so close to the primary, we still get the sub-escape-velocity mergers for a much denser main body.

In turn, while we lack sufficient information regarding the full catalogue of potential shapes formed through sub-escape-velocity mergers, the samples from Figs. 3, 5 and 6 cover a significant fraction of them and identifies that observed contact binary shape types such as LAB (Selam-like), OAB (Itokawa-like) and SAB (Apophis-like) could be formed in orbit around rubble-pile primaries similar in size and density to Ryugu, Bennu, Didymos and Dinkinesh.

Fig. 8 Impact velocity in terms of mutual escape velocities for the case of the merger data in Table 1 and moonlets in Table 2, setting a static bulk density of 1190 kg m−3. The red square represents the values for the W24 merger.

4.1 Post-merger evolution and tidal distortion

After formation, the moons show susceptibility to deformation due to tidal forces from the primary via disruption and mass shedding or displacement of boulders, leading to homogenisation, elongation or all of them. Having explored three cases for the configuration in Tables 1 and 2 with varying distance from the primary barycentre, r₀ = 2.43R_primary, 1.1r₀ = 2.673R_primary and 1.2r₀ = 2.916R_primary, we observed variations in the level and nature of deformation, also for different shape types. Where several of the more initially elongated bilobate satellites fail structurally for the 1.1 r₀ case in Fig. 5 and undergo catastrophic disruption, other shapes retain all their mass. Hence, using analytical predictors for the tidal behaviours of these objects does not suffice, as initial shape and state of rotation alter the outcome even at similar distances from the primary, which is in agreement with the findings of Burnett et al. (2024), who investigate tidal dissipation in rubble-pile asteroid satellites. To better our understanding of the tidal disruption regimes for a given shape and particle composition, we performed a set of 45 h simulations for the case of the nominal α = 0°, where we altered its distance from the primary in increments of 0.01 r₀ between 1.01 r₀ and 1.25r₀. The resulting aggregates at the final time step for cases initiated up to 1.19r₀ from the primary have been plotted in Fig. 9. The geometrical and orbital properties of the largest remnant at the end of each simulation are provided in Table D.4.

Notably, there seems to exist a discrete distance where the initially bilobate structure of the aggregate can stay intact. Inside this distance, nearly all particles originally belonging to moon-let A get tidally stripped from the moon, leaving the original moonlet B behind with a few additional grains that are either retained or re-accreted. From the distance of 1.08r₀ to 1.14r₀, only about 10% of moonlet A is shed, while the formed moon retains all of its mass from 1.15r₀ and onward. This phenomenon can be seen in Fig. 10, which shows the a/b values for each moon at the end of the simulations. At post-merger distances of r < 1.08r₀, the moon undergoes a catastrophic tidal disruption event, leaving behind a Dimorphos-like OS on an eccentric orbit, e_orb ≳ 0.4. Once more, we see examples of how our DEEVE estimates for a/b fail for very small bodies, as all these objects exhibit clear OS shapes from a visual evaluation, but only three of them have values matching that of Dimorphos. When instead only 10% of the mass is lost in a mild tidal disruption event, we see a sharp increase in a/b, as the final body highly resembles the reference W24 merger, which evolved at r₀ with the disc still present. Even for the same type of disruption, there are distinct variations in their geometric features, as the moons evolved at 1.08r₀ and 1.09r₀ exhibit LAB features. Further out from the primary, the objects become more compact and homogeneous with a dip in a/b for 1. 11 r₀, before the bilobate structure re-emerges at 1.13r₀. To understand this trend, we take a look at the rotational and orbital periods at 48 h for all moons, shown in Fig. 11. Here, the y-axis shows the periods at the end of each simulation for different initial post-merger distances from the primary. Note that the orbital periods are highly similar for cases that undergo mild and no disruption because the former migrated outwards after losing mass (see Table D.4). With the reduction in a/b in Fig. 10, there is a corresponding increase in rotational period in Fig. 11, indicating that a non-synchronous rotational state leads to a less prolate shape. From Table D.4, we observe that it is mainly changes in b that cause this decrease, as both a and ϕ stay similar compared to other cases in this regime. The changes in b can be attributed both to tidal distortion in the form of stretching and re-accretion of shed material along this principal axis. Analysing the change in mass over time for each tidal disruption scenario in Fig. 12, we also see how the moons undergo mass loss at different points in their dynamical evolution, showing the distance-dependence of the tidal strength. Given that all bodies have the same initial angular velocity, this comes from the fact that each moon is at a different stage of its rotation when the mass loss occurs. This, in turn, is highly dependent on the level of tidal dissipation in the aggregate material, which greatly decreases with distance from the primary. Studying the yaw angle, ψ, (computed as the rotation around the principal axis in the equatorial plane of the primary) at the time step where each body first undergoes mass loss, we see a clear pattern where the bodies at 1.08r₀ and 1.09r₀ have ψ-values closer to zero than for the subsequent cases (see Fig. C.3 in Appendix C). Hence, the tidal forces are only strong enough at distances 1.06r₀ ≤ r ≤ 1.10r₀ to synchronise the rotation of the objects on the timescales considered here (48 h), causing the pull to be parallel to the longest principal axis of the aggregate, a, leaving it elongated and creating a distinct neck between the lobes. In Fig. 13, we show the change in yaw over time for no-disruption scenarios. An initial offset of the libration amplitude is expected, due to the rotational period being kept static while the orbital period changes with distance, but the evolution over time is non-linear, highlighting the importance of tidal dissipation for the orbital and structural evolution of each moon.

At 1.14r₀, the elongation of the moon becomes more pronounced again as tides are strong enough to stretch out the body, while also stripping off some mass. The a/b value reaches its peak at 1.15r₀, where all mass is retained, but tidal forces are still sufficiently strong to elongate the body. Beyond this point, tides gradually grow weaker but also act at larger angles relative to the principal axis, a. The change in libration amplitude with distance over time reinforces the complicated nature of the aforementioned spin-state of each moon and how it depends on tidal disruption history and tidal strengths at the new post-merger positions of the moons. Although the exhibited pattern calls for deeper investigation beyond the scope of this work, the sensitivity of the spin-orbit state points towards the high tidally dissipative nature of rubble-pile satellites indicated from the results of Burnett et al. (2024). Crucially, from this analysis, we learn that the sub-escape-velocity nature of mergers in this regime highly affects the final moon shapes. Due to the soft impact and rigidity of our moonlets, there is limited initial displacement of boulders. Hence, structural failure always occurs at the same places, the innermost edge of moonlet A closest to the primary and the connection point between the two moon-lets. This leads to a predictable behaviour, where the entirety of moonlet A is stripped off at r < 1.08r₀ and only the tip of moon-let A is lost at 1.08r₀ ≤ r < 1.15r₀. In the latter region, the tides also serve to create the aforementioned distinct neck between the two moonlets, most visible in the cases 1.08r₀, 1.09r₀ and 1.13r₀. We further note that distinct necks can also form without the moon undergoing any tidal disruption events, as the case of α = 225° in Figs. 3 and 5, where an OAB structure has been reshaped into a LAB without any mass loss involved.

Fig. 9 Shapes and orientations of aggregates after 45 h of post-merger evolution, based on the case of α = 0° in Fig. 3. Each case represents a different shift in distance from the primary barycentre after the merger has occurred, between r0 and 1.19r0.

Fig. 10 Principal axis ratio a/b for the largest remnant after 45 h of simulation at different shifted distances from the primary barycentre post merger (in r0 = 2.43Rprimary) for the case of α = 0°. The dashed line shows the estimated a/b value of Dimorphos (Daly et al. 2024), while the dashed-dotted line shows the final a/b value for the reference W24 merger, evolved with the debris disc still present.

Fig. 11 Orbital (red square) and rotational (black dot) periods for the largest remnant after 45 h of simulation (48 h total) when initiated at different shifted distances from the primary barycentre for the case of α = 0°.

Fig. 12 Change in mass over time for the simulation scenarios leading to the shapes in Fig. 9 that undergo tidal disruption. The bodies undergoing catastrophic disruption have dashed lines, while the mild disruption scenarios have solid lines. The colour of each line represents its initial post-merger distance from the primary.

Fig. 13 Post-merger change in yaw angle over time for the no-disruption cases of Fig. 10. The colour of each line represents the initial distance in r0 = 2.43Rprimary. The data has been smoothed using a cumulative moving average with a window width of 30 time steps to reduce noise caused by small shifts in the body-fixed frame.

4.2 Particle number and size distribution

With the considerable improvement in resolution achieved with the GPU-acceleration scheme of GRAINS, increasing the particle number by almost a factor of 20, we could more closely study the effect of having a steeper SFD, leading to a considerably larger fraction of smaller particles compared to the relatively low-resolution configuration of W24. At a first glance, coarser, lower-resolution distributions offer an approximation of deformation, as we see clear similarities when comparing the shapes in Figs. 3 and 6. Moreover, it suffers comparable levels of disruption to our high-resolution moons, where most of moonlets A and A2 are lost in the respective cases. However, the increase in stability where the largest remnant of the high-resolution aggregate retains more mass shows the importance of resolving particle-particle contacts on a smaller scale when studying tidal deformation.

Recent studies on the behaviour of granular media in tidal disruption scenarios primarily investigate the effects of altering particle distribution properties for rubble piles consisting of spherical elements (Zhang & Lin 2020; Zhang & Michel 2020), approximating granular contact mechanisms with a linear spring-dashpot system along with a non-spherical shape parameter (Schwartz et al. 2012; Zhang et al. 2017). The few studies that use polyhedral particles either opt for a monodisperse (Marohnic et al. 2023) or uniform SFD (Movshovitz et al. 2012). It is difficult to meaningfully compare our results with the aforementioned DEM studies because of the spread in: problem scale, ranging from encounters with stars and terrestrial planets to small asteroids, the variation in particle size distributions, shapes and packing, as well as the contact methods and the numerical capabilities available at the time each study was conducted. Even so, we address a few points and key differences that could prove useful for future investigations of rubble-pile stability. Firstly, while brought up as likely to be important, studying the effect of particle SFD for rigidity was beyond the scope of these works. Since then, a series of articles using SPH methods find that boulder packing and size distribution can have a crucial effect on the stability of rubble-pile objects during high-energy collisions (Raducan et al. 2022; Raducan et al. 2024; Raducan et al. 2024). Moreover, for our purposes, observations of Dimorphos’s surface indicate a boulder SFD with a mean and maximum diameter of 1.4 and 16 m, respectively, with a best fit power-law index of −3.4 ± 1.3 aligning with similar observed SFDs from Didymos, Itokawa, Eros, and Toutatis that all have indices below −3 (Pajola et al. 2024). Secondly, we once more emphasise that rubble-pile satellites in binary asteroid systems are likely more dissipative than previously believed (Burnett et al. 2024) and that polydisperse SFDs can facilitate propagation of mechanical waves (Marti et al. 2024) in granular media and reduce macroporosity due to improved packing (Ferrari & Tanga 2020). These conclusions align with our results, as there is an evident change in the amount of deformation that occurs during mergers for our high-resolution case. Thirdly, we highlight the fact that the mentioned studies have been focusing on disruptions caused by close encounters with stars or planets. It might be that the difference in scale explains why they detected a seemingly weaker dependence on resolution for aggregate stability. The consequent improvement in stability caused by increasing N_p from 10³ to 10⁴ is possibly too weak relative to the strength of the tidal forces in these environments, while highly relevant when the primary is a small asteroid. A difference in mass loss of 0.05M₀ when comparing the low- and high-resolution disruptions is small in the scope of the system. Yet, this amount of mass can be decisive when attempting to accurately capture deformation. Finally, it is critical to acknowledge the uniqueness of the structure considered in our tidal study. Not only is it highly elongated, but it is also of a bilobate nature with a weak neck. Investigating a more structurally homogeneous target, we would likely not have identified the notable shift in dynamical behaviour under the influence of tidal forces.

Thus, we argue that it is crucial that DEM studies use models with polydisperse particle distributions and methods that can capture the mechanical effects of irregularly shaped grains to get a full picture of the reshaping of rubble-pile targets in asteroid environments from mergers and tidal forces. With all the complexities at hand, this calls for a deeper investigation beyond the scope of this work into the effects of varying particle SFDs, aggregate bulk density and packing.

4.3 Sub-escape-velocity merger formation scenarios for Dimorphos and Selam

While we still have limited direct observations of asteroid satellites, the recent imaging of Dimorphos by the DART spacecraft (Daly et al. 2023; Chabot et al. 2024) and Selam by Lucy (Levison et al. 2024) has provided us with invaluable information. Given that they both orbit small, rapidly rotating top-shaped asteroids, theirexistence agrees with a debris disc formation scenario. From the results and analysis in W24, there were strong hints that sub-escape-velocity mergers can generate both oblate and bilobate shapes. This has been further reinforced in this work. Additionally, we have shown that OS and Selam-like LAB satellites do not necessarily need to form directly from mergers, as they can end up with similarshapes from post-mergerdynamical events. Dimorphos-like objects can consistently form from catastrophic tidal disruptions, as long as the merger occurs sufficiently close to the primary. In our simulations, these objects end up on eccentric orbits. However, given the influence of a disc on the orbital evolution of our moons (see Fig. 2), we cannot exclude the possibility that the presence of a disc can circularise the orbit. Moreover, from Fig. 11 we see that the final rotation of the largest remnant in these catastrophic disruption events is dependent on whether it occurred when its progenitor was tidally locked with the primary. From measurements of Dimorphos’s rotational state before the DART impact, it appeared to be in a synchronous state (Rivkin et al. 2021; Agrusa et al. 2022; Naidu et al. 2024). For example, for the case of 1.05r₀, the rotational period of the body appears to match the corresponding period for a circular orbit at approximately 3Rprimary, which is where Dimorphos is situated in the Didymos system. Hence, disregarding any secular synchronisation of the rotation for a moment, given that the tidal disruption distance would be dependent on the progenitor’s bulk density, we should be able to predict where a possible merger and subsequent disruption must have occurred for the satellite to have formed via this scenario. This would then be made possible with improved measurements of Dimorphos’s internal properties from the Hera mission (Michel et al. 2022). Numerical experiments of impacts into a Dimorphos-like target further find that the bulk density of the object is likely lower than that of Didymos (Raducan et al. 2024), which also speaks for this tidal disruption formation model, as this would shift the disruption limit outwards, further from the primary.

Predicting the formation history of Selam is a more challenging task. The contact binary satellite is surprisingly massive with an estimated m_moon = 0.06 M_primary and orbits Dinkinesh with a semi-major axis of 9 Rprimary. What is more, the geological features are highly distinct, with a narrow neck connecting two equal-sized lobes with diameters of 210 m and 230 m (Levison et al. 2024). One of the lobes has a noticeable ridge that appears to span around it, inclined at an angle of approximately 50^° relative to the orbital plane. It was suggested by Levison et al. (2024) that this ridge could have been formed by accretion of material from a debris disc. Alternatively, it could have formed by mergers of similar-sized moonlets (Raducan et al. 2025). In order to produce Selam from a single mass-shedding event, the debris disc must be of substantially higher mass compared to the cases explored in W24 and Agrusa et al. (2024). Its contact binary structure could then be formed via a mild tidal disruption after most of the disc mass has been bound in a single, merged moon similar to our α = 0° case. From our analysis in Sect. 4.1, we know that such a bilobate structure would fail at two specific points, the closest point to the primary and the contact area between the moonlets. Hence, there could be a discrete distance, such as 1.08r₀ in Fig. 9 for the α = 0° body, at which the tidal forces will strip off the innermost edge of the most massive moonlet, leaving it less prolate while at the same time pulling the lobes apart and creating the neck between them. Such an event is likely to be followed by migration further out into the system, where the satellite would be less susceptible to homogenisation by tides, retaining its LAB shape. It remains to be investigated whether a disruption followed by circularisation via angular momentum exchange with a surrounding disc could be a key mechanism in placing Selam at its current position at 9 R_primary, or if higher order effects such as binary YORP are necessary for it to migrate this far out (Jacobson et al. 2014). Such a migration would occur over a period of approximately 10⁶−10⁷ years (Merrill et al. 2024). Still, within this single mass-shedding formation scenario, it would be non-trivial to explain the characteristic lobe ridge. We cannot resolve such a feature in our current debris disc simulations, as that would require more particles. To reach impact velocities that could create it via a merger, the impact angle would have to be much larger than those observed in W24. For this to occur, the largest moonlet would have to undergo a tidal disruption, stripping off most of its mass and leaving it less prolate while putting it on a highly eccentric or unbound orbit where it collides with a lower-mass moonlet (which tend to have more oblate structures). It could alternatively be accelerated by some other dynamical effect, such as Lindblad resonance with the remaining disc, as proposed by Raducan et al. (2025). Given that the discs in W24 and Agrusa et al. (2024) fully disperse over just a few days, we argue that this would be an unlikely mechanism in our model. For a formation history involving two or more mass-shedding events, we can imagine that two separate OS satellites were generated via tidal disruptions, one-by-one, leading to migration and a merger beyond 3R_primary with a velocity near υ_esc, similar to the pyramidal regime scenario discussed in both Levison et al. (2024) and Raducan et al. (2025). Whichever scenario is true for the mysterious contact binary satellite, it is evident that its dynamical history is highly complex and bound to teach us more about the chaotic nature of binary asteroid system formation.

5 Conclusions

In a first approach to explore which dynamical mechanisms determine moon shapes resulting from sub-escape-velocity mergers, we have performed several sets of polyhedral N-body simulations. We find that each moonlet undergoes little initial deformation due to the impact because of its sub-escape-velocity nature. The orientation of the prolate and most massive moon-let in our configuration highly affects the final shape of the formed moon, creating a wide range of post-merger structures that we classified as long-axis bilobate (LAB), off-axis bilobate (OAB), and short-axis bilobate (SAB). Further, we studied the stability of the formed moons, letting them dynamically evolve around their Ryugu-like primary for up to 48 h at various distances from the primary. Depending on their shape and rotational state, the bodies undergo different types of deformation via tidal disruption as well as tidal distortion that can significantly alter their structures. As a result, some of the moons become more homogeneous and lose their bilobate nature, while it is retained or even emphasised in other cases. The structurally fragile LABs sometimes undergo catastrophic tidal disruptions, becoming stripped of more than half of their mass and leaving behind an oblate spheroid similar in shape to Dimorphos. For the fiducial merger case, we observed that the severity of the mass loss occurs in discrete regimes of distance, where the innermost lobe is always fully disrupted closer-in, while only around 10% of the mass is lost in the next regime. Here, the tidal pull is generally directed along the longest principal axis of the body, leading to an accentuated neck between its two lobes. We believe this mechanism to be a potential explanation for the possible contact binary shape of Selam. Further out, the entire structure stays intact and any deformation of the body depends on its initially non-synchronous spin-orbit state. Employing a higher-resolution SFD, which is more true to the observed boulder diameters on rubble-pile asteroid surfaces, leads to more initial deformation from the merger, as having more small boulders in the aggregate smooths out contacts between the larger constituents. Nevertheless, the initial orientation of the most massive moonlet still highly affects the final shape of the merged object, and we observed similarities between the low-and high-resolution outcomes.

With these results in mind, along with the initially elongated shapes of the largest moonlets born in a binary asteroid formation scenario, we argue that it is necessary to take the non-sphericity of moonlets into account when modelling mergers in asteroid debris discs. Moreover, the use of irregularly shaped particles and polydisperse SFDs reproduces key granular mechanics that are crucial when capturing dynamical effects such as tidal distortion and disruption. Not only do these processes alter the shape of a rubble-pile satellite, but the timing of a tidal disruption event also highly affects its resulting orbit, determining how far it migrates outwards and its spin-orbit state. As shown in this work, even a shift of 1% in distance from the primary can lead to significant changes in the structural evolution of these objects. Furthermore, the change in dynamical behaviour for a merged satellite in the presence of a debris disc and the potential sub-sequent accretion of mass adds another layer of complexity to this problem that must be studied to properly understand the origin of atypical asteroid satellite shapes. Even so, with more data on the internal structure, particle size distribution, and bulk density of Dimorphos from the Hera mission (Michel et al. 2022) and future missions involving in situ observations of asteroid satellites such as Selam along with more in-depth numerical studies of the debris disc model and tidal evolution of rubble-pile satellites, we will be able to further constrain their formation scenarios and dynamical history and obtain a more satisfying answer to how these fascinating objects came to be.

Data availability

The movies are available at https://www.aanda.org

Acknowledgements

J.W. and F.F. acknowledge funding from the Swiss National Science Foundation (SNSF) Ambizione grant no. 193346. M.J. acknowledges support from SNSF project no. 200021_207359. We extend our gratitude to Harrison F. Agrusa, Sabina D. Raducan, Po-Yen Liu and Richard Cannon for helpful discussions and input throughout this project.

References

Appendix A Rotations and quaternions in GRAINS

In Chrono, the orientation of each particle is tracked using unit quaternions, which are four-dimensional vectors that have the following form: $$\mathit{q}=q_0+q_1\mathit{i}+q_2\mathit{j}+q_3\mathit{k},$$(A.1)

where q₀, q₁, q₂, q₃ are real numbers and i, j, k are basis elements fulfilling i² = j² = k² = i jk = −1. Moreover, they always satisfy |q| = 1. As a result, each particle’s rotational state can conveniently be represented by just four numbers, with q = {1, 0, 0, 0} being the initial state. Another notation q = (s, u) is based on the fact that a quaternion consists of a scalar s and an imaginary vector u.

Assuming we want to rotate a given particle by an angle α around its centre-of-mass in an x, y, z Cartesian coordinate system, this can be achieved using quaternion multiplication. This is possible due to the fact that the product of a quaternion q by its conjugate q* = (q₀ − q₁i − q₂ j − q₃k), is a quaternion without its imaginary part as $\mathit{q}{\mathit{q}}^{*}=\left(q_0^2+q_1^2+q_2^2+q_3^2\right)$. Hence, it follows that a quaternion describing the rotational state of particle b, p_b, can be transformed by another quaternion, q, representing some rotation. Given a rotation α around some general unit vector n, this quaternion has the form q = {cos(α/2), n_x sin(α/2), n_y sin(α/2), n_z sin(α/2)}. This operation leads to a new transformed quaternion: $${{\mathit{p}}^{\prime }}_b=\mathit{q}{\mathit{p}}_b{\mathit{q}}^{*}.$$(A.2)

Since Chrono uses unit quaternions, the norm of p_b remains unchanged. Furthermore, this treatment can also be applied to vectors representing position, velocity and angular velocity, which follows from the fact that some vector υ can be represented by the quaternion p_υ = (0, υ). Applying Eq. (A.2) then yields the transformed entity ${{\mathit{p}}^{\prime }}_{\upsilon }=\left(0,{\mathit{\upsilon }}^{\prime }\right)$. For more details regarding quaternion properties and their use in Chrono, see Tasora (2020).

Appendix B Bilobate DEEVE fitting

In order to determine the DEEVE of each lobe after the merger and estimate their semi-major axes, we implemented a simple algorithm. First, we assumed that each particle belonged to the same parent moonlet pre- and post-merger and computed the barycentre of each lobe, s_A and s_B. Second, we iterated over each particle and evaluated their position relative each barycentre regardless of which lobe is their initial parent to account for mixing. For a given particle i, if the distance d_A = |r_i − s_A| was less than d_B = |r_i − s_B|, it belonged to lobe A after the merger and vice versa. After repeating this process for all particles, we calculated the new barycentre and total moment of inertia for each lobe. Given the soft impact and the little-to-no mixing, one iteration of this algorithm was enough to give good visual DEEVE fits.

Appendix C Complementary plots

First, we show the resulting shapes of the merged moons from Fig. 3 after 45 h of additional simulation (48 h in total) in Fig. C.1. Note that most of the bodies have suffered catastrophic tidal disruption, with α = 300° being in the middle of its disruption at the end of the simulation. We also include a plot depicting the corresponding shapes for a case where they were evolved at 1.2r₀ after 48 h of total simulation in Fig. C.2.

The next plot in Fig. C.3 shows the yaw angle for the distance dependency study in Sect. 4.1 when a body first undergoes mass loss. A majority of the bodies with synchronous rotation in Fig. 11 were also tidally locked when the disruption event occurred. For the bodies with r > 1.14r₀, we simply plotted the yaw angle at the end of the simulation, which explains the variance in rotational period for these objects.

Fig. C.1 Shapes of the bodies in Fig. 3 after each system has been evolved for an additional 45 h (48 h in total).

Fig. C.2 Same as Fig. 3 but after putting each body on a circular orbit with a semi-major axis of 1.2r0 and evolving each system for and additional 45 h (48 h in total), where r0 is the distance from the primary at 3 h.

Fig. C.3 Yaw angle, ψ, at the first time step after mass loss for each case in Fig. 10. Note that the no-disruption cases show the yaw angle at the final time step.

Appendix D Moon property tables

All Tables

All Figures

Fig. 1 Visualisation of how new initial conditions are generated by rotating the most massive moonlet counterclockwise around its barycentre by an angle, α, here set to 45°. On the left is the original configuration from the W24 merger, with a vector â0 parallel to the aggregate’s semi-major axis. We performed the rotation for each particle, ending up with a new aggregate with its semi-major axis parallel to the vector â1.

In the text

	Fig. 2 Evolution of semi-major axis (top), eccentricity (middle), and inclination (bottom) for the moons formed in the W24 merger with (solid blue) and without a disc present (dashed orange).
In the text

Fig. 3 Shapes and orientations of each formed moon after 3 h of simulation time after a sub-escape-velocity merger between the two moonlets in Table 2. The value α, represents the counter-clockwise rotation relative to the initial position of moonlet A in Fig. 1 around its barycentre. The particles originally belonging to moonlet A have been coloured in blue, while the corresponding particles from moonlet B are orange. The black arrow in each plot is pointing in the direction of the primary.

In the text

	Fig. 4 Angle, γ, between the longest principal axis of the DEEVE of moonlet A (blue); a; and the vector rAB, which is the positional vector between the barycentres of moonlet A and moonlet B (orange).
In the text

	Fig. 5 Same as Fig. 3 but after putting each body on a circular orbit with a semi-major axis of 1.1 r0 and evolving the system for an additional 45 h (48 h in total), where r0 is the distance from the primary barycentre at 3 h.
In the text

	Fig. 6 Higher-resolution versions of a few of the cases given in Fig. 3, with the moonlets consisting of 35 729 particles in total. The moonlet properties are given in Table 4.
In the text

	Fig. 7 Final shapes for the case of α = 0° in Fig. 6 after an additional 21 h of simulation (25 h in total) for a bulk density of 626.8 kg m−3 (left) and 1190 kg m−3 (right).
In the text

	Fig. 8 Impact velocity in terms of mutual escape velocities for the case of the merger data in Table 1 and moonlets in Table 2, setting a static bulk density of 1190 kg m−3. The red square represents the values for the W24 merger.
In the text

	Fig. 9 Shapes and orientations of aggregates after 45 h of post-merger evolution, based on the case of α = 0° in Fig. 3. Each case represents a different shift in distance from the primary barycentre after the merger has occurred, between r0 and 1.19r0.
In the text

Fig. 10 Principal axis ratio a/b for the largest remnant after 45 h of simulation at different shifted distances from the primary barycentre post merger (in r0 = 2.43Rprimary) for the case of α = 0°. The dashed line shows the estimated a/b value of Dimorphos (Daly et al. 2024), while the dashed-dotted line shows the final a/b value for the reference W24 merger, evolved with the debris disc still present.

In the text

	Fig. 11 Orbital (red square) and rotational (black dot) periods for the largest remnant after 45 h of simulation (48 h total) when initiated at different shifted distances from the primary barycentre for the case of α = 0°.
In the text

	Fig. 12 Change in mass over time for the simulation scenarios leading to the shapes in Fig. 9 that undergo tidal disruption. The bodies undergoing catastrophic disruption have dashed lines, while the mild disruption scenarios have solid lines. The colour of each line represents its initial post-merger distance from the primary.
In the text

	Fig. 13 Post-merger change in yaw angle over time for the no-disruption cases of Fig. 10. The colour of each line represents the initial distance in r0 = 2.43Rprimary. The data has been smoothed using a cumulative moving average with a window width of 30 time steps to reduce noise caused by small shifts in the body-fixed frame.
In the text

	Fig. C.1 Shapes of the bodies in Fig. 3 after each system has been evolved for an additional 45 h (48 h in total).
In the text

	Fig. C.2 Same as Fig. 3 but after putting each body on a circular orbit with a semi-major axis of 1.2r0 and evolving each system for and additional 45 h (48 h in total), where r0 is the distance from the primary at 3 h.
In the text

	Fig. C.3 Yaw angle, ψ, at the first time step after mass loss for each case in Fig. 10. Note that the no-disruption cases show the yaw angle at the final time step.
In the text