Rapid formation of binary asteroid systems post rotational failure: a recipe for making atypically shaped satellites

Binary asteroid formation is a highly complex process, which has been highlighted with recent observations of satellites with unexpected shapes, such as the oblate Dimorphos by the NASA DART mission and the contact binary Selam by NASA’s Lucy mission. There is no clear consensus on which dynamical mechanisms determine the final shape of these objects. In turn, we explore a formation pathway where spin-up and rotational failure of a rubble pile asteroid lead to mass-shedding and a widecircumasteroidaldebrisdiskinwhichthesatelliteforms.Usingacombinationofsmooth-particle hydrodynamical and N -body simulations, we study the dynamical evolution in detail. We find that a debris disk containing matter corresponding to a few percent of the primary asteroid mass extending beyond the fluid Roche limit can consistently form both oblate and bilobate satellites via a series of tidal encounters with the primary body and mergers with other gravitational aggregates. Principally, satellites end up prolate (elongated) and on synchronous orbits, accreting mainly in a radial direction while tides from the primary asteroid keep the shape intact. However, close encounters and mergers can break the orbital state, leading to orbital migration and deformation. Satellite–satellite impacts occurring in this regime have lower impact velocities than merger-driven moon formation in e.g. planetary rings, leading to soft impacts between differently sized, non-spherical bodies. The resulting post-merger shape of the satellite is highly dependent on the impact geometry. Only moons having experienced a prior mild or catastrophic tidal disruption during a close encounter with the primary asteroid can become oblate spheroids, which is consistent with the predominantly prolate observed population of binary asteroid satellites.


Introduction
The presence of atypically shaped satellites in binary asteroid systems has been an trending topic as of late (Madeira and Charnoz, 2024;Agrusa et al., 2024).First, close-up observations were made of the spherically oblate Dimorphos orbiting its main body Didymos by NASA's DART (Double Asteroid Redirection Test) spacecraft, which then proceeded to impact the body on September 26th 2022, altering its orbit and likely its shape in a historical event (Daly et al., 2023a;Chabot et al., 2024;Rivkin et al., 2021;Raducan et al., 2024a).More recently, images of a contact binary satellite, now named Selam, were provided by NASA's Lucy mission during a flyby of the asteroid Dinkinesh (Levison et al., 2021(Levison et al., , 2024)).Observations report that ∼ 15% of the near-Earth asteroid (NEA) population have satellites, a number that increases to ∼ 65% for rapidly rotating bodies with a diameter of at least ∼ 300 m (Margot et al., 2002;Pravec et al., 2006).A large majority of the observed satellites are prolate (elongated) in form (Pravec et al., 2016), making the shapes of both the oblate Dimorphos and bilobate Selam surprising.Notably, while significant parts of these bodies have yet to be imaged, which adds a level of uncertainty to their inferred shapes, lightcurve measurements of asteroid systems are biased against detecting satellites that are more oblate (Pravec et al., 2016;Agrusa et al., 2024).Consequently, oblate spheroidal satellites could be more prominent than expected from current observational data.
The primary bodies in both systems are small, with Didymos having a volume-equivalent diameter of ∼ 760 m (Barnouin et al., 2023) and Dinkinesh a similar effective diameter of ∼ 720 m (Levison et al., 2024).Moreover, the two display the typical spinning-top-shape of rapidly rotating asteroids (Pravec et al., 2006;Pravec and Harris, 2007;Benner et al., 2015).The high rotation rate of these objects might be primarily caused by the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect, where a weak net torque acts on asteroids due to anisotropic re-emission of absorbed sunlight over long periods of time caused by their irregular shape (Rubincam, 2000).Many of these asteroids are so-called 'rubble piles' and are gravitational aggregates consisting of a large number of smaller constituent particles with littleto-no bulk tensile strength and moderate porosity.These objects have observed spin-rates exceeding the theoretical disruption limit for fluid materials, indicating that they have shear strengths consistent with angles of friction around 40 • , which yields behaviours similar to terrestrial granular materials (Scheeres et al., 2010;Sánchez andScheeres, 2014, 2016;Walsh, 2018).
The observed spin-rates of asteroids with diameters larger than 200 m and smaller than 40 km show that very few bodies are rotating with periods shorter than ∼ 2.2 h, indicative of a critical spin limit (Pravec and Harris, 2000;Pravec et al., 2002).In turn, there must exist some mechanism to prevent these bodies from spinning faster.It was postulated early on after the identification of this critical rotational period that YORP could be spinning rubble pile asteroids up beyond this limit, leading to rotational failure and disruption, which would in turn reduce the rotational angular momentum of the body and decrease its rotation rate (Rubincam, 2000).This mechanism was then further tied to the possible formation of asteroid binaries (Vokrouhlický and Čapek, 2002;Bottke et al., 2002Bottke et al., , 2006;;Walsh et al., 2008;Harris et al., 2009).
In this paper, we explore a formation scenario for binary asteroid systems where spin-up of an initially spheroidal rubble pile asteroid leads to a chaotic, instantaneous massshedding event, generating a wide debris disk and explore how the dynamical evolution of the disk post rotational failure can produce asteroid satellites.In particular, we focus on the various dynamical mechanisms that govern the shapes of moons in this regime and constrain the different formation pathways, both qualitatively and quantitatively, that can produce the atypical shapes of Dimorphos and Selam.In order to do so, we employ a combination of spin-up simulations of rapidly rotating rubble pile asteroids performed with a smooth-particle hydrodynamics code that provides initial conditions for a gravitational N-body code with non-spherically shaped particles that evolves the disk until we obtain a satellite of significant mass.In section 1.1, we introduce and summarise the previous work on binary asteroid formation.In section 2, we go through the two different numerical codes used in our simulations and the interface we have created to connect them.In section 3 we show the results from our main model, based on a Ryugu-like body with similar properties and in section 4, we apply our model to the Didymos system and compare how the outcome changes with a smaller primary body.We then discuss the effects of disk geometry, moonlet merger properties and particle size frequency distribution in section 5, which is followed by conclusions in section 6.

Binary asteroid formation
While the narrative of spin-up via YORP leading to structural failure of rubble pile surfaces and global landslides is heavily supported by both theoretical and numerical work along with observations, there is yet no consensus in the community regarding just how these mass-shedding events lead to the formation of binary asteroid systems, especially when considering the existence of atypically shaped satellites such as Dimorphos or Selam.The timescales for the formation process are poorly constrained, where e.g.surface analyses of Didymos and Dimorphos indicate that the main body is between 40 and 130 times older than its satellite (Barnouin et al., 2023).In turn, it is difficult to determine how much mass needs to be shed over how long a period to successfully form a stable satellite.

Continuous mass-shedding
A first attempt to explicitly tie equatorial mass-shedding to binary formation was made by Walsh et al. (2008), who used a numerical hard-sphere discrete element model (HS-DEM) to perform N-body simulations of the slow shedding of mass from a cohesionless, rapidly rotating primary spun up by YORP.The primary itself was modelled using one thousand rigid, similarly sized spheres arranged into a spheroidal aggregate using hexagonal closest packing (HCP).This particle configuration yielded an angle of friction of around 40 • , providing the shear strength needed to withstand internal deformation until surface shedding could occur.As the rotation rate of the main body was gradually increased to imitate the process of slow spin-up via the weak torques from YORP and approached the critical spin limit, Walsh et al. observed how particles began moving from the poles towards the equator which came to serve as the origin of the subsequent mass loss.In turn, the primary became oblate and obtained the characteristic top-shape observed for so many rapidly rotating asteroids.Mass was then ejected continuously through shedding events from the equator of the primary and eventually, with enough mass in orbit, the debris began quickly accumulating into a satellite.A followup study investigating a broader set of initial conditions with higher resolution was done in Walsh et al. (2012), reaffirming the correlation between this formation pathway and the observed population of binary asteroid systems.
While the formed satellites have orbits and masses that align with the observed population, the two studies did not explore their final shapes.Due to the nature of tidal forces for granular aggregates, most of bodies created in this manner must first begin accumulating near the theoretical fluid Roche limit, as most bodies formed interior of this distance would be tidally disrupted before reaching a high enough mass (Holsapple and Michel, 2006;Sharma, 2009).Initially having an orbit near the Roche limit, the tidal forces acting on the aggregate should render it prolate (Porco et al., 2007;Tiscareno et al., 2013) and in synchronous rotation, which aligns with observations of binary satellites on close orbits (Pravec et al., 2016).As we will show in this work, a rubble pile moon formed without undergoing any events of deformation does indeed exhibit this pattern, which was also indicated by other recent studies from Madeira and Charnoz (2024) and Agrusa et al. (2024).
Hence, the long-term formation process explored in Walsh et al. (2008Walsh et al. ( , 2012) ) can account for the predominantly prolate population, but offers no satisfying explanation for the oblate shape of Dimorphos, which has estimated semimajor axes of ≈ 87.90×86.96×57.16m (Daly et al., 2023b).The body further exhibits geological properties of a rubble pile (Barnouin et al., 2023), which aligns with results from impact simulations (Raducan et al., 2024a,b), meaning it cannot simply be a large monolithic boulder from a past, chaotic dynamical event.A detailed review of the dynamical state of Dimorphos before and after the DART impact can be found in Richardson et al. (2024).To obtain further data on the post-impact state and the geological properties of Dimorphos, the European Space Agency's Hera mission is sending a spacecraft which will arrive in 2027 to perform a close-up analysis of both Didymos and Dimorphos, as well as in situ measurements of the latter with CubeSats (Michel et al., 2022).This will greatly aid the endeavour to constrain the internal structure of Dimorphos and in turn its dynamical history.
As previously mentioned, polydisperse aggregates react differently to external forces and strain than those consisting of similarly sized particles.Moreover, the packing of material and their effective shapes can greatly affect the dynamics of granular aggregates, which was shown by e.g Zhang and Michel (2020), who investigated tidal disruption of rubble piles with different packing, particle shapes and numerical treatments for particle-particle contacts.More specifically, they found that models based on HSDEM and HCP packing produce different dynamical results than the more complex SSDEM treatment with random closest packing.This comes from HSDEM being unfit to simulate granular systems where the particles sustain long-lasting contacts (e.g.Sánchez and Scheeres, 2016).Furthermore, as pointed out by Agrusa et al. (2024), new inhomogeneities caused by global deformation of the surface of spinning asteroids can highly affect the net torques caused by YORP throughout the satellite formation process, causing large variations in the amount of mass that is being shed over long periods of time (Statler, 2009;Cotto-Figueroa et al., 2015;Zhou et al., 2022).

Contact binary fission
Another formation scenario was explored by Jacobson and Scheeres (2011a), based on the 'rotational fission' of a contact binary asteroid due to long-term spin-up via YORP, where the bilobate object is disrupted and separated into two components and the larger one is the primary while the smaller is the secondary (Scheeres, 2007(Scheeres, , 2009)).This yields a chaotic and instantaneous shedding event where mass ejected into the system is less prone to escape before being accumulated into a stable satellite, which is another issue that long-term YORP-driven formation scenarios must take into account (Jacobson and Scheeres, 2011a).Under the assumption that the contact binary is initially arranged to satisfy its global minimum energy configuration, they generated each contact binary body as two fused ellipsoids with their longest axes aligned.They also used a second type of configuration, where the primary is an ellipsoid while the 'secondary' consists of two spheroids, allowing for a 'secondary fission' event where the created satellite also can undergo fission due to spin-orbit coupling, leading to a ternary (triple) system.The introduction of this secondary event allowed for the fissioned body to dissipate energy in an alternative way for systems with a low mass ratio between the secondary and primary body (<0.2), resulting in positive free energy.When the secondary body was kept monolithic for these initial conditions, it would then ultimately always escape the system.
Because of the implementation approximating the granular aggregates as consisting of two or three components, they were able to perform very long-term dynamical evolution (up to 1000 years) of the systems post rotational fission using a high-order integration scheme, accounting for inelastic collisions, the non-spherical gravitational potentials of each component and energy dissipation due to mutual body tides using a semi-analytical implementation.Nevertheless, despite the advanced treatment of the system dynamics capturing influential long-term dynamical effects, the choice of using rigid bodies as representations for each component may not capture the full complexity of granular mechanics.As we have already established, tidal forces highly affect rubble piles forming within the theoretical fluid Roche limit, causing deformation, mass-shedding or total disruption and the resilience of a rubble pile to these effects depends strongly on internal packing, cohesion, as well as the size frequency distribution and shape of its particles (Zhang and Michel, 2020).For a granular aggregate to survive within the theoretical tidal disruption limit, it would be necessary for its material to have a friction angle larger than zero and particle interlocking effects (Movshovitz et al., 2012), significantly higher bulk density than the resulting primary or some level of cohesion (or perhaps a combination of all these effects), providing the gravitational aggregate with higher internal mechanical strength such that it can withstand disruption (e.g.Holsapple and Michel, 2008;Sánchez and Scheeres, 2016).For example, as shown by Agrusa et al. (2022), Dimorphos would be tidally disrupted if it formed from a rotational fission event with Didymos, which makes such a formation scenario implausible given its seemingly low level of cohesion (Raducan et al., 2024a,b) and a bulk density similar to that of Didymos (Daly et al., 2023b).
The contact binary fission scenario has since its introduction by Jacobson and Scheeres (2011a) been further explored in several studies (Jacobson et al., 2016;Boldrin et al., 2016;Davis and Scheeres, 2020;Ho et al., 2022).Moreover, a similar mechanism for binary asteroid system formation was tied to the existence of equatorial cavities on rapidly rotating asteroids, where Tardivel et al. (2018) suggested that the ejection of a large chunk of material could lead to the formation of a satellite, assuming the material would have some level of cohesion which would prevent it from being tidally disrupted immediately after the fission event.

Debris disks
More recently, a third possible formation scenario has been suggested that combines the two previous ideas.As mentioned, several studies on spin-up of rubble piles have found that the presence of cohesion and/or friction caused by the angle of friction of the granular material can lead to higher resistance to rotational failure (Hirabayashi et al., 2015;Sánchez andScheeres, 2016, 2020;Zhang et al., 2017Zhang et al., , 2018Zhang et al., , 2021;;Tardivel et al., 2018;Sugiura et al., 2021).In turn, as aggregates with higher mechanical strength can spin up to more rapid rates prior to mass shedding, the ejected material will have higher momentum and travel farther out into the system onto wider orbits, generating a debris disk rather than just an equatorial ridge where particles can accumulate into a satellite over time.In turn, the slow massshedding from movement of surface particles towards the equator during rotational failure observed in Walsh et al. (2008Walsh et al. ( , 2012) ) is linked with the instantaneous mass-shedding of Jacobson and Scheeres (2011a), overcoming the main issues with these theories, namely the unpredictable nature of YORP and escaping ejecta, along with the tidal disruption for the contact binary fission scenario, preventing early formation of satellites.
The instantaneous mass-shedding was first explicitly tied to the formation of asteroid satellites in Hyodo and Sugiura (2022), who dynamically evolved a debris disk around a topshaped primary using three-dimensional smoothed-particle hydrodynamics (SPH) simulations.The initial conditions for their simulations were taken from the study of rotational failure by Sugiura et al. (2021), who also used SPH to study the deformation of rapidly rotating rubble piles.This type of numerical hydrodynamical method has been used extensively to investigate disruption and deformation of rubble pile bodies due to high-energy impacts (e.g.Benz and Asphaug, 1999;Jutzi et al., 2008Jutzi et al., , 2022;;Jutzi, 2014Jutzi, , 2015;;Raducan et al., 2022;Sugiura et al., 2018).In their simulations, Sugiura et al. applied a slow angular acceleration (10 −10 rad s −2 ) to a spherical, Ryugu-like body consisting of ∼ 25, 000 similarly sized particles with random packing.Using a new treatment where they combine cohesion and friction angle into one single parameter in SPH (introduced in section 2.1), they found that rotational failure for bodies with some level of cohesive strength in their model leads to debris disks with masses around 10% of the primary mass after ∼ 30 h of simulation, with their outer edge inside the theoretical fluid Roche limit.Hyodo and Sugiura (2022) then used the same model and numerical implementation as in Sugiura et al. (2021) and proceeded to propagate the system with debris disks with masses of 10 − 20% of the primary for ∼ 700 h, studying the formation of moons in the debris disk.By the end of their simulations, they found that several moons had formed with a combined mass of ∼ 4% of the primary.This is significantly more massive than Dimorphos and satellites of similar systems which are ∼ 1% of the primary mass (Pravec et al., 2016(Pravec et al., , 2019)).In this work, we use a similar SPH-based model (see section 2.1) to investigate the spin-up and subsequent mass-shedding of rubble pile asteroid, adding a more detailed investigation of the subsequent disk evolution, tracking the shape and deformation of formed satellites over time.
In an accompanying paper to Hyodo and Sugiura (2022), Madeira et al. (2023) used 1D hydrodynamical simulations with the code HYDRORINGS (created to investigate the formation of satellites in the rings of Saturn in Charnoz et al., 2010;Salmon et al., 2010) to study the viscous evolution of a debris ring around Didymos and how it can be tied to the formation of Dimorphos.In this study, the authors started their simulations directly when a disk has formed, or at the onset of mass shedding with an analytic prescription for depositing material into the disk.At the start of their simulations, the disk was always confined within the fluid Roche limit and its material then migrated over time due to viscous spreading, taking at least a year to reach this distance from the primary's surface.Any mass moving beyond the fluid Roche limit was assumed to be a satellite and turned into a spheroidal body that was subsequently tracked using an analytical treatment for its orbit, incorporating tidal effects and ring-satellite torques.The disk population was set to be monodisperse with its constituents having diameters of 1 m and equal density.
For the case with an instantaneously formed disk, assuming an initial disk mass equal to 25% of  Didymos , they found that a satellite with 0.93 Dimorphos can grow within just a few days.It then migrates outwards, away from the disk due to tidal effects and dynamical resonance with its material and slowly grows via accretion of migrating material, reaching an equal mass to that of Dimorphos after more than 40 years.When material is deposited into the disk over time, they found that long deposition timescales lead to the formation of many similarly sized moonlets that migrate outwards and merge, ultimately forming one massive satellite.This behaviour, where similar mass bodies migrate through a debris disk and merge is referred to as the growth in the pyramidal regime (Charnoz et al., 2010;Crida and Charnoz, 2012).
They revisit their previous work and expand on their model in the follow-up article Madeira and Charnoz (2024).In order to avoid issues where ejected material escapes the system before it can be accreted by a satellite (Jacobson and Scheeres, 2011a), they decreased the timescale for the formation of Dimorphos in the pyramidal regime by adjusting the duration for shed mass to spread to the Roche limit, while also altering the geometry of the resulting debris disk such that it extended to 1.5 primary radii from Didymos surface.With this alteration, they showed that their model can form a satellite with a mass of Dimorphos from a ring of only 0.04 Didymos , significantly smaller than what was indicated in the initial study.Forming in the pyramidal regime, the largest satellite undergoes a series of mergers with impacts of ∼ 2 − 3 mutual escape velocities that boost its mass over time.Using the results from Leleu et al. (2018) in a qualitative comparison, who investigated mergers between similarsized moonlets in the pyramidal regime, the authors suggest that their model can reproduce the shape of Dimorphos, as well as the contact binary satellite Dinkinesh.
The mechanism of Dimorphos forming in a lower mass, wide disk that is not narrowly confined within a region near the primary was initially suggested by Agrusa et al. (2024).Basing their model on an instantaneous, single massshedding event detailed in Zhang et al. (2017), they dynamically evolve debris disks with a few percent (2-4%) of the primary mass in a region within 1.5 primary radii from the main body surface, near the theoretical effective Roche limit for a granular aggregate with a friction angle of 35 • (Holsapple and Michel, 2006).They evolve each system using the gravitational N-body code pkdgrav (Richardson et al., 2000;Stadel, 2001).While the code is based on spherical particles it uses SSDEM to handle contacts and allows for interparticle friction, along with plastic twisting and rolling friction, approximating the behaviour of granular aggregates consisting of material with a given friction angle (Schwartz et al., 2012;Zhang et al., 2017).Providing a proof of concept simulation, detailing the evolution of a rapidly spinning aggregate accelerated by YORP with a friction angle of 40 • consisting of ∼ 90 000 particles with diameters between 4 and 16 m (consistent with observations of boulders on Dimorphos surface (Pajola et al., 2023)), they demonstrate the full process of mass-shedding near the equator of the object leading to the formation of Dimorphos-like satellite in about five days.For the remainder of the simulations, they instead let the primary be a rigid, hollow and oblate object consisting of many bound spheres with properties similar to Didymos (as a solid body).Limiting the minimum particle size to 6 m in diameter, they could perform over 100 simulations spanning 100 days as the total particle numbers varied between ∼ 4 200 and ∼ 8 400.Also varying the friction angle of the material along with the total mass of the disk, they obtain a large volume of simulations with varying parameters.
The results of the study indicate that the formation process of binary asteroid satellites is highly chaotic and unpredictable.Nevertheless, they conclude that a disk mass of 0.02 Didymos can produce moons with masses and diameters similar to that of Dimorphos, with the process becoming more consistent for a disk mass of 0.03 Didymos .Assessing the size and orbital properties of their satellites, they find that they most often end up in a state of synchronous rotation with prolate shapes, likely due to the outer boundary of their disk being near the theoretical Roche limit.Yet, a handful of the simulations produce the more atypical oblate shapes and some of the satellites end up with very prolate shapes.
In this work, we aim to further explore the details of the dynamical mechanisms that determine the shape of asteroid satellites by performing state-of-the-art N-body simulations using irregularly shaped particles and studying the mechanics that deform granular aggregates in this formation regime, focusing on analysing disk geometry, tidal encounters and impact properties for mergers.

Method
In an attempt to accurately model the complex dynamical mechanics related to the deformation of rubble pile asteroids and subsequent satellite formation, we have opted to divide the process into two stages.The first is characterised by a deformation event caused by spin-up via e.g.YORP.Numerical hydrodynamical codes, such as SPH, are a good choice for modelling deformation as they use very high resolution, and excel at tracking structural reshaping (e.g.Raducan et al., 2024a).Moreover, as the implementation evaluates contact mechanics over a smoothed continuum rather than discrete particle-particle interactions, it becomes computationally cheaper than N-body simulations.The next stage is represented by the dynamical evolution of the ejected mass settled into a debris disk.This regime is modelled using a gravitational N-body code.A number of different studies have previously combined SPH and N-body in this manner to investigate rubble pile dynamics, mainly related to asteroid collision, fragmentation and disruption (Michel et al., 2001(Michel et al., , 2002(Michel et al., , 2003(Michel et al., , 2004;;Durda et al., 2004Durda et al., , 2007;;Benavidez et al., 2012;Ferrari et al., 2022).
In order to transition between the different codes a handoff is required (e.g.Ballouz et al., 2019).Said operation is not trivial due to two main issues: (i) Since SPH particles do not represent discrete particles, but constituents of a smoothed continuum, they tend to overlap in the physical space.In turn, it is not possible to translate them one-toone as physical particles in a gravitational N-body setting; (ii) SPH simulations require a high resolution, i.e. large number of particles (on the order of ∼ 10 6 ), to properly resolve e.g.mass-shedding or deformation.However, using that many particles offers little improvement in capturing key physical effects in an N-body simulation and simply leads to a large increase in elapsed real time until it is resolved.To facilitate the hand-off between SPH and N-body, we have developed an interface that identifies clusters of SPH particles and groups them as single, angular fragments with physical properties based on their constituents.
We now go through the details of the two different codes and proceed to introduce the developed interface that performs the hand-off.Finally, we summarise the simulation setups.

Smooth particle hydrodynamics code
In this study, we use the Bern SPH code to simulate the rotation and breakup of our asteroid.It is a shockbased grid-free numerical hydrodynamics code created to simulate various geological materials.Originally developed by Benz &Asphaug (Benz andAsphaug, 1994, 1995) to model collisional fragmentation of rocky bodies, the code has been expanded with additional physics for a more accurate representation of materials.Some features of this code include various equations of state, modelling of porous (Jutzi et al., 2008) and granular materials, pressure-dependent strength (Jutzi, 2015), and a tensile fracture model making the code highly suitable to simulate and analyse asteroids with diverse material properties.The simulations produced have been tested extensively and are in close agreement with experimental results (e.g.Jutzi, 2015).A recently developed fast integration scheme allows for the modelling of lowvelocity granular flow (Jutzi et al., 2022;Raducan et al., 2024a).
For the modelling of friction and cohesion, we follow the approach by Sugiura et al. (2021) who showed that interparticle cohesion has a significant effect on the rotational failure modes of rubble-pile asteroids.To implement this effect into their SPH simulations, they opted to use a parameter referred to as the effective friction angle,  ef f .The parameter was introduced as a simplification to reduce the number of variables that the deformation modes could depend on when representing shear strength of a cohesive granular material as  = tan() + , where  is the true friction angle,  is the confining pressure and  is the cohesive strength.Assuming that the effective friction angle can account for cohesion between the granular particles in their simulations, the expression was rewritten as  = tan( ef f ).For reference, a value of  ef f = 0 • means the granular media will behave like a fluid.Sugiura et al. found that for cases where  ef f ≤ 60 • , the rotation spheroid will undergo internal deformation, while surface landslides begin to occur for values  ef f ≥ 70 • .They also note that the landslides are axisymmetric when spinning up the initial, spheroidal rubble pile over the course of a few days, meaning the body deforms into an axisymmetric top-shaped object.For longer timescales, leading to slower spinups, above one month, the landslides are localised and result in a nonaxisymmetric shape for the body.We note that this treatment for cohesion is still under some scrutiny (see Agrusa et al., 2024) and will need further study.

GRAINS -an N-body code for angular particles
To investigate post rotational failure dynamics of the asteroid system, we employ the gravitational N-body DEM (Discrete Element Method) code GRAINS (Ferrari et al., 2017(Ferrari et al., , 2020)).The software has been directly developed for the purpose of studying granular mechanics in rubble pile asteroids and allows for the use of non-spherical, angular particles.Using the c++ libraries of Chrono::Engine (Chrono, Tasora et al., 2016) in combination with a GPU acceleration scheme, GRAINS can perform numerical integrations of complex systems with thousands of particles.The use of angular particles allows us to replicate physical effects that are more difficult to capture when using spheres to represent the shape of each body, such as interlocking, off-centre collisions and particle spin orientation (Korycansky andAsphaug, 2006, 2009).Physical effects characteristic for granular media can be approximated using sphere-based codes such as pkdgrav with the inclusion of a linear spring-dashpot method with multi-directional friction and a shape parameter that gives a statistical measure of a non-spherical shape (Schwartz et al., 2012;Zhang et al., 2017).It was shown by Ferrari and Tanga (2020) using GRAINS that the use of non-spherical particles increases our ability to capture key dynamical effects in simulations of granular mechanics.A similar study for pkdgrav has been performed by Marohnic et al. (2023) who also concluded that the use of non-spherical particles, here in the form of aggregates of bounded spheres, affects the dynamics of rubble piles.As for resolving contacts between particles in GRAINS, there exist capabilities for using both non-smooth (impulse-based) and smooth (force-based) implementations.In this work we have opted to use only the Smooth-Contact Method (SMC) of Chrono, which behaves similarly to softsphere contact methods implemented in other N-body DEM codes (Sánchez and Scheeres, 2011;Schwartz et al., 2013) and is optimal for resolving low-velocity and long-lasting contacts.
Regarding the SMC method for computing contacts (i.e.force-based), it is based on a two-way normal-tangent spring-dashpot system (Ferrari et al., 2020).For this problem we use a non-linear Herzian model which is described in detail in Fleischmann et al. (2016).There are several parameters that govern surface interactions including static and sliding friction; restitution; adhesion; stiffness; and damping.We kept these parameters at their default values which have been used to verify Chrono against the behaviour of granular material in laboratory studies (Tasora et al., 2016).Energy dissipation during contact between the particles is mainly determined by the coefficient of restitution set to make collisions inelastic.Moreover, the adhesion coefficient is kept at zero, meaning that aggregates formed in our simulations are kept together purely via self-gravity.Note that the way these coefficients relate to granular behaviour in space is an ongoing field of study and will be further constrained in the future.

SPH to N-body interface
In our interface, each aggregate is created using a friends-of-friends algorithm (e.g.Huchra and Geller, 1982) with a dynamic linking length based on the combined SPH smoothing length of a particle pair, , , given by  link = 0.5(ℎ  +ℎ  ).The smoothing length is the characteristic radius of the kernel function and defines the range of influence from neighbours interacting with a given particle .For our spin-up model, this value varies between 20 m for the disk of a Didymos-like primary and 25 m for the Ryugulike primary case.To avoid outlier values in the smoothing lengths, which can occur for isolated single SPH particles and lead to unrealistic fragments, we require that each smoothing length is within 5% of the median smoothing length value in the system.While there are other prevalent clumping algorithms, such as DBScan (Ester et al., 1996) or K-means, which has been used in combination with machine learning to identify clusters in SPH simulations (Sakong et al., 2019), friends-of-friends is better suited for this problem due to the fact that it does not rely on a pre-defined number of fragments (K-means), static search distance for the fragment or struggle with large variations in density (DBscan).Hence, our implementation is completely dynamic and only relies on the smoothing length of the SPH particles without the need for additional input that could significantly affect the outcome of the hand-off.To briefly summarise the algorithm, it superimposes a cubic lattice on a three-dimensional computational box and generates a linked list by mapping the sorted values of the SPH particle coordinates along the axes of the box to the cells of the lattice.Afterwards, we choose a particle, , and check its distance to all potential neighbours inside the lattice which are within one smoothing length, ℎ  in every direction.If the distance between two particles fulfils |⃗   − ⃗   | ≤  link , they are part of the same fragment.
After it has been generated, each fragment gets assigned physical properties based on its constituents.Summing up the mass, we determine the body's barycenter and compute the corresponding position and velocity relative to the origin.We also determine the resulting orbital angular momentum of each fragment and compare it to the total corresponding sum for each SPH particle that has been used to generate it.If the total change in angular momentum in the system is greater than 10%, we reject the solution and attempt to tweak the input parameters to generate a more physically accurate output.
The ability of GRAINS to perform integrations with angular particles allows us to further improve the hand-off by capturing complex granular dynamics during the evolution of the disk.Hence, instead of simply merging the fragment constituents into a sphere, we can assign convex or concave shapes for each fragment based on the distribution of their SPH particles.This is accomplished using a geometrical algorithm based on Delaunay triangulation, known as shape (Edelsbrunner and Mücke, 1994), which has been used in previous studies of asteroid dynamics such as Ferrari and Tanga (2020).The method was also used in the SPH to N-body interface by Ballouz et al. (2019) in the scope of asteroid family formation, but just for the largest remnant in the system while we have expanded the use of -shapes to all granular aggregates in the system.

𝛼-shape implementation
To elaborate, the geometric algorithm allows us to determine the concave shapes of fragments, with their individual SPH particles acting as vertices.We require that the -shape has a density that falls inside a pre-defined range based on a target density,  target , given by (1 − ) target ≤  f rag ≤ (1 + ) target , where  ∈ (0, 1).For our purposes, we found that a tolerance of  = 0.05 worked well.If the density of a fragment is too small, the most distant vertex of the shape is removed and a new -shape is created.Note that the fragment still retains the same physical quantities, such as mass and velocity, which are calculated relative to the centreof-mass of its SPH particle members.Shapes with less than 20 constituents are instead given a randomised convex hull by filling a box with a volume based on the target density value with 30 randomly placed vertices.
Due to this approach, the high resolution of the SPH simulation is preserved while removing any overlap between particles and improving performance of the N-body code.The main drawback of -shape when implemented in interfaces such as ours has been that it needs input based on a characteristic size for individual SPH particles, which makes little sense physically as they represent infinitesimal fluid elements.Moreover, an additional input referred to as  is required, which sets the concavity of the final threedimensional shape.It can be thought of as the squared radius of a ball rolled over the vertices of a point set which yields the final surface of the shape.A small value for  increases the concavity, providing a higher resolution of the vertex distribution while the shape becomes convex when  → ∞.To avoid choosing the values of these parameters arbitrarily, we use an improved approach where  is computed individually for each fragment based on the minimum value needed for a set of vertices to form a single coherent shape.This method is available in the Alpha_shapes_3 library (Da et al., 2023) of the Computational Geometry Algorithms Library (CGAL, The CGAL Project, 2023).Furthermore, CGAL has two modes for its -shape computation: generalised and regularised.While we leave the details of these implementations to the library documentation, the regularised version will exclude any singular facets from the shape, meaning the algorithm automatically identifies and excludes outlier members of the fragment from its final geometrical form.In turn, using the regularised mode ensures that we do not accidentally combine two separate fragments as one.The formal definition for the optimal value,  opt , is that it must satisfy the following two conditions: (i) Each data point in the set is contained within or on the boundary of the regularised alpha shape; (ii) The alpha shape has a number of solid constituents that is equal to or smaller than the specified number of components in the final output (Da et al., 2023).Since this provides us with the highest possible resolution for the resulting shape, the final product is not necessarily representative of real asteroid surfaces.To end up with a smoother output, we chose to multiply  opt with a factor of 10 to serve as the initial guess for each fragment shape when performing the handoff.Artificially inflating the  value like this only marginally affects the density and volume, but removes artifical crevices and craters on the surfaces which can produce undesired interactions between fragments when resolving contacts in

GRAINS.
All-in-all, the use of CGAL in our interface greatly facilitates our modelling of granular mechanics in rubble pile asteroids.

SPH simulation setup
The SPH spin-up simulations were done for two separate cases, the first being a Ryugu-like body and the second Didymos-like.The two bodies differ in size and estimated bulk density, where Ryugu has a value of  0 = 1190 kg∕m 3 (Watanabe et al., 2019) and Didymos  0 = 2400 kg∕m 3 (Daly et al., 2023a), allowing us to investigate whether the bulk density affects the mass and shape of the secondary asteroid in the system post rotational failure.We note that the density value for Didymos was chosen before the latest values of 2760 ± 130 kg∕m 3 from Naidu et al. (2023) were published.For both scenarios, we generate spheres consisting of 50 165 particles.The sizes of the spheres correspond to 500 m and 400 m for Ryugu and Didymos, respectively.We used the fast integration scheme (Jutzi et al., 2022;Raducan et al., 2024a) which allows us to carry out the SPH simulations for sufficiently long times.The spheres were initially rotating with a period of ≈ 3.6 h.The final spinup was then performed by using an angular acceleration of ∕ ∼ 10 −9  −2 , until rotational failure occurred.
For each main body case, we used two different values for the effective friction value.Given that the goal was to investigate cases where a distinct debris disk is formed post rotational failure, we opted for values such that  ef f > 50 • , allowing us to avoid scenarios where the spin-up only results in internal deformation.
The SPH simulations were run until 3 × 10 5 s to ensure that significant mass had been lost from the primary such that it could settle in a resulting debris disk.

N-body simulation setup
Analysing several debris disks produced in the SPH simulations, we could narrow down the spread of properties to generate initial conditions for the N-body simulations.Sticking to the case of  ef f = 60 • , the average mass of the debris disk ended up being around ∼ 5% of the initial primary mass, radially extending about 1 km (2 main ) from the inner edge with a thickness of ∼ 200 m (0.4 main ).Assuming that particles within the range  main <  p <  main + 100 m would get reaccreted onto the main body, we set the inner edge of the debris disk as  main + 100 m.An example of a disk containing SPH particles compared with the angular particles post hand-off is shown in figure 1.The number density of the disk post hand-off can be found in figure 2.
At the point of the hand-off for each SPH simulation, the expansion of the disk had mostly halted from dissipation due to gravitational interactions with the primary and the rest of the shed matter, along with inelastic particle-particle contacts.As a result, using the velocities for each particle computed by clusterfinder, we found that their orbits were effectively circularised when passing the initial conditions to GRAINS.Hence, when generating disks for the standalone N-body simulations we simply gave each particle a Keplerian velocity for a circular orbit and noticed no apparent difference in the behaviour of the disk during its evolution compared to the hand-off simulations.As for the individual spin of particles, there is no way of inferring any angular velocity from the SPH data so we opted to not add any artificial spin at the onset of each N-body simulation.
For the resulting particle sizes from ten different spin-up simulations, we found that the minimum boulder diameter was ∼ 20 m, which corresponds to single, isolated SPH particles.The mean diameter varied between 25.3 m and 25.9 m for the different disks.The lower limit is inconsistent with observed sizes of rubble on asteroids, where the minimum diameters on for example Dimorphos surface appear closer to ∼ 1 m (Pajola et al., 2023;Robin et al., 2023).This is a direct result of the resolution of the SPH simulation determined by the number of particles in the primary, which sets a minimum mass and size for any isolated grains produced during the hand-off procedure.The exponential SFD we used to generate random particle distributions for the standalone N-body simulations satisfies where   = 1∕( mean −  min ).While using a distribution based on the mean and minimum diameter values from the SPH simulations could replicate the corresponding, steep SFD, this would diverge from SFDs derived using observational data.In turn, we opted for a flatter distribution with a smaller minimum value for the diameter.Noting that an average diameter of 25 m also significantly increased the number of particles up towards ∼ 10 4 , which would severely slow down the simulations, we also increased the  mean value.To keep the lapsed real time down for resolving the simulations, the final SFD was set to have  min = 5 m and  mean = 30 m.This distribution also sometimes generates particles with diameters close to 100 m, which is much larger than the maximum size of particles observed on the surface of Dimorphos.However, the maximum boulder size observed on the surface of Didymos is 93 m with a lower limit of ∼16.5 m set by the resolution of images from the DART mission (Pajola et al., 2023).The particles generated using clusterfinder can also reach these values when there has been significant clumping during the mass shedding phase with maximum values of 54.5 m for the Ryugu-like disk (see section 3.1) and 125.8 m for the Didymos-like systems (see section 4).Nevertheless, given that the SFD of particles on Didymos likely has a much lower minimum size than indicated by the low spatial resolution images of its surface along with the sizes observed on Dimorphos, it is clear that the resolution of our SPH simulations introduces a possible limitation in our N-body modelling.Consequently, we emphasize that not each body used in the simulations corresponds to an individual grain, but also clusters of particles that are represented as one monolithic body.Given that our SPH model includes some level of cohesion built into the effective friction angle implementation, clustering occurs even early on during the mass-shedding phase in accordance with Tardivel et al. (2018).The level of clustering does not change significantly if we choose an earlier or later point to perform the hand-off.The fact that these clusters survive within the fluid Roche limit could be explained if in the text, as well as the probability density of the particle SFD obtained from the hand-off f60_R2 as a histogram with 20 bins, which is our example simulation discussed in section 3.2.The SFDs have been generated using the previously defined exponential function in equation 1.The blue curve is our closest approximation to observed sizes on Dimorphos surface, while the green matches the particles from the f60_R2 handoff.The orange curve shows the much flatter distribution we use to randomly generate particles for our standalone N-body simulations.
they represent large chunks of matter, or small aggregates kept together by e.g.particle interlocking or the presence of fine-grained material ≲ 10m providing cohesive strength to the body (e.g.Sánchez and Scheeres, 2014).The overall potential presence of finer grain material in the system after the mass-shedding event would have little consequence for our simulations as they would be lost from aggregates during mergers or tidal disruptions and subsequently ejected from the system due to solar radiation pressure (Ferrari et al., 2022).In turn, we do not include any cohesion for the Nbody simulations.Moreover, our implementation does not yet include any fragmentation for these larger boulders in our system, causing them to stay intact no matter the stress they endure during the course of the simulation due to e.g.spin, tides or collisions, which might introduce unphysical behaviour.The implications of our SFD choice are further discussed in section 5.3.Three different SFDs generated with our exponential distribution, as well as a probability density of particle diameters obtained from an example hand-off discussed in section 3.2, have been plotted in figure 3. The steepest function (blue) uses a minimum and mean that yields distributions that have particles with diameters between 1 and ∼ 16 m, which are the minimum and maximum diameters observed on Dimorphos surface (Pajola et al., 2023).The green curve with  min = 21 m and  mean = 26 m approximates the SFD of particles generated from the SPH to N-body handoff.Finally, the orange curve shows the chosen SFD for randomly generating debris disks for our standalone N-body simulations.
Furthermore, the shape generation scheme involves randomly generating  vert vertices in a square box with a side of  p and then computing the resulting convex hull.This often creates particles that are coarser than observed grains on the surface of Dimorphos and other rubble pile asteroids (Robin et al., 2023).This trend and its consequences for our simulations is also discussed in section 5.3.
As for the primary body, we opted to use solid, monolithic bodies in order to reduce the total number of particles to track during the simulations1 , significantly increasing performance.For Ryugu, we used a simple sphere of radius  Ryugu = 500 m, a bulk density of 1190 kg∕m 3 (Watanabe et al., 2019) and corresponding moments of inertia for a spheroid.In the second case where we tried to recreate the Didymos-Dimorphos system, we opted to use a derived preimpact shape from the DART mission for the main body with a radius of around  Didy = 390 m.The corresponding shape also has pre-defined estimates for the moments of inertia, which we fed into the simulation while also giving it a density of 2400 kg∕m 3 (Barnouin et al., 2023).For the primary rotation rate, we opted to use the pre-impact value of 2.26 h measured by Pravec et al. (2006).
Our initial conditions post rotational failure show disks that are more radially and vertically extended than those dynamically evolved in Agrusa et al. (2024).In their work, the disks extended radially for 1.5 main and 0.05 main (20 m) in each direction along the -axis, making them narrower by half a primary radius and much flatter.This can be attributed to the nature of the rotational failure of the spinning primary asteroid.The mass-shedding in Agrusa et al. (2024) occurs right by the equator, while the primary in our model ejects mass in a wider region that extends ∼ 200 m above and below the equator.The two models differ mainly in the fact that they use N-body simulations to perform the spin-up for a Didymos-like shape while our model is based on SPH and uses spherical aggregates to approximate the shapes of Ryugu and Didymos.Moreover, the resolution of the spinup simulations in Agrusa et al. ( 2024) is higher at ∼ 90 000 particles while we use ∼ 50 000.With the many differences between the two simulation methods, it is not trivial to pinpoint what could be the major cause of this deviation in disk geometry as it would require an in-depth analysis.In the future, we aim to compare spin-up simulations between our SPH code and GRAINS to better understand any differences in the disks they generate.However, for the purposes of this article, we direct the focus to the disk evolution post rotational fission and we compare our results to Agrusa et al. (2024) in more detail in section 5.5.

Numerical setup
For the GRAINS simulations, we opted for a Barzilai-Borwein solver (Barzilai and Borwein, 1988), which is an iterative gradient descent method which offers good convergence, using a tolerance of 10 −4 and a maximum of 300 iterations to reach a solution.The time step was kept below 0.5 s, going as low as 0.01 s for complex debris disks to ensure that contacts between particles could be resolved accurately.
In order to detect the fragments in the GRAINS simulations of the debris disk evolution, we opted for a modified version of the clusterfinder algorithm from section 2.3 with a linking length of  link = max(  ,   ).In this method, we skipped the shape determination scheme and simply saved the -shapes or convex hulls representing each aggregate constituent and used them to generate a composite shape.Moreover, due to the angular nature of each grain, we determined the moment of inertia of the aggregate using the parallel axis theorem.With this value at hand, we proceeded to compute the rotational and orbital angular momenta for each detected cluster at each time step in the simulation.

Results: Ryugu scenario
We begin by introducing the results from the SPH spinup, going through the disk formation process and proceed by presenting the results from two specific hand-offs and the subsequent dynamical evolution of said disk, which ultimately leads to the formation of a satellite.After that, we show results from five less computationally demanding standalone N-body simulations.

Spin-up and rotational failure of rubble pile
All of the spin-up cases, where we let  ef f = 60 • , behave similarly as the failure of the rotating primary body occurs at around 8 × 10 4 s (in the final spin-up phase; see section 2.4).The mass shedding is localised, leading to chunks of matter getting ejected near the equator, which is consistent with previous studies of rotational failure (Hirabayashi et al., 2015;Tardivel et al., 2018;Ferrari and Tanga, 2022).However, due to the fast rotation of the body, the resulting flow of matter quickly becomes symmetric.Moreover, because of the slight increase in strength that the body gets with this effective friction angle (and correspondingly higher rotation rate at failure), the matter flowing into orbit carries significant momentum with velocities greater than the orbital speed and travels several primary radii before it settles into an orbit due to energy dissipation caused by gravitational interactions with the rest of the disk and contacts between particles.. Figure 1 shows an example of the hand-off result for a Ryugu-like case simulation with  ef f = 60 • that takes place 38.8 h (1.4 × 10 5 s) after the onset of the spin-up simulation.
The resulting disk has a width of ∼ 1 km (2 main ) measured from its inner edge at 600 m from the centre of the primary and a thickness of 200 m.The disk has a mass of 0.055 main containing 2233 particles and the fragment size distribution computed from the volume equivalent diameter fulfils  frag ∈ [20.9, 50.0] m with an average size of  mean = 25.9The distribution of disk masses, number of particles and fragment sizes for the post rotational failure debris disks created by spinning up a Ryugu-like primary using SPH.
m.The corresponding parameters for the remaining spin-up cases with the same initial properties for the primary can be found in table 1. Due to the chaotic nature of the rotational failure, each resulting disk has slightly different properties.

Dynamical evolution post SPH to N-body hand-off
We performed two hand-off simulations for two different disk mass cases, namely runs f60_R2 and f60_R6 from table 1.We use the former as an example and go through it step by step to highlight the different physical mechanisms involved.We have provided a movie of this simulation in the supplementary media repository2 , which we encourage the reader to view while going through this example.The formation process for a satellite body in this setting is largely hierarchical.As the disk settles around the primary and we perform the hand-off, the distribution of azimuthal velocities causes particles farther out from the primary to move at lower velocities.Along with gravitational collapse, this quickly leads to the formation of denser filaments in the disk over the course of a few orbits, creating a spiral pattern in the disk due to Keplerian shear.This aligns with the behaviour of pure particle circumplanetary disks generated by giant impacts (e.g.Ida et al., 1997;Kokubo et al., 2000).The converging trajectories of the particles in these dense regions lead to rapid formation of smaller gravitational aggregates with more than 10 constituents, reaching masses around 0.5% of the primary mass after only about 24 hours.In turn, the spiral patterns in the disk quickly diminish, meaning angular momentum transfer is thereafter largely dominated by particle-particle collisions (Takeda and Ida, 2001).Snapshots from the example simulation can be found in figure 4, while the growth over time for the largest remnant in the system, with mass  lr , is shown in figure 5 (cyan line).As the trajectory of growth indicates, there is significant fluctuation of the mass of the largest remnant, which is not necessarily always the same object.This is due to the accretion or loss of individual grains from the aggregate or because of another   body becoming the most massive remnant in the disk.While the initial growth is slow and largely dominated by the steady and slow process of grain accretion, the subsequent growth after 18 h is driven by mergers of aggregates and examples of such events can be seen at 18, 29, 37, 40, 57 and 73 h where there are significant increases in the mass.The two more significant mergers occur at 37 and 73 h, events that we refer to as R2a and R2c, respectively.Notably, there is also a sharp decrease in mass occurring shortly after the 40 h mark, an event we will call R2b.This is due to the largest remnant in the system getting too close to the fluid tidal disruption limit and losing much of the mass it had previously gained during event R2a and a smaller subsequent merger.Furthermore, said remnant was also on a tidally locked synchronous orbit, with respect to the primary, before this mass shedding event took place.The tidal forces stripping away a chunk of mass effectively altered the trajectory due to loss of angular momentum, putting the aggregate on a more eccentric orbit and sending it farther out into the system.The specific orbital and rotational angular momentum for the largest remnant in the system have been plotted against time in the bottom graph of figure 6.The corresponding changes in orbital elements can be seen in figure 7.In addition, the figure displays the distance from the primary's barycenter over time.The close encounter with the primary occurring during event R2b functions as a scattering event, inducing a spin for the aggregate and placing it on a seemingly stable circular orbit for the next 15 hours.At that point, around the 55 h mark, the gravitational pull from the primary once more brings the aggregate into a close encounter, leading to another merger at 57 h and a subsequent scattering event, albeit a weaker one.This Figure 6: Correlation between distance to the primary and the angular momentum of the satellite in f60_R2.In all plots, the specific events of interest referred to in the text as R2a, R2b and R2c have been marked with dotted lines.Top: the distance to the barycenter of the primary for the largest remnant in the debris disk over time.Bottom: specific orbital (blue) and rotational (orange) angular momentum over time for the largest remnant.Not that the value displayed does not necessarily always correspond to the same object, e.g. for time steps before 18 h.
effectively put it on a collision course with the secondlargest remnant in the disk.These two aggregates finally merged at the 73 h mark during the R2c event, creating the final satellite consisting of 869 grains with a total mass of  moon = 0.028 main .To evaluate the shape of the body, we used the dynamically equivalent equal-volume ellipsoid (DEEVE) measure, determining the ellipsoidal axes from the principal moments of inertia of the body.As shown by Agrusa et al. (2024), this measure provides a good estimate for the true extents of an ellipsoidal body that is an aggregate of many smaller constituents but leads to systematically larger values than the true physical extents.First, we compute the total inertia tensor for the composite object using the Huygens-Steiner theorem and then proceed to evaluate the tensor's eigenvalues to obtain the principal moments of inertia , , .The correlation between , ,  and the DEEVE axes , ,  are as follows: Given that the axes , ,  are oriented along some orthogonal coordinate system  ′ ,  ′ ,  ′ determined by the specific rotation of an object that aligns the angular momentum and angular velocity vectors, an oblate shape in the orbital plane corresponds to ∕ values close to unity.This assumption holds as long as there is little pitch and roll for the asteroid.In turn, we define the orbital plane axes such that  >  while we assume that  is pointing out of the plane.The ∕ value then determines how flat the body is, i.e. how elongated it is perpendicularly to the orbital plane.We found that ∕ = 1.27 and ∕ = 1.58, meaning the resulting satellite is less oblate and flatter than Dimorphos, which has a shape fitting the values ∕ = 1.06 ± 0.03 and ∕ = 1.47 ± 0.04 (Daly et al., 2023b).
The final shape from a top-down and edge-on perspective can be seen in figure 8.We note that the sharp features and the size of the individual grain constituents do not align with the observed particles on Dimorphos surface and may generate unrealistic details (Robin et al., 2023;Pajola et al., 2023).
Studying the simulation in detail, we found that the last merger of the largest remnant occurred in a state of rotation where the aggregate was spinning with a period of around 9 h, whereas its orbital period was 13 h.Due to the spinning nature of the target in the merger, i.e. the most massive of the Figure 7: The changes in orbital elements for the largest remnant in the binary system during its formation in simulation f60_R2.In both plots, the specific events of interest referred to in the text as R2a, R2b and R2c have been marked with dotted lines.Top: the change in semi-major axis over time for the largest remnant in the system.Middle: the corresponding change in eccentricity over time.Bottom: the distance to the barycenter of the primary.
two bodies, the final shape of the resulting body was highly dependent on the orientation of the longest axis of the target, along with the location and impact angle of the collision.In this scenario, the prolate impactor merges with the prolate target such that their major principal axes are parallel and side-to-side, creating a more oblate final shape, changing the DEEVE axis ratio from ∕ ∼ 1.8 to a value closer to 1.2.Since the minor axes in the orbital plane align and become longer in this scenario, we will continue to refer to this type of merger as short-axis, while using the name long-axis for the opposite type that occurs when the two bodies merge at the edges of their longest axis, generating more prolate shapes.
Taking another look at figure 6, there is another significant dip in specific angular momentum right before the final merger which is due to a third scattering by the primary.As can be seen, the orbit of the largest remnant stabilises after this close encounter and the fluctuation in angular momentum is dampened, ultimately placing it on an orbit with a semi-major axis of 2.52 main with an eccentricity expected to land somewhere in between 0 and 0.09.
The second simulation, f60_R6, went through a slightly different growth trajectory than the largest remnant of the example simulation, as can be seen when comparing the two tracks in figure 5.While the formation process of the largest remnant is also driven by mergers, it did not undergo any disruptions and takes more than 20 h longer to reach its mass at the final time step of 0.023 main .This can be attributed to the lower mass and fewer particles of the initial disk post hand-off.Furthermore, it ended up with a more prolate shape than our example simulation, having DEEVE axes ratios of ∕ = 1.34 and ∕ = 1.42.
Here, the shape was also largely affected by the final impact.As a smaller prolate cluster spiralled in towards the primary, it was catastrophically disrupted and an oblate remnant was scattered outwards, ending up on a trajectory that converged with that of the largest remnant.The oblate impactor then underwent a head-on short-axis merger with the spinning target close to the centre of its major principal axis with a small impact angle.This resulted in a less prolate target, with its DEEVE axis ratio going from ∕ ∼ 1.4 to 1.34.Moreover, the nature of the impact resulted in a significant slowdown in the rotation rate of the target, placing it on a tidally locked and synchronous orbit with a period of 15 h.

Standalone N-body simulations for Ryugu-like systems
To compare with the hand-off simulations, we also performed a series of pure N-body simulations where we generated a system with a debris disk orbiting an asteroid and let it propagate dynamically.The total mass was set to be equal to the estimated mass of Ryugu (Watanabe et al., 2019).Using the SPH simulation results as a basis, we distributed the mass such that 5% of the total mass was put in the debris disk and the rest in the primary.
We ran five simulations for at least 120 h and propagated them further if there were still many aggregates in range for mergers or if the orbit of the largest remnant was deemed unstable and likely to decay or become hyperbolic.Due to the random nature of the debris disk population synthesis, each disk ended up with a different number of particles from the same grain size distribution.From runs 1 through 5, the disk consists of 5458, 4926, 5530, 5465 and 5301 particles.Hence, allowing for a lower minimum particle diameter effectively more than doubled the number density of particles despite the disks having similar masses to the hand-off cases.
The mass growth for the five different cases is displayed in figure 9. From the trajectories, it is clear that within our simulation framework, an aggregate can reach a relative mass similar to that of Dimorphos on average ∼ 50 hours after the formation of the debris disk and in as little as 24 hours.This is made possible by early mergers that drastically boost the mass of each largest remnant, in turn increasing their ability to accrete both other clusters and single grains, promoting mass growth.A similar mechanism can be seen in planet formation, where collisions between migrating planetesimals in protoplanetary disks can greatly promote their ability to accrete pebbles, increasing their growth rate (e.g.Liu et al., 2015;Wimarsson et al., 2020).However, just like for the hand-off case (see section 3.2) there were many situations where an aggregate would shed particles during close encounters with the primary body in the system.In all of the five cases, there were significant disruption events that either prolonged the growth process or halted it.Looking at the magnitude of the most critical mass shedding event in each case, it is clear that the systems with the smallest relative changes in mass reached the largest final masses, i.e. runs 1 and 2 which both have satellites two times our target mass.In runs 1 and 2 there were only a few mild disruption  events where a smaller fraction (< 0.5 lr ) of the body's total mass was stripped off from tidal interactions with the primary, which slowed down the growth rate.While the largest remnant in run 1 still went through another merger, significantly increasing its mass, that of run 2 only grew via accretion of single grain or small clusters making it take an additional 60 hours until it had regained a substantial fraction of the mass lost in its disruption event, occurring at the 60 h mark.The accretion rate still being efficient for the secondaries of runs 1 and 2 can largely be attributed to them being more massive post tidal disruption.In turn, they more easily reaccreted single grains and clusters surrounding them in the disk.
In figure 10, we have plotted the relative masses of the aggregates in each system against their distance from the system barycenter at the final time step.The coloured markers represent the bodies that are massive enough and Figure 10: The distribution of masses and distance from the barycenter of the primary body at the final time step of each standalone N-body simulation.The largest remnant in each system has colour, while the grey markers represent the remaining aggregates that either are too low mass or too far out in the system (i.e.outside 5 km from the primary) to count as the most stable and massive satellite in the timescales we consider.The cyan and brown points marked with upside down and sideway triangles represents the result from the hand-off simulations described in section 3.2.
close enough to be considered the most stable and massive satellite within the scope of this article, while the grey markers show the remaining aggregates with masses greater than 0.001 main .Looking at the system configuration for these three cases, we can see an explanation as to why the growth rate slowed down.For both runs 1 and 2, the minor satellites are on much wider orbits and will likely not undergo any more mergers with the largest remnant, meaning all the smaller aggregates in the vicinity of the largest remnant have already been absorbed by the final Figure 11: The final shapes at the end of each simulation for the Ryugu-like scenario.Each ellipsoidal axis has been evaluated from the corresponding dynamically equivalent equalvolume ellipsoid to the moment of inertia of each aggregate.Note that only the coloured markers represent bodies that have a large enough mass and are situated within 5 km of the primary body.They can be compared to the final masses and positions in the different systems from figure 10.We note that these values are likely slightly larger than the those obtained from the principal axes corresponding to the actual physical extents of the bodies (Agrusa et al., 2024).
time step and growth mainly occurs through single-grain accretion.
In simulations 3, 4 and 5, the largest remnants reached relative masses above the target of 0.01 main but ended up going through catastrophic disruptions during close encounters with the primary and lost more than half their mass.The aftermaths of these tidal events differ significantly however, as the largest remnants in runs 4 and 5 retained their orbit close to the primary within the debris disk and reaccreted much of the shed mass before they each ultimately encountered another remnant formed in the tidal event and underwent a critical merger, once more putting their relative mass at a level above 0.01 main .On the contrary, the corresponding body in system 3 did not regain any of its lost particles and most of the mass ended up in orbits separated from the main debris disk.Moreover, there is again a correlation between the amount of mass held in additional aggregates still present in the system and where they are positioned with respect to the largest remnant.Looking at the final configuration for the aggregates in system 3 (see figure 10), all of the remaining clusters have similar masses between 0.003 main and 0.004 main and the largest remnant is separated from the other two, meaning it is unlikely to go through subsequent mergers and further growth as the system continues to evolve.

Satellite shapes
The distribution of the ∕ and ∕ values for the DEEVE axes of each aggregate in the pure N-body simulations and the hand-off case can be found in figure 11.The marker and colour of each data point follow the same scheme as in figure 10 and the corresponding values for Dimorphos from Daly et al. (2023b) have also been added.Only the largest remnant from run 5 with an axis ratio of ∕ = 1.08 displays a similar type of oblate shape as Dimorphos when comparing them visually.Three of the outcomes, two of them being hand-off cases, end up with similar ∕ values but are more elongated than Dimorphos and the rest have prolate bodies more similar to the majority of observed asteroid satellites.The final body of run 4 stands out as it has a value of ∕ = 1.02 and a much longer -axis yielding ∕ = 2.26.
To be able to better understand the formation process and evolution of the axis ratio ∕ over time, we have plotted them for each pure N-body case in figure 12.Note that the many changes in mass and deformation through tidal interactions, along with fluctuation between which principal axis is larger have made the data noisy.In turn, we applied a one-dimensional uniform filter to smooth it out and make it easier to interpret.Hence, some of the values might not be correctly displayed but the overall trends remain the same.One of the main features of the evolution tracks is large spikes, where the ∕ values increase drastically.Comparing the behaviour with the mass growth patterns in figure 9, these features can be compared with significant changes in mass, such as disruption events or mergers.While the large deformations in runs 1, 2, 3 and 5 are due to close encounters with the primary, leading to tidal disruptions, this is not the case for run 4. Instead of its main satellite being reshaped due to mass shedding at 160 h, the change into a highly prolate shape came from a merger with another aggregate, indicating that the largest remnant is, in fact, a bilobate satellite.We show the shape of the final aggregate in figure 13 along with the corresponding -wrap 3 for the configuration of grains.Due to a significant relative velocity between the bodies at the time of the merger (equal to 0.84 of the mutual espace velocity), the long-axis impact resulted in a subtle bilobate shape that is only evident because of narrow concavities at the centre of the body.These results agree with Leleu et al. (2018), who showed that the relative velocity and impact angle in a merger of aggregate bodies greatly affect the final shape.We discuss this further in section 5.2.
The evolution of the largest satellite shape in run 5 also provides key insight into what determines the final form.In the first 70 h of the simulation, it followed the same behaviour as the rest of the simulations where accretion of grains and small aggregates, along with tidal deformations, generated the characteristic prolate shape of asteroid satellites.However, the previously discussed catastrophic disruption that began around 75 h into the simulation (see figure 9), led to a drastic change in shape that differs in nature from the other cases.Instead of the tides deforming 3 An -wrap is a geometrical algorithm available in CGAL (The CGAL Project, 2023) which is also based on Delaunay triangulation that uses iterative refinement.It is more appropriate for capturing detailed features of a coherent body and can be imagined as the shape generated when shrinkwrapping a sheet on top of a collection of facets or points.the body such that some mass was lost from both longaxis edges, leaving a more oblate shape behind, here the disruption was so strong that only the outermost part of the body survived.The close encounter and strong disruption event resulted in the single remnant getting ejected from the system, ending up on an orbit 10 km from the primary, as can be seen in figure 10.As a result, due to the main satellite being practically disintegrated, there was a hierarchical shift in which aggregate was the most massive.As the disruption ejected a barrage of grains farther out into the system, some of the mass got accreted by the new largest aggregate, which then steadily grew and reached a relative mass of > 0.01 main after another 75 h.There are two more significant events of deformation, as it first experienced a short-axis merger right before the 150 h mark and then continued to accrete aggregates that have ∼ 10% of its mass.The final deformation that can be observed at around 180 h of figure 12 was induced by a close encounter with another aggregate which was subsequently ejected farther out into the system, which then allowed the largest remnant to relax and settle at its final shape, being an oblate but a less flat version of Dimorphos with axis ratios of ∕ = 1.08 and ∕ = 1.32 and a slightly higher relative mass of 0.014 main .

Formation patterns for satellites in a radially extended disk
Looking at the formation histories and final shapes of the largest satellites, it is clear that the process is highly chaotic, which agrees with the findings of Agrusa et al. (2024) and can lead to different shapes and number of satellites just from a small change in the initial conditions of the simulation.In order to make more sense of the formation mechanisms involved and find conditions that favour the formation of oblate satellites, we formulate three general patterns for our type of debris disk that will aid us in the analysis: • Insignificant disruption history (IDH): The main satellite forms farther out in the system, which results in little-to-no disruption caused by tidal interactions with the primary and its growth is steady throughout the simulation, where it accretes single grains and merges with other, smaller aggregates.Given its distance from the primary, the body is less likely to become tidally locked, meaning that its angular rotation rate is mainly decided by impacts and mergers that either stop or induce rotation, depending on the impact angle and whether they are short-or long-axis.
The mass ratio, and in turn diameter ratio, between target and impactor,  impactor ∕ target is likely small.An example of this formation pathway is f60_R6 in figure 5.
• Mild disruption history (MDH): Being closer to the primary, the largest remnant is more susceptible to tidal disruption and undergoes one or a few such events where the tidal forces deform the aggregate, elongating the body until it structurally fails and loses Table 2 The distribution of disk masses, number of particles and fragment sizes for the post rotational failure debris disks created by spinning up a Didymos-like primary using SPH.
less than half its mass.Post disruption, the body regains some of its lost mass through single-grain accretion and mergers leading to a gradual change in shape.
Given the mass shedding of the largest remnant, the ratio  impactor ∕ target during mergers is likely greater than for the insignificant disruption history scenario, leading to more significant fractional changes in  lr ,  lr and ∕ values.The rotation rate of the body is kept low due to tidal locking but impacts and mergers can induce higher rates.Hand-off run f60_R2 along with N-body runs 1 and 2 are examples of this type of formation history.
• Catastrophic disruption history (CDH): Due to close encounters with the primary, the largest remnant undergoes strong tidal interactions and loses more than half its mass, which significantly slows down the growth process.The remaining body tends to be less prolate, but due to mergers having a  impactor ∕ target close to unity, the fractional change in  lr ,  lr and ∕ can be very large, causing significant deformation.In some cases, the disruption leads to a complete halt in growth after ejecting a significant portion of the disk mass from the system.N-body runs 3, 4 and 5 display these types of behaviour.
In general, short-axis mergers favour oblate shapes while long-axis mergers often lead to prolate shapes and the creation of contact binaries, which is further explored in section 5.2.We will henceforth refer to the different formation history archetypes as IDH, MDH and CDH, respectively.

Results: Didymos scenario
Here, we perform the same type of analysis as we did in the previous section but for a Didymos-like primary, using the newly established formation patterns to describe the evolution of each system.The main difference from the general Ryugu-like setup is that each grain, along with the primary, has a higher density, which leads to on average smaller particles in the system.
For the Didymos-like primary, we performed five SPH simulations, again with  ef f = 60 • and the properties of the resulting disks can be found in table 2. There is a clear difference in the number of particles and total mass shed from the primary post rotational failure which can be seen in the average disk mass being 0.068 main compared to the 0.045 main for the Ryugu-like primary.While the topic of rotational failure still needs dedicated research and any analysis related to this mechanism and our SPH simulations in particular are far from conclusive, there are still patterns here worth mentioning.Given that the density of the primary is higher and the radius is smaller, the structural integrity of the body is stronger than for the Ryugu-like primary.This aligns with the density dependency of the critical spin rate for rubble pile asteroids which can be inferred from lightcurve amplitude using  crit ∕3.3 h = √ (1 + Δ)∕, where Δ ∼ 2.5 log(∕) is the lightcurve amplitude from observed magnitude converted to the axis ratio ∕ (Walsh, 2018).While the debris disks for the Ryugu-like scenario were principally symmetrical, this is not the case for all debris disks in the Didymos-like case.Only two out of five are symmetrical, namely f60_D0 and f60_D4 which also have smaller maximum particle diameters.We show an example of a nonaxisymmetric disk in figure 14.From the righthand plot which displays the output from the clusterfinder algorithm, it is evident that the disk is denser towards the bottom of the plot where it contains two large boulders with large diameters of 83 and 93 m.This shows us that the nature of the mass shedding for this body was more localised, with large chunks of particles being ejected from its equator (as postulated by Tardivel et al., 2018), whereas the symmetric disks have a history of global symmetric landslides.Moreover, the distinct extended features on the left side of the debris disk also confirm that the rotational failure of the primary was initially more localised in this scenario, starting at a specific location, allowing the corresponding shed mass to spread farther out in the system than the rest of the debris.

Dynamical evolution post hand-off for Didymos-like systems
We evolved two of the hand-off systems dynamically, choosing one symmetric and one nonaxisymmetric disk in simulations f60_D0 and f60_D2, respectively.Other than their shapes, the geometry of the disks from each simulation is similar to that of those formed from the rotational failure of the Ryugu-like bodies.They radially extend for about 1 km (2.5 main ) from the surface of the primary and have a thickness of ∼ 200 m, with a few individual particles being farther out for the nonaxisymmetric case, as can be seen in figure 14.The resulting growth tracks for the two systems can be found in figure 16.Just like in the case of the Ryugulike system, there is one case of IDH (f60_D0) and one of MDH (f60_D2).While there was a mild disruption for the case of f60_D0 around 30 h into the simulation, it had little consequence for the growth rate and final shape of the body.Comparatively, the formation process of the largest remnant in f60_D2 got severely slowed down by its disruption events.Moreover, the body ended up on an orbit less than a primary radius from the main body after the first instance of mass shedding.This explains why it mainly grew via accretion of individual grains and small clusters from 60 hours onward until it underwent a second disruption event at 125 h, leading to fragmentation and migration outwards for the largest remnant.Looking at the change in the ∕ axis ratio over time for the two cases in figure 17, we observe how large of an effect the final disruption event had on the resulting shape of the body.The tidal forces split the moon into two smaller aggregates along its longest axis that eventually recombine in a short-axis merger, making the final body oblate with axis ratios of ∕ = 1.14 and ∕ = 1.48.We note that the largest satellite might still undergo further mergers, as two additional aggregates remain in the system on similar orbits.Yet, as the system is unresolved after 148 h and it illustrates the importance of mergers for obtaining oblate shapes, we opted not to continue the dynamical evolution while recognising that the final shape might be different for this satellite when the debris disk has been cleared of material.
The growth patterns for each largest remnant in the two systems are distinctly different, which could potentially be tied to disk geometry, given that the two disks had similar masses and number of particles.For all the symmetric disk scenarios, the growth trajectories are largely similar up until the 50 h mark of the simulations, the sole exception being f60_R6 which has less mass and fewer particles than the rest which significantly slows down the growth process.However, for the nonaxisymmetric disk in f60_D2 growth slowed down after the 20 h mark.This could also be tied to the amount of mass lost from the system early on, i.e. single grains and aggregates that end up with hyperbolic orbits and have semi-major axes larger than 10 km.By the 20 h mark, system f60_D2 had ejected 0.02 disk which was two times as much mass as ejected by f60_D0.By the 40 h mark, this value had increased to 0.04 disk which was three times as much.The values intersect once more towards the end of the simulations, but the effect of the early ejection of mass can be significant.Despite the amount seeming small compared to the largest remnant masses at the time, it still corresponds to several small aggregates that could have increased the accretion efficiency of the largest remnants in the system by boosting their mass and diameter.Nevertheless, this is a complex problem where the specific geometry of the overdensities in the disk may also have a large influence on the growth rate.In turn, the effect of disk asymmetries will be a topic for future research as our current understanding is inconclusive and our SPH model is still under development.

Pure N-body simulations for Didymos-like systems
For the five pure N-body simulations, we once more opted to use the same disk mass of 0.05 main as in the Ryugu-like case.The SFD for the generated grains was also kept the same for consistency.The higher density along with the lower total disk mass led to a significant decrease in number of particles, averaging 2454 across the five simulations.Examining the growth tracks for each of the simulations, which are presented in figure 18, three of the formation histories are IDH (runs 1, 2 and 3), while run 4 is an MDH The algorithm has generated -shapes and convex hulls with volumes matching the densities and masses of SPH particle clusters.Zooming in will show the resulting 3D shapes with black grid lines marking the edges that separate the faces in blue.
case and run 5 is a CDH where the satellite gets completely disintegrated.The final masses and positions relative to the barycenter have been plotted in figure 19.Lastly, the resulting shapes for the aggregates at the end of the simulations are given in figure 20 while the evolution of the elongation in the orbital plane, i.e. the ∕ axis ratio, for the largest remnants can be found in figure 21.Compared to the Ryugulike system cases, the final bodies generally have higher relative mass, with three out of seven simulations generating secondaries with masses greater than 0.03 main .Notably, the aggregate most similar to Dimorphos from run 4 is an MDH scenario, growing much more slowly and steadily after its disruption event, primarily through single-grain accretion with one last merger setting its final mass at 0.017 main .
With principal axis ratios of ∕ = 1.10 and ∕ = 1.38, it exhibits the same type of flattened oblate shape that we are looking for.Studying the evolution of the ∕ values for the secondaries grown in Ryugu-like and Didymoslike cases in figures 12 and 21, a key difference between them is that the aggregates are initially less prolate for the latter.Overall, strong tidal interactions with the primary are less common due to its smaller fluid Roche limit since  Roche ≈ 2.44 main ( main ∕ p ) 1∕3 (Roche, 1847), where  p is the density of a particle in orbit, which explains why we have more IDH scenarios with these initial conditions for the geometry of the disk.We note that this expression does not take into account the oblateness of the shape model we use for Didymos and address this further in 5.1.Looking at the growth tracks for the bodies once more, comparing them to the corresponding evolution of ∕, it becomes clear that deformations in this case are primarily driven by mergers rather than tidal forces.With this in mind, we can conclude that both runs 1 and 2 have produced bilobate satellites, albeit significantly less prolate ones than the case of the largest moon in Ryugu-like run 4.We study these merger cases in more detail in section 5.2.Properties such as final mass, orbital elements and geometrical values for each largest remnant in our simulations, along with the corresponding Dimorphos values can be found in table 3. We signify each type of shape for the aggregates as oblate spheroid (OS), prolate spheroid (PS) and bilobate (BL).

Discussion
Based on our results and analysis in sections 3 and 4, we can identify initial disk conditions and dynamical evolution histories that facilitate the formation of oblate spheroids such as Dimorphos.Working backwards from how we obtain our aggregates with ∕ < 1.15, they all went through at least one critical merger that led to significant deformation.As previously mentioned, a formation path where a body is formed via slow accretion of material periodically shed from the primary body in the system logically favours a prolate shape.This also seems to apply to cases where the body rapidly forms via mergers in a debris disk after a gravitational instability has occurred and then mainly grows through the accretion of smaller aggregates and single grains during the later stages of its formation.This comes from the fact that if the body does not undergo any mergers with larger aggregates and remains close to the fluid Roche limit, it will Table 3 The properties of the largest moons at the final time step of each simulation along with the corresponding values for Dimorphos, which are highlighted in bold.The shape types correspond to oblate spheroid (OS), prolate spheroid (PS) and bilobate (BL).As for the formation history archetype definitions, see section 3.3.2.
in most cases end up on a tidally locked, synchronous orbit.This aligns with the results of Agrusa et al. (2024), who investigated the formation of binary asteroid systems using a large number of simulations starting from a more narrow and thinner debris disk with different masses around a Didymoslike primary.Much like in our model, in the cases where they ended up with an oblate Dimorphos-like secondary, the later stages of formation was primarily driven by tidal disruption events and/or mergers with favourable geometry.The same results were found by Madeira and Charnoz (2024), who expanded on their 1D ring-satellite model study Madeira et al. (2023) to consider cases where mass is deposited into the debris disk over different timescales.When the mass shedding is instantaneous, the satellite that forms does not undergo any critical mergers that significantly could alter its shape, rendering it prolate.
Looking at the formation histories for the largest remnants in our simulations and their corresponding shapes, it is clear that a body growing slowly via IDH with few to no critical mergers during the later stages of its formation The change in principal axis ratios ∕ for the corresponding DEEVE shape of each largest remnant over time for the hand-off simulations for the Didymos-like system.Due to excessive noise, the data has been smoothed using a uniform one-dimensional filter to facilitate interpretation.Note that this may skew the individual data points from their actual value.The evolution prior to 25 h has been omitted due to the naturally chaotic behaviour of the formation process.Comparing the tracks to those representing the largest aggregate mass over time in figure 16, there are clear overlaps between large mass shedding or accretion events and spikes in ∕ values.
would favour a prolate shape, such as e.g.Ryugu run 2. Most of our simulations go through shape-defining deformation events such as mergers and tidal disruptions, which create the many discontinuities observed in figures 5, 9, 16 and 18 which display the growth tracks for each largest remnant.Again, focusing on e.g.Ryugu run 2 after the 75 h mark, the shape of the bodies that reach a high mass early on will Figure 19: The distribution of final masses for all major aggregates in the simulations for the Didymos-like case plotted against their distance from the barycenter.The largest remnant in each system has colour, while the grey markers represent the remaining aggregates that either are too low mass or too far out in the system (i.e.outside 5 km from the primary) to count as the most stable and massive satellite in the timescales we consider.The cyan and brown points marked with upside down and sideway triangles represents the result from the hand-off simulations.
be more weakly affected by each subsequent merger as each aggregate it collides with will most often have a lower mass.Hence, while the formation method is highly chaotic, a debris disk that favours mergers between bodies of similar mass should have a higher likelihood of forming an oblate spheroid.To further this claim, even though similar mass mergers far from guarantee an oblate shape, as the outcome is highly dependent on impact velocity and angle (Leleu et al., 2018), tidal locking in our simulations is usually broken by accretion events that induce a higher spin-rate for the largest aggregate and push it out into a wider orbit.In Note that only the coloured markers represent bodies that have a large enough mass and are situated within 5 km of the primary body.They can be compared to the final masses and positions in the different systems from figure 19.We note that these values are likely slightly larger than what we would obtain from the principal axes corresponding to the actual physical extents of the bodies (Agrusa et al., 2024).
turn, the subsequent accretion of single particles and small aggregates for the spinning satellite will be isotropic and often further reduce its ∕-value, making the body more oblate as time passes and the disk is slowly depleted of material.
An example that reinforces this idea is Didymos run 3.In the same manner as Ryugu run 2, it does not undergo any significant mergers after the 75 h mark.However, when studying the trend of the ∕ value for the largest remnant in this system, it tends towards a more oblate shape as it accretes single grains and smaller aggregates.This is due to its final major merger that broke its synchronous rotation, causing a yaw rotation back and forth.As a result, the single grain and small aggregate accretion no longer had a preferred direction along its longest axis, making tides less efficient at retaining the prolate shape.Over the later stages of its evolution, the aggregate went from having an ∕ value of ∼ 1.4 to 1.29.Nevertheless, at the final time step, almost the entire disk was depleted of material showing us how the mechanism where single grain accretion makes a satellite more oblate over time depends on the remaining disk mass and the initial ∕ value of the body.
Based on these observations, we will proceed with the discussion to eventually motivate which properties a debris disk should have to increase the probability of forming an oblate spheroid satellite.After that, we go through which developments could improve our model and, as a result, our understanding of this complex problem.

Debris disk properties and tidal influence
Judging from our different simulation setups, including hand-off cases and N-body simulations for Ryugu and Didymos-like bodies, the initial geometry of the disk seems to highly affect the formation history and the final shape of the main satellite, more so than the debris disk mass, the number of particles and their size distribution.This is also consistent with the outcome for Agrusa et al. (2024), who found that the final secondary shape and formation process for disks with masses of 0.02 main , 0.03 main and 0.04 main are highly similar, with the lowest mass disk producing more prolate bodies than the other cases.The main difference for these cases is the slower growth rate for lower disk masses, while a higher disk mass leads to greater final relative masses, semi-major axes and eccentricities of the formed secondaries.A more radially extended disk, with its edge situated beyond the fluid Roche limit for the primary body, allows for a higher survival rate of aggregates which leads to a higher number of bodies in the system, more mergers and more possibilities for deformation.The type of chaotic rotational failure that we have explored also results in higher growth rates than for previous studies such as Walsh et al. (2008), where the mass shedding was slow and continuous rather than instantaneous as a result of their implemented spin-up model.
We have also observed that the initial shape (at ∼ 24 h into the simulation) of the largest remnant before any tidal events or mergers is less prolate for a Didymos-like disk, which according to our SPH simulations extends farther beyond the fluid Roche limit than the disk created in the Ryugu-like rotational failure event.The edges of the disks in these cases are located at ∼ 3 main from the primary center-of-mass in the Ryugu case and ∼ 3.5 main in the Didymos case.Moreover, the moons formed in our model have a mean semi-major axis of ∼ 3.03 main by the final time step, which aligns with observations of binary asteroid systems (Monteiro et al., 2023) and migrate moderately over the course of the simulations unless they experience a strong collision or tidal event, as can be seen from an example in figure 7.In turn, since our aggregates primarily grow through mergers with other aggregates, angular momentum and energy transfer mainly occur via soft particle-particle contacts and lead to quick migration outwards beyond the fluid Roche limit for aggregates formed inside of it, which for our cases would be located at ∼ 2.44 main .Since the debris disk ends up more radially extended for the Didymos case relative to the primary radius, the proto-moon will most often form close to the Roche limit or outside of it.As a result, any deformation due to tides is less likely to occur throughout the growth process, which again explains why a majority of the formation histories for the Didymoslike systems are IDH, while most Ryugu-like systems are MDH or CDH.This provides a different pathway for forming Dimorphos-like bodies since single-grain accretion and mergers with lower mass moonlets for an aggregate that is initially less prolate and faster rotating than a tidally locked body will be isotropic and make the shape more oblate over time.We discuss the effect of the initial shape of the proto-moon for deformation via mergers in the next section.Despite this trend for the Didymos-like systems, an overall pattern for our simulations is that only formation histories that are MDH or CDH can produce oblate spheroids similar to Dimorphos in mass and form, emphasising the importance of close encounters with the primary to deform the body from its initially prolate state before it enters another stage of accretion which leads to the final shape.This also indicates that Dimorphos could form over longer timescales than explored in this work after the debris disk has been more or less cleared of matter.In this scenario, a prolate body undergoes orbital decay through tidal interactions with the primary and then has a catastrophic close encounter which results in mass shedding and scattering, leaving behind a less massive but oblate shape on a wider orbit, a phenomenon we have observed in several of our simulations before the disk is depleted.While the theoretical fluid Roche limit provides a convenient qualitative measure for where the tidal interactions will occur throughout our simulations, using a static limit for the tidal disruption for Rubble pile asteroids fails to capture the intricacies of aggregate dynamics (Holsapple andMichel, 2006, 2008) and has been explicitly numerically modelled and investigated in several (Richardson et al., 1998;Walsh and Richardson, 2006;Movshovitz et al., 2012;Schunová et al., 2014;DeMartini et al., 2019;Zhang and Michel, 2020), primarily regarding tidal disruption of rubble piles during close encounters with planets, but also with white dwarfs (e.g.Li et al., 2021).The effect of using irregularly shaped grains on the integrity of the rubble pile during a tidal encounter was the topic for Movshovitz et al. (2012), who in the context of rubble pile asteroids emphasised that particle interlocking, dilatancy 4 and more complex packing of granular media contributes to the structural strength of the aggregate, making it more difficult to disrupt than a collection of spherical particles.With this in mind, it is unsurprising that the aggregates in our simulations often undergo close encounters with the primary, having orbital pericenters closer than  Roche without being disrupted, as can be seen when comparing figures 5 and 7 before the critical merger at 72 h.Furthermore, we also expect that this contributes to the total number of aggregates available for collisions in our debris disks, as they will not be disrupted as easily during each stage of the formation process and that this number is greater than for a disk consisting of spherical particles, which we aim to confirm in a future study.As we have already established, the more aggregates that exist in the system, the likelier it is for the largest remnant to have a geometrically favourable merger that gives it a more oblate shape.

Asymmetries for the primary body and disk
The problem of identifying the tidal disruption limit becomes more complicated when introducing asymmetries for the primary shape.The model we use for Didymos is the official pre-derived impact mesh for the DART mission which is an oblate, nonaxisymmetric spheroid (Barnouin et al., 2023).Not only does such a shape affect the position of the theoretical fluid Roche limit, but its asymmetric inertia tensor can also lead to extra angular momentum transfer to the disk and affect its evolution (Takeda and Ida, 2001).The asymmetries on the surface of an asteroid induce chaotic diffusion of the eccentricity of particles in low orbits (Sicardy, 2020;Rollin et al., 2021;Madeira et al., 2022;Ribeiro et al., 2023), which is expected to affect the geometry and velocity of impacts between particles5 .In order to investigate whether or not the effect is significant to the formation process, we intend to expand upon our hand-off scheme to also provide an shape with a corresponding inertia tensor for the resulting contour of the primary post rotational failure and compare the evolution of two disks with identical initial conditions but different primary shapes, one that is spherical and one that is based on the aforementioned nonaxisymmetric post fission body.This is a non-trivial problem, as the primary often has formed temporary clumps of SPH particles along its equator at the time step for the hand-off which would not be representative of the final shape when the primary aggregate has settled after its initial deformation.
As previously discussed in section 4.1, the evolution of a pure particle debris disk seems to be affected by initial nonaxisymmetries, such as the one in figure 14.While we lack a statistically significant volume of simulations to confirm this phenomenon, the unique growth track pattern for the corresponding disk (see figure 16) provides us with a strong indication that asymmetries lead to significant changes in the satellite formation patterns and warrants further investigation, especially given the lack of literature dedicated to this problem for pure particle disks and circumasteroidal disks in general.A critical difference with the asymmetric disks is that they already have large fragments of SPH particles at the hand-off point, which explains why the largest grain diameter is at least two times greater than for the symmetric disks (see table 2).As previously mentioned, while the initial sizes of these grains diverge compared to observed boulder sizes on Dimorphos (Pajola et al., 2023;Robin et al., 2023), they represent aggregates that already have formed during the mass shedding event.The presence of these large grains may have a significant effect on the shape and mass of the final moon as they are monolithic and will not get disrupted by any tidal effects or impact events.In turn, for future handoff implementations, converting these larger grains into aggregates consisting of smaller constituents with smoother shapes and sizes could prove to be important.

Particle number dependency
As shown by Kokubo et al. (2000) and Takeda and Ida (2001), the number of particles has little effect on the evolution of circumplanetary particle disks confined within the Roche limit for values 10 3 ≤  ≤ 10 4 .However, Sasaki and Hosono (2018) found that disk behaviours and moon formation histories can differ for very high resolution Nbody simulations with  > 10 5 .Nevertheless, we cannot argue with certainty that this would be the case for our model, as these scenarios used spherical particles to investigate a problem on a different scale than ours.Moreover, even when using an SFD more true to observations with e.g. min = 3 m and  mean = 15 m, this would not yield a number greater than ∼ 2.5 × 10 4 .Nevertheless, in the future we aim to perform simulations with higher particle numbers, restricted only by computational limitations, to verify whether or not this directly affects the disk evolution and formation processes for the moons in our model.

Effect of merger geometry for aggregates
Studying the ∕ value evolution for our simulations in figures 12, 17 and 21, it is evident that we have several types of mergers that affect the shape of the largest remnant in different ways.More importantly, we can distinguish between mergers that lead to a more oblate body and those that lead to a more prolate body.To better understand which merger properties govern the final shape of the body, we investigate the ∕ value of the largest remnant before and after critical mergers identified by comparing the ∕ value evolution over time with the mass growth tracks in figures 5, 9, 16 and 18.For consistency, we fix which of the axes are  and  prior to the impact, meaning that we do not alter this definition even if there is a change such that  <  post merger.
In figure 22, we have plotted several impact properties of a few select critical mergers from our simulations.The selection criterion was that the merger had to be the final significant reshaping event for the largest remnant during its formation history.In each plot, the grey marker represents the ∕ value before the merger, while the coloured marker is the same ratio afterwards by the final time step of the simulation when the resulting aggregates have settled.In the top left plot, we show the impact velocity, i.e. | ⃗  target − ⃗  impactor |, normalised by the mutual escape velocity given by Here, the target body is always referring to the largest of the two.Notably, the impact velocities we observe in The DEEVE axis ratio ∕ value before (grey) and after (coloured) the critical merger which determines the final shape of the considered secondaries at the final time step of our simulations, compared with four different impact properties.Arrows pointing upwards generally mean that the body became more prolate (elongated) after the merger, while an arrows pointing downwards indicates that a deformation to a more oblate shape.Note that we fix the definition of each axis as  or  before the merger, which is why the ∕ axis ratio can be below unity for some cases (e.g.Didymos run 1).Upper left: the ∕ axis ratios and the corresponding impact velocity divided by the mutual escape velocity.Upper right: same axis ratio change plotted against the impact angle between the two velocity vectors.Lower left: the axis ratio compared to the mass ratio of the impactor and target.Lower right: the distance of the merger from the barycenter of the primary in terms of main radii.
our simulations are much lower than for previous studies exploring deformation through mergers such as Leleu et al. (2018), who investigated collisions between initially spheroidal moons in the rings of Saturn for impact velocities between  impact ∕ esc values between 1 and 4, leading to growth in the pyramidal regime.We attribute this to the fact that our work explores a different regime entirely, where bodies on similar orbits near the primary undergo comparably low-energy gravitational encounters, leading to soft mergers rather than the violent, higher-velocity impacts that occur for the migrating moonlets in the rings of Saturn and cause significant deformation.This is quantified in the impact velocity being approximately proportional to  −1∕3 , where  = ( impactor +  target )∕ main (Leleu et al., 2018).
Furthermore, the dissipative nature of the dynamics in our dense disks with numerous inelastic particle-particle collisions will also contribute to the low velocities observed as it effectively dampens any orbits with larger eccentricities and inclinations, keeping the system from getting more dynamically excited.Leleu et al. found that collisions between gravitational aggregates at these velocities lead to the characteristic ridges and deformations that can be seen from images of e.g.Pan and Atlas.Hence, if Dimorphos formed via a series of higher-velocity mergers, we would expect to observe distinct ridges and lineaments on its surface, yet it seems largely homogeneous (Barnouin et al., 2023) with no spatial preference for boulder density and sizes (Pajola et al., 2023), indicating a different formation pathway was involved 6 .Moreover, for lower impact velocities less than 0.75 esc , Leleu et al. get unstable bilobate objects that end up separating because of tidal or Coriolis forces due to the close proximity to Saturn.While we cannot speak for the long-term stability of our objects, it is clear that our lowest velocity mergers below this value lead to more oblate shapes rather than generation of more prolate bodies.Intuitively, a low-velocity merger between two spherical bodies at a low impact angle should produce bilobate aggregates, which also is seen in other studies such as Jutzi and Asphaug (2015).
The dissimilarity in outcomes must accordingly be related to the prolate pre-impact shapes with ∕ > 1.2, as well as the orientation of the longest axes of the two merging bodies.Indeed, the differences between the pyramidal and our hierarchical regime are emphasised by the initial shape of the aggregates in our mergers, as well as their masses which can be seen in the lower left plot.The largest mass ratio  impactor ∕ target is 0.84, while the mean value is 0.51 and we observe that lower mass ratios lead to smaller effective changes in ∕ despite having a significant spread in  impact ∕ esc .In the upper right plot, the change in axis ratio has instead been plotted against the absolute angle between the velocity vectors of the two objects, denoted .Seeing this angle is smaller than 23 • for all our cases, we can deduce that the orbits of the target and impactor are widely similar prior to the merger, which also aligns with the low impact velocities we find.Furthermore, there is no trend as to which impact angles produce which shapes or as to how far from the primary centre-of-mass the merger occurs,  merger , which has been plotted in the bottom right plot.Note, however, that the far-in locations of the mergers for runs f60_R2 and Ryugu run 4 do not coincide with their orbits by the final time step as each of these mergers is soon followed by a migration outwards farther away from the primary.
The reason why  cannot directly predict the post-merger shape of the body is due to the fact that it tells us nothing of how the asymmetric target is rotated relative its trajectory and if it collides with the impactor at its short-axis or long-axis.To verify this idea, we define the angle  p between the longest principal axis pre-merger, , and the relative velocity vector.An angle of  p = 0 or 180 • then corresponds to a relative velocity vector completely aligned with , while a value of 90 • tells us that the relative velocity vector is parallel to the shortest axis, .We plot this value against the change in ∕ in figure 23.There is a correlation between a reduction/increase in ∕ value of the largest remnant and the  p value of the merger, indicating that short-axis/long-axis mergers indeed favour the creation of more oblate/prolate bodies.We study the two outlier cases, Didymos runs 2 and 4 which are both long-axis mergers per our definition but still show a reduction in ∕.The former does not become an oblate spheroid, but rather a prolate, bilobate body where the pre-merger axis  instead is the 6 We note that any ridges or lineaments on the surface of Dimorphos could have been homogenised since the coalescence of the satellite due to exchange of material from Didymos to Dimorphos (Yu et al., 2019;Madeira et al., 2024) or isotropic single-grain accretion during the later stages of its formation history.longest axis.As for Didymos run 4, the impactor is a highly elongated but low-mass aggregate on an eccentric orbit that undergoes a close encounter with the primary, enters the Hill radius of the target and is subsequently slowed down before it is tidally disrupted and accreted, which makes it behave unlike any of the typical short-or long-axis mergers we observe.
While the orientation and shape of the impactor will also affect the final structure of the largest moon in our systems, the generally low mass ratio between the two objects renders the orientation and initial shape of the target far more important.This becomes evident when comparing figure 23 with figure 24, where we instead have plotted the angle between the relative velocity and the largest principal axis of the impactor.Only two of the impactors have distinct short-axis mergers with  p,impactor close to 90 • while five of the corresponding mergers generate more oblate bodies.Instead, their orbits have the largest effect on the outcome of the merger.Looking at the cases which generate more prolate bodies and contact binaries such as Ryugu run 4, as well as Didymos runs 1 and 2, this becomes evident when combining their distribution of impact angles with the corresponding  p values.Here we see that the longaxis mergers Ryugu run 4 and Didymos run 1 have the largest impact angles as well, showing that the impactor has a dissimilar trajectory to the target.Indeed, investigating the merger in detail confirms that the impactor in Ryugu run 4 has a wider orbit than the target and approaches from the outside, while the opposite is true for Didymos run 1 where the impactor instead has an orbit closer to the primary and intersects the orbit of the target from inside.
From these results, there is an indication that a formation model more likely to generate low-velocity, short-axis mergers is efficient at forming oblate bodies.Tying this to the geometry of the disk, a wider initial radial extension similar to those generated by giant impacts onto planets (e.g.Ida et al., 1997;Kokubo et al., 2000) facilitate the fast formation of multiple aggregates.Due to the hierarchical nature of the process, one aggregate will quickly grow more massive than the other satellites in the system and often becomes prolate and tidally locked on a synchronous orbit.Undergoing one or more mergers with other satellites (depending on the mass ratios) under favourable circumstances such as low impact velocity, similar trajectories and  p values close to 90 • will then reduce the elongation of the body significantly and can in some cases generate oblate spheroids.

Bilobate satellites
Further recalling our analysis from the section 5.1, the IDH, MDH and CDH scenario outcomes are dependent on the position of the Roche limit relative to the outer edge of the disk and are as a result system-specific.Because the strength of the tidal forces largely affects the initial shape prior to mergers, forming most aggregates farther away from this limit will generate less prolate bodies that are less likely to have synchronous orbits, which is the case for our Didymos-like systems.For the corresponding IDH cases where growth is dominated by isotropic accretion of single grains and low-mass aggregates, the moonlet will become more oblate over longer timescales.However, for systems where there still are additional massive satellites, it follows that the long-axis mergers will be likelier, which in turn increases the probability of forming bilobate bodies.Hence, for our rotational failure model, primaries with smaller radii and higher densities with Roche limits closer to the surface should be more likely to have contact binary moons.This is an important topic for further research, especially given the recent discovery of the bilobate satellite Selam of the (152830) Dinkinesh system (Levison et al., 2024) by the Lucy mission during a flyby (Levison et al., 2021).Dinkinesh has a distinct top-shape indicative of rotational failure and an effective diameter of 720 ± 24 m making it smaller than both Didymos and Ryugu.The satellite Selam is orbiting the primary at a distance of ∼ 9 main and has a distinct bilobate structure where the components have diameters of 212 ± 21 m and 234 ± 23 m (Levison et al., 2024).Despite the much greater separation between primary and secondary, the Lucy images show that the lobes of Selam align radially with the primary, indicative of a tidally locked orbit.In turn, there is a question of whether or not Selam formed closer to the primary and then migrated outwards via dynamical effects such as binary YORP (BYORP, Ćuk and Burns, 2005) or tidal interaction (Jacobson and Scheeres, 2011b) or if it formed farther out in the pyramidal regime according to the model of Madeira and Charnoz (2024).Moreover, the analysis from Levison et al. (2024) also shows that the mass ratio between the satellite and primary, assuming equal density for both objects of  = 2400 ± 350 kg∕m 3 , is 0.06 which is significantly higher than any of the satellites we can form during the course of our simulations (see table 3) as we assume a lower initial disk mass.Nevertheless, judging from our spin-up simulations of Didymos in table 2 and the work of similar studies (e.g.Sugiura et al., 2021;Hyodo and Sugiura, 2022), it is entirely possible for primaries with higher bulk densities to shed such large masses via rotational failure when introducing some weak cohesion to the primary.Hence, by performing followup SPH simulations with a primary based on Dinkinesh with higher resolution, we can get a better understanding of whether or not forming such a massive contact binary satellite is achievable with the circumasteroidal disk model.Given the low velocity of our impacts and the fact that our model can form weakly bilobate bodies even with lowresolution simulations, we believe this to be the case.
With the pronounced equatorial ridge of Dinkinesh in mind, a history of deformation and mass shedding due to rotational spin-up is very likely.In turn, we expect that the existence of Selam cannot be explained with satellite formation via giant impacts (Durda et al., 2004(Durda et al., , 2007) ) alone.Yet, we can imagine a scenario where a combination of rotational failure and a grazing giant impact can lead to a similar type of debris disk investigated in this work and Agrusa et al. (2024), especially since Dinkinesh displays nonaxisymmetry and a significant concavity that could originate from a collision.

Particle size distribution, shape and resolution
One important caveat of the SFD in our model and the low resolution of our hand-off simulations is that it leads to unrealistically large boulders in the system and an absence of smaller grains.As shown in a series of recent publications related to the DART mission impact, the choice of boulder packing, cohesion and SFD highly affects the dynamics of rubble piles and whether or not they are disrupted or reshaped from high-energy collisions (Raducan et al., 2022;Raducan et al., 2024a,b).Combining observations from the DART mission (Li et al., 2023;Dotto et al., 2024) and SPH simulations of the DART impact, Raducan et al. (2024b) found that Dimorphos is expected to be a cohesionless rubble pile with a boulder packing of less than 40% on the surface and slightly below it.This accounts for boulders with  p > 2.5 m, meaning most of the mass is constituted by smaller grains and boulders, which aligns with observations of Dimorphos surface (Pajola et al., 2023).The volume of boulders in Dimorphos, their packing and SFD will be further constrained with the upcoming Hera mission (Michel et al., 2022), which will revisit the binary system satellite once more post DART impact and determine how the body was reshaped.
In our simulation, the presence of smaller particles has been omitted and their mass is simply accounted for by using larger clusters of particles represented by boulders with densities similar to the bulk density of the primary for computational purposes.Including particles with diameters less than 5 m would put a high computational strain on our simulations that cannot currently be handled with the resources at our disposal as the number of particles in each system would increase drastically, way above 10 4 .Nevertheless, results from Raducan et al. (2024b) and previous studies such as Zhang et al. (2017Zhang et al. ( , 2021) ) show that large differences in particle sizes affect the stability and structural dynamics of rubble pile asteroids.For our purposes, excluding the small meter-sized particles from our simulations can generate undesired behaviour in the contact mechanics and artificially high mechanical strength for the inner structure from particle interlocking between our large, highly angular particles.This becomes particularly important when considering the high rate of tidal interactions and mergers we observe.While our SFD is polydisperse for the standalone N-body simulations, it may be that smaller grains with diameters below 5 m could smooth out contacts between larger boulders, making the aggregates more fluid-like and thereby more susceptible to tidal disruption and deformation during mergers.Further, a more fluid-like behaviour would increase the time it takes for the newly formed aggregates to settle in a specific shape as particle interlocking effects would have less influence on the mechanical integrity of the body.With weaker mechanical strength, dynamical effects such as tidal forces, impacts, spin and self-gravity would have a more dominant role in the determination of the final configuration of the constituents.Keeping these potential caveats in mind, the specific shapes of the moons formed in our simulations may be unphysical, which is why we have opted to solely rely on DEEVE estimates for their geometrical properties.Moreover, another consequence of altering our SFD in this manner could be that it facilitates the creation of a distinct neck between two moonlets during the formation of a bilobate object.The low impact velocities in our model would allow for longer gravitational interactions prior to the contact and smaller loosely bound grains could flow towards the other body, creating the characteristic neck between them, seemingly present between the lobes of Selam (Levison et al., 2024).However, if rubble piles get their mechanical strength from the presence of larger boulders beneath their surfaces, interlocking effects between granular particles could potentially also provide the cohesive effect needed to keep structural integrity for the smaller constituent in contact binaries (Meyer and Scheeres, 2024).We aim to investigate the validity of our SFD and how it affects the stability of our aggregates in a follow-up study using higher-resolution simulations of select mergers from our simulations and tidal encounters with the primary.
With these model characteristics in mind, we also note that there are no evident dissimilarities in outcomes and dynamical behaviour of the hand-off and pure N-body systems in this work, which use drastically different SFDs.

Length of simulations and system stability
The heavy computational task of resolving contacts between polyhedra not only limits the number of particles we can use, but also how long we can realistically simulate a system, let alone several.Therefore, we had to settle for resolving each system a maximum of ten days post rotational failure.While we can capture and assess the main dynamical mechanisms involved in generating oblate and bilobate shapes, we cannot speak for the long-term stability of our aggregates and their orbits.There are many reasons a potential moon could end up being deformed or even disrupted.Since a majority of the debris disks in our systems have yet to be cleared of matter, further mergers with aggregates still remaining in the system that migrate over time and end up intercepting the orbit of the largest satellite can lead to further growth and reshaping.Furthermore, while unlikely, a lower mass object on an eccentric orbit that undergoes a close encounter with the primary body can be scattered, directly collide with the largest remnant and disrupt it or cause significant deformation due to the high-velocity impact.Orbit decay through long-term tidal interactions with the primary or scattering events with other aggregates could also lead to further mass shedding events and deformation due to tidal disruption.
Despite the mass available in our systems and the number of aggregates that interact gravitationally, most systems remain dynamically cold in the sense that there the orbital inclinations are kept small, causing a large majority of the close encounters to lead to mergers with low impact velocities.In consequence, changes in orbits mainly occur through energy dissipation from impacts, tidal effects and gravitational scattering with the primary.The mergers also lead to breaking of the often initially synchronous orbits of the largest remnant, leaving only a few of the satellites in such a state.One possible result is that the perturbed aggregate then enters a state of non-principal axis rotation, limiting or breaking any potential for migration via BYORP (Ćuk and Burns, 2005;Quillen et al., 2022).Depending on the inertia ratios of the satellite and its collisional history, it can also enter a state of tumbling (Agrusa et al., 2021).For bodies that retain their synchronous orbit, BYORP becomes important as a mechanism for migration and if multiple satellites remain in the system, mean motion resonances can further affect the orbital evolution.Long-term dynamical effects of binary asteroid systems are complex and cannot be captured in a meaningful way over the course of just a few days of dynamical evolution and we refer the interested reader to Agrusa et al. (2024) who analysed the orbital stability and rotational state of 96 simulations lasting 100 days for a similar, debris disk based binary asteroid formation model.

Comparison to Agrusa et al. (2024)
While there have been numerous previous studies on asteroid binary formation (Walsh et al., 2008(Walsh et al., , 2012;;Jacobson and Scheeres, 2011a;Tardivel et al., 2018;Hyodo and Sugiura, 2022;Madeira et al., 2023;Madeira and Charnoz, 2024), the work most related to ours is Agrusa et al. (2024) (henceforth referred to as A24).In this section, we summarise any differences and similarities between the work brought up throughout the text with the previous discussion in mind.The two models can consistently form oblate spheroid satellites in wide debris disks around similar primary bodies over timescales on the order of days, a few % of the cases in A24 and 3 out of 14 simulations for this work.Judging from the distribution of DEEVE axis ratios for the satellites formed in A24 (see their figure 12a), the simulations also seem to have resulted in a few bilobate objects.The dynamical mechanisms that determine the final shape are highly similar in both models as oblate spheroids form in cases where the largest remnant in each system undergoes a catastrophic re-shaping event due to a tidal encounter and/or a critical merger.Therefore, the overall trends for the dynamical evolution of the disks appear similar, independent of differing SFDs, particle shapes and contact physics.However, while it is difficult to argue that our implementation, on average, forms more oblate spheroids because of the low number of simulations we could carry out, it would align with previous points made during the discussion.Given how mergers and tidal encounters with the primary body are the main mechanisms that determine the shape of a satellite in this regime, this would mainly be attributed to the number of aggregates that simultaneously can exist in the system throughout the evolution of the disk, especially during the later stages.This number is governed by how efficiently a given disk can create aggregates and their mechanical strength after formation.With this in mind, we proceed to identify the main differences between the models in our work and A24 that could cause our aggregates to form more easily and survive for longer.
As mentioned in sections 2.2 and 2.5, the methods used to perform the spin-up simulations are based on two distinct codes and setups, which accounts for the variation in initial conditions for the subsequent dynamical evolution of the resulting debris disks.Crucially, the disks simulated in A24 were symmetric, less radially extended (by up to 1 main ), thinner and less massive (max 0.04 main ).Hence, it is clear that the type of mass-shedding event considered in A24 was less chaotic and more confined to areas around the equator of their primary.As discussed in section 5.1, the geometry of the debris disk appears to affect the number of aggregates available in the system, as the ability of an aggregate to survive the tidal forces depend on where they form in relation to the fluid Roche limit.The thickness of the disk should also affect its efficiency at forming satellites, as a thicker disk would let an aggregate accrete isotropically for longer before it starts accreting only within the plane, which is less efficient (analogous to 2D and 3D accretion rates in planet formation (e.g.Johansen and Lambrechts, 2017)).Thus, a satellite forming in a thick disk would be more efficient at clearing its vicinity of matter, reach higher masses and survive for longer before it ultimately merges with another aggregate.The longevity of a given aggregate is also dependent on its mechanical strength.Given that the use of non-spherical particles, as compared to using spherical particles with approximated granular behaviour, has been shown to increase the rigidity of rubble pile objects (Ferrari and Tanga, 2020) and ability to withstand tidal disruption (Marohnic et al., 2023), it is likely that our use of angular particles with rich geometric diversity and non-linear contact mechanics also contribute to a higher mechanical strength.This allows aggregates to more easily survive under the influence of strong tidal forces in our model.
Another main difference between the two studies comes from the disk masses considered.The satellites in A24 with masses similar to that of Dimorphos all occur in cases where the disk mass is 0.03 main or less, which offers an explanation for why our model generates more massive largest remnants.While the oblate spheroid satellites from our results all have lower masses than their more prolate counterparts (< 0.02 main ), the systems they form in are not fully resolved by the final time step as mentioned above.Therefore, it would be meaningful to run additional SPH spin-up simulations with a smaller effective friction angle to get initial conditions for less massive debris disks.Moreover, by running standalone N-body simulations with lower initial disk masses and altering their disk geometry to be more similar to the disks in A24, we would be able to further narrow down the driving factors that increase the frequency of systems with atypically shaped satellites.

Conclusions
In light of recent close-up observations of satellites in binary asteroid systems that diverge from the typical elongated shape, such as the oblate spheroid Dimorphos targeted by the DART mission (Rivkin et al., 2021), or the contact binary Selam of the Dinkinesh system detected by the Lucy mission (Levison et al., 2021(Levison et al., , 2024)), we have aimed to narrow down the mechanisms involved in the formation process of these objects in this work.In turn, we performed short-term dynamical evolution of circumasteroidal debris disks formed around rapidly spinning asteroids.Taking initial conditions obtained from SPH simulations of mass shedding through rotational failure, we used a newly developed tool based on geometrical algorithms to perform numerical hand-off to a granular N-body code that captures complex contact mechanics and gravitational interactions between irregularly shaped particles.Investigating two different systems, based on the asteroids Ryugu and Didymos, we carried out two hand-off simulations and five standalone N-body simulations for each case.Out of the 14 simulations, we end up with three systems with oblate spheroid secondaries and three systems with bilobate satellites by the final time steps.
We find that a radially extended disk wider than 2 primary radii favours mergers and deformation via tidal interactions with the primary body.Growth is rapid, hierarchical and primarily driven by mergers, allowing satellites between 0.014 main and 0.032 main to form in just a few days.As most moonlets initially end up in orbits near the fluid Roche limit, they are most often prolate and on synchronous orbits that are subsequently perturbed by mergers and close encounters with the main asteroid.In our model, only systems where the largest satellite undergoes mild or catastrophic disruption (losing less or more than half its mass, respectively) during close encounters with the primary, making it less prolate, are capable of forming satellites with Dimorphos-like shapes.
Analysis of critical mergers during our simulations shows that the impact velocities in this formation regime are below one mutual escape velocity, which is significantly lower compared to previous studies which have mainly been focused on merger-driven moon formation in planetary rings (e.g.Leleu et al., 2018).This hints at the existence of a different regime for mergers during the formation of asteroid moons which results in homogeneous bodies without the significant ridges and deformations caused by high-velocity impacts.The low-energy mergers facilitate the formation of bilobate satellites with two lobes connected by a thin neck, which we expect to find when performing higher-resolution simulations of some impacts.Moreover, we find that the initial shape and the orientation of the longest principal axis of the most massive body in the merger highly influence the post-merger shape.Given the bodies are initially prolate before impact, mergers where the relative velocity vector is parallel with the shortest axis lead to more oblate shapes while mergers where the same vector is more closely aligned with the longest axis result in more prolate bodies and bilobate objects.
Due to the heavy computational task of resolving contacts between irregularly shaped particles, we had to settle for a shorter simulation time than desired of at most 10 days and an SFD not allowing for particles smaller than 5 m in diameter.In turn, we cannot draw any conclusions regarding the stability of our systems and the potential presence of higher-order dynamical effects such as BYORP, mean motion resonance between multiple satellites and tumbling.Nevertheless, our results provide valuable insight into the formation mechanisms that govern deformation of the initially prolate bodies asteroid satellites that make up a majority of the population (Pravec et al., 2016).Judging from the chaotic formation pathways in our simulations, oblate spheroids and bilobate satellites may be more prevalent in asteroid binaries than indicated by lightcurve measurements, which are biased against detecting satellite shapes that are more oblate (Pravec et al., 2016;Agrusa et al., 2024).
While computationally expensive, the added internal mechanic strength from the contact mechanics in our N-body simulations caused by particle-particle interlocking, dilatancy and complex random packing of the irregularly shaped grains better our understanding of how granular aggregate satellites behave during mergers and close encounters with their primaries.When complete tidal disruption is more difficult, more aggregates will survive the earliest stages of formation, allowing for more mergers and reshaping events.However, recent studies such as Raducan et al. (2024b) have emphasised that the structural integrity of aggregates is highly dependent on the distribution between large boulders and smaller particles with diameters less than 2.5 m.Furthermore, the observed boulder diameters on the surface of Dimorphos go down to 1 m and up to 16 m, while we include much larger particles serving as clusters of smaller constituents.Given that the presence of smaller particles can smooth out contacts during deformation events such as mergers and tidal interactions, making the body behave more like a fluid, we recognise that our aggregates may be too rigid and will further investigate the effect of particle SFD choice for our implementation.With upcoming missions such as Hera (Michel et al., 2022), which will send a spacecraft to revisit the Didymos system in 2027 after the DART impact and determine the resulting shape of Dimorphos and its internal structure, the boulder volume fraction and SFD can be further constrained, allowing us to significantly improve our formation model.

Figure 1 :Figure 2 :
Figure1: Left: the position of the SPH particles at the time of the hand-off for a symmetric debris disk formed from a Ryugu-like primary.Right: the final output of clusterfinder.The algorithm has generated -shapes and convex hulls with volumes matching the densities and masses of SPH particle clusters.Zooming in will show the resulting 3D shapes with black grid lines marking the edges that separate the faces in blue.

Figure 3 :
Figure3: This figure depicts three different SFDs mentioned in the text, as well as the probability density of the particle SFD obtained from the hand-off f60_R2 as a histogram with 20 bins, which is our example simulation discussed in section 3.2.The SFDs have been generated using the previously defined exponential function in equation 1.The blue curve is our closest approximation to observed sizes on Dimorphos surface, while the green matches the particles from the f60_R2 handoff.The orange curve shows the much flatter distribution we use to randomly generate particles for our standalone N-body simulations.

Figure 4 :
Figure4: Snapshots from the hand-off simulation f60_R2 where a debris disk around a Ryugu-like primary quickly forms a satellite over the course of 73 hours.The final shape of the moon is largely determined by the merger between the two largest remnants in the system, occurring during the period separating the two final snapshots at 64.2 and 78.3 h.The two prolate aggregates collide side-to-side and form a more oblate shape.

Figure 5 :
Figure 5: Growth over time after the formation of the debris disk, , for the largest remnant in the system with mass  lr , scaled by the mass of the primary body  main .The dashed black line indicates the target mass ratio between Dimorphos and Didymos, assuming equal densities for the two bodies.

Figure 8 :
Figure 8: The final shape for the 0.028 main aggregate consisting of 869 grains that formed during the dynamical evolution after mass shedding by rotational failure for a Ryugu-like body in simulation f60_R2.Left: the aggregate as seen from a top-down perspective.Right: an edge-on view of the final satellite.

Figure 9 :
Figure 9: The same plot as in figure 5 but for the case of five separate standalone N-body simulations of a Ryugu-like system.

Figure 12 :
Figure 12: The change in principal axis ratios ∕ for the corresponding DEEVE shape of each largest remnant over time for the five Ryugu pure N-body simulations.Due to excessive noise, the data has been smoothed using a uniform one-dimensional filter to facilitate interpretation.Note that this may skew the individual data points from their actual value.The evolution prior to 25 h has been omitted due to the naturally chaotic behaviour of the formation process.Comparing the tracks to those representing the largest aggregate mass over time in figure 9, there are clear overlaps between large mass shedding or accretion events and spikes in ∕ values.

Figure 13 :
Figure 13: The final shape of the largest satellite in run 4 of the standalone Ryugu N-body simulations.Left: the configuration of individual grains in the aggregate at the final time step of the simulation.Right: the corresponding -wrap shape for the same configuration of particles.

Figure 14 :
Figure14: Left: the position of the SPH particles at the time of the hand-off for an nonaxisymmetric debris disk formed from a Didymos-like primary.Right: the final output of clusterfinder.The algorithm has generated -shapes and convex hulls with volumes matching the densities and masses of SPH particle clusters.Zooming in will show the resulting 3D shapes with black grid lines marking the edges that separate the faces in blue.

Figure 15 :
Figure 15: Left: a histogram showing the number density of the disk obtained after performing the hand-off for the simulation f60_D2 in the -plane.Corresponds to the same distribution seen in figure 14.Can be compared to the Ryugu example case with a more symmetric distribution in figure 2. Right: the same metric but for the -plane.Note the differing scales between the and -axis of the figure.

Figure 16 :
Figure16: Growth over time after the formation of the debris disk, , for the largest remnant in the system with mass  lr , scaled by the mass of the primary body  main .The dashed black line indicates the target mass ratio between Dimorphos and Didymos, assuming equal densities for the two bodies.

Figure 18 :
Figure 18: The same plot as in figure 16 but for the case of five separate standalone N-body simulations of a Didymos-like system.

Figure 20 :
Figure20:The final shapes at the end of each simulation for the Didymos-like scenario.Each ellipsoidal axis has been evaluated from the corresponding dynamically equivalent equalvolume ellipsoid to the moment of inertia of each aggregate.Note that only the coloured markers represent bodies that have a large enough mass and are situated within 5 km of the primary body.They can be compared to the final masses and positions in the different systems from figure19.We note that these values are likely slightly larger than what we would obtain from the principal axes corresponding to the actual physical extents of the bodies(Agrusa et al., 2024).

Figure 21 :
Figure21: The same plot as in 17 but for the standalone N-body simulations for a Didymos-like primary.The graphs can be compared to the change in mass over time in figure18to identify major mass shedding events and mergers.
Figure22: The DEEVE axis ratio ∕ value before (grey) and after (coloured) the critical merger which determines the final shape of the considered secondaries at the final time step of our simulations, compared with four different impact properties.Arrows pointing upwards generally mean that the body became more prolate (elongated) after the merger, while an arrows pointing downwards indicates that a deformation to a more oblate shape.Note that we fix the definition of each axis as  or  before the merger, which is why the ∕ axis ratio can be below unity for some cases (e.g.Didymos run 1).Upper left: the ∕ axis ratios and the corresponding impact velocity divided by the mutual escape velocity.Upper right: same axis ratio change plotted against the impact angle between the two velocity vectors.Lower left: the axis ratio compared to the mass ratio of the impactor and target.Lower right: the distance of the merger from the barycenter of the primary in terms of main radii.

Table 1
The DEEVE axis ratio ∕ value before (grey) and after (coloured) the critical merger which determines the final shape of the considered secondaries at the final time step of our simulations, compared against  p which is the angle between the longest principal axis and the relative velocity vector for the merger.A short-axis merger corresponds to  p ∼ 90 • while a long-axis merger would mean a value close to 0 or 180 • .