Direct N-body Simulations of Satellite Formation around Small Asteroids: Insights from DART’s Encounter with the Didymos System

We explore binary asteroid formation by spin-up and rotational disruption considering the NASA DART mission's encounter with the Didymos–Dimorphos binary, which was the first small binary visited by a spacecraft. Using a suite of N-body simulations, we follow the gravitational accumulation of a satellite from meter-sized particles following a mass-shedding event from a rapidly rotating primary. The satellite’s formation is chaotic, as it undergoes a series of collisions, mergers, and close gravitational encounters with other moonlets, leading to a wide range of outcomes in terms of the satellite's mass, shape, orbit, and rotation state. We find that a Dimorphos-like satellite can form rapidly, in a matter of days, following a realistic mass-shedding event in which only ∼2%–3% of the primary's mass is shed. Satellites can form in synchronous rotation due to their formation near the Roche limit. There is a strong preference for forming prolate (elongated) satellites, although some simulations result in oblate spheroids like Dimorphos. The distribution of simulated secondary shapes is broadly consistent with other binary systems measured through radar or lightcurves. Unless Dimorphos's shape is an outlier, and considering the observational bias against lightcurve-based determination of secondary elongations for oblate bodies, we suggest there could be a significant population of oblate secondaries. If these satellites initially form with elongated shapes, a yet-unidentified pathway is needed to explain how they become oblate. Finally, we show that this chaotic formation pathway occasionally forms asteroid pairs and stable triples, including coorbital satellites and satellites in mean-motion resonances.


INTRODUCTION
Binaries make up ∼15% of the near-Earth asteroid (NEA) population (Margot et al. 2002;Pravec et al. 2006).This fraction increases to ∼65% for fast-rotators greater than ∼300 m in diameter (Pravec et al. 2006).Given the relatively short median dynamical lifetimes of NEAs of about 10 Myr (Gladman et al. 2000), this high binary fraction implies an efficient formation mechanism that can maintain a steady-state population or formation before exiting the main asteroid belt.Among NEA binaries, the primary component is typically small, less than ∼5 km and with a moderately sized secondary, usually less than ∼60% of the primary's diameter.Furthermore, the two bodies are typically on close orbits, separated by less than ∼10 primary radii.Typically, the primary is rapidly rotating and has a near-spherical shape with an equatorial bulge (Pravec et al. 2006;Pravec & Harris 2007;Benner et al. 2015).In addition, the secondaries of these systems frequently have an elongated shape and are typically in synchronous rotation when the secondary is on a close, near-circular orbit (Pravec et al. 2016).The specific angular momentum of these systems is often close to the angular momentum of a sphere having the same total mass rotating at its critical disruption limit, indicating that the creation mechanism of these binaries must be related to some sort of fission or mass-shedding process (Pravec & Harris 2007;Pravec et al. 2010).
The Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect is a process in which the absorption and re-emission of solar radiation imparts a small torque on irregularly-shaped bodies (Rubincam 2000).The strength of the YORP effect is highly dependent on the size and shape of the body and its heliocentric distance, and is thought to be the dominant mechanism that spins-up small asteroids leading to rotational disruption and the formation of a binary (Vokrouhlický & Čapek 2002;Bottke et al. 2002).For a review of YORP and the formation of small binary asteroids, see Walsh & Jacobson (2015).
The recent Double Asteroid Redirection Test (DART) mission was the first spacecraft to visit (albeit briefly) a small binary asteroid (Rivkin et al. 2021).The primary, Didymos, has a flattened top-shape with a volume-equivalent diameter of ∼760 m and a short spin period of 2.26 h (Naidu et al. 2020;Daly et al. 2023b).The satellite, Dimorphos, has an approximately circular orbit with a semimajor axis of ∼3 primary radii (Naidu et al. 2022;Scheirich & Pravec 2022).Dimorphos was thought to be in synchronous rotation prior to DART's kinetic impact (Richardson et al. 2022), although this could not be determined directly.Both bodies have surfaces consistent with being gravitational aggregates or "rubble piles" (Barnouin et al. 2023), although any direct measurements of their interiors will have to wait for the arrival of the European Space Agency's Hera mission in 2027 (Michel et al. 2022).Dimorphos has a volume-equivalent diameter of ∼150 m, corresponding to a secondary-to-primary size ratio of ∼0.2.If Didymos and Dimorphos have equal bulk densities (which not known), then they would have a mass ratio of ∼0.01.Dimorphos appears to have a remarkably oblate shape, in which its semi-major axis, a, and semi-intermediate axis, b, are nearly identical in length and roughly ∼1.5 times longer than its semi-minor axis, c (i.e, principal rotation axis) (Daly et al. 2023b).Additionally, the presence of rocks perched on boulders with slopes no higher than ∼35 • suggests that Dimorphos's surface may be cohesionless with a friction angle of ∼35 • (Barnouin et al. 2023).A ∼35 • friction angle is also derived from an independent method based on the angularity of boulders on Dimorphos's surface (Robin et al. 2023).Although there could be some strength at depth, some recent models of DART's impact are consistent with a cohesionless interior (Raducan et al. 2023).
The shape of Dimorphos came as a surprise, as previously observed systems (either via radar or lightcurve) have found that the secondary tends to have a more prolate, or elongated, shape in which a > b > c, where a, b, and c are the body's respective semimajor, semi-intermediate, and semi-minor axes.Based on Dimorphos's physical extents (which are well-fit to a uniform ellipsoid), a preliminary estimate puts Dimorphos's elongation at a/b = 1.02 ± 0.03 (Daly et al. 2023b) which has since been updated to a/b = 1.06 ± 0.03 (Daly et al. 2023a).Nevertheless, no previously published studies employing lightcurve-or radar-based shape determination of similar binary systems have found a satellite with a/b ≲ 1.1 (Ostro et al. 2006;Naidu et al. 2015;Becker et al. 2015;Pravec et al. 2016).The principal motivation for this study was to understand how Dimorphos might have formed with its present shape, but as we will show later, the results are broadly applicable to all small binaries formed by rotational disruption.
In Sections 1.1 and 1.2, we summarize previous work on binary formation and YORP spin-up.Then we introduce pkdgrav in Section 2 and show a full end-to-end example simulation demonstrating the formation of a binary system via a single YORP-driven mass-shedding event in Section 2.1.Then, based on this simulation and constraints derived by DART, we introduce the simulation set-up and initial conditions used for this study in Section 2.2.The results from the numerical simulations are shown in Section 3 followed by conclusions in Section 4.

Previous work
There is general agreement that many small binary asteroids are likely created via YORP-induced spin-up in which a satellite forms when the primary exceeds its critical spin limit.However, there is still some debate about precisely how the binary system forms, which we address here.Walsh et al. (2008) modeled the formation of a binary using a rubble-pile asteroid model consisting of thousands of spherical particles in which angular momentum is slowly added to the system as a proxy for YORP-induced spin-up.Through the spin-up process, the primary reshapes, forming an equatorial bulge, and the secondary is gradually built in orbit via gravitational accumulation of material shed from the primary's equator.This idea was revisited in Walsh et al. (2012) with a broader simulation suite at higher resolution, finding that that the resulting top-shaped primary and secondary properties are broadly consistent with the observed population.
However, this model suffers from a few issues.First, it has been shown that the magnitude and direction of the YORP torque is highly-sensitive to small changes in the primary's shape, meaning that each landslides, mass shedding, and natural impacts could change the strength and direction of the YORP effect.This could make YORP spin-up effectively a random walk process and may significantly increase the amount of time required to build a secondary in orbit (Statler 2009;Cotto-Figueroa et al. 2015;Zhou et al. 2022).This effect may significantly decrease the efficiency of such a formation process, making a scenario in which the secondary forms from a single rotational disruption event more attractive 1 .On a related note, Jacobson & Scheeres (2011) point out that a gradual process of satellite formation suffers from the fact that single particles can escape the system before the satellite is fully formed.The argument is that escape can occur before the next shedding event because the system has positive free energy, and the single (or multiple) shed particles would escape before the next YORP cycle.Second, in the Walsh et al. (2008Walsh et al. ( , 2012) ) works, the primary was initialized in a hexagonal close-packed (HCP) arrangement and particle contacts were handled using the Hard-sphere Discrete Element Method (HSDEM), which is ill-suited to simulate particles undergoing persistent contact with multiple other particles (see Murdoch et al. 2015;Sánchez & Scheeres 2012, 2016).Although this model captures the general process of binary formation, the HCP packing and HSDEM contact physics may affect the precise details of the mass shedding and satellite formation.Jacobson & Scheeres (2011) presented an alternative theory for binary formation in which the secondary arises from a single large rotational disruption event, dubbed "rotational fission" (Scheeres 2007a).There seems to be some confusion about these two ideas in the literature, and we attempt to point out here that the "gravitational accumulation" idea of Walsh et al. (2008) and the "rotational fission" idea of Jacobson & Scheeres (2011) are similar, yet fundamentally distinct.Jacobson & Scheeres (2011) present a model that follows the spin-orbit coupled dynamical evolution of a binary system following a fission event.Their model assumes the initial rubble pile consists of two ellipsoidal components (that will later become the primary and secondary).Angular momentum is slowly added to the single body and at fast spin rates, the long axes of the two constituent ellipsoids become aligned while still at rest on one another, owing to this being the only stable equilibrium configuration for two ellipsoids (Scheeres 2007a).When the critical spin limit is exceeded, the bodies then "fission" and begin to orbit each other.Jacobson & Scheeres (2011) integrate the fully-coupled spin and orbital evolution of the binary system, accounting for their nonspherical gravitational potentials along with a treatment for tidal evolution.For mass ratios < 0.2, they argue that the postfission system has positive free energy (i.e., the sum of the bodies' kinetic energies and the mutual potential), meaning that there are one of two possible outcomes: the secondary must either escape the system or it must fission itself to form (temporarily) a triple system.These secondary fissions occur when the secondary is gravitationally torqued by the primary such that it exceeds its stability limit and is then split into two separate ellipsoids (whose shapes are generated randomly).In most of their simulations the third body re-impacts the primary, but it can also escape the system or re-impact the secondary.This ongoing fission process, along with tidal evolution, will eventually lead to a state with negative free energy, ultimately allowing the binary system to be stable.
This model has the advantage that it only requires a single rotational disruption event, thus avoiding the "stochastic YORP" problem.This model broadly reproduces the characteristics of the NEA population, including the existence of asteroid pairs.A follow-up study that includes a more sophisticated asteroid population model shows that the asteroid fission theory can reproduce many of the characteristics of binary asteroids (Jacobson et al. 2016).Recently, 1 Recent work shows that all twelve of the available and reliable YORP detections show that that the asteroid's spin rate is increasing in time ( Ďurech et al. 2023).This could due to YORP having an underlying preference to increase the spin rates of asteroids rather than decrease them (Golubov & Krugly 2012).If there really is an underlying preference for spin-up rather than spin-down, then the idea of stochastic YORP does not present a significant issue to binary formation models.
the rotational fission model has been extended to include 3D dynamics, and many of the conclusions still largely hold true (Boldrin et al. 2016;Davis & Scheeres 2020;Ho et al. 2022).
The Jacobson & Scheeres (2011) suffers from one potential weakness; in order to allow for secondary fission events, which are necessary to achieve stability, they invoke that the secondary itself is a rubble pile.This is of course a sensible assumption, however, it presents a major issue; at the very moment of the initial fission event, when the secondary (which is itself a rubble pile!) detaches from the primary, it must tidally disrupt, as it is lying within the primary's Roche limit (Holsapple & Michel 2006;Sharma 2009) rather than continue to evolve as a coherent body.This dilemma can be avoided in cases where the secondary's bulk density is much higher than the primary's such that the Roche limit lies within the primary itself, or if the secondary has a small amount of cohesion to prevent a tidal disruption (e.g.Holsapple & Michel 2008;Sánchez & Scheeres 2016;Tardivel et al. 2018), or if the secondary has a bilobate shape in which the neck region could fail before the rest of the body (Hirabayashi & Scheeres 2019), or some combination of all three.However, neither of these three contingencies seem viable for Dimorphos to form purely by fission if it is a rubble pile with little-to-no cohesion and a friction angle of ∼35 • (Barnouin et al. 2023;Raducan et al. 2023).If such a body were to rotationally fission from Didymos's equator, it would immediately undergo tidal disruption (e.g.Agrusa et al. 2022).
In this paper, we present an alternate formation path that incorporates the ideas of both the Walsh et al. ( 2008) and Jacobson & Scheeres (2011) models.In our model, the secondary forms via gravitational reaccumulation, similar to the Walsh et al. (2008) theory, however, all the required mass is shed in a single impulsive event, like Jacobson & Scheeres (2011).We will demonstrate that this simple model of accumulation from debris shed from a single rotational disruption event can lead to stable binary systems and secondaries having a wide range of shapes.
Recently, Madeira et al. (2023) considered the gravitational accumulation of Dimorphos from debris produced by Didymos using a 1D ring-satellite model (HYDRORINGS) (Charnoz et al. 2010(Charnoz et al. , 2011;;Salmon et al. 2010).In this model, material migrates outwards via viscous spreading until it reaches the Roche limit, at which point it is immediately converted to a moonlet.Assuming the initial debris ring is narrowly confined at the primary's surface, they find that it takes ∼1 yr for the ring to spread to the Roche limit, after which Dimorphos would accrete rapidly, reaching ∼90% of its current mass within days.This process requires ∼25% of Didymos's mass to be put into orbit to form a satellite with ∼1% of Didymos's mass.Once Dimorphos forms, they estimate that the ring would take thousands of years to disappear.In terms of timescales and the required amount of mass, the predictions of this model differ significantly than what we present in this work, most likely due to different assumptions in the initial conditions and necessary simplifications of the 1D model.A significant body of literature suggests that a rubble pile undergoing rotational failure and mass shedding, due to its nonzero friction and/or cohesion, will shed material onto much wider initial orbits (e.g.Hirabayashi et al. 2015;Sánchez & Scheeres 2016, 2020;Zhang et al. 2018;Hyodo & Sugiura 2022) rather than within a narrowly confined region, including several studies focused on Didymos itself (e.g., Yu et al. 2018;Zhang et al. 2017Zhang et al. , 2021)).This means that satellites can start forming much quicker, without having to wait for a ring to spread to the Roche limit.Also, the Madeira et al. (2023) study assumes that all ring particles are 1 m in diameter, whereas Dimorphos's contains boulders at least as large as 16 m (Pajola et al. 2023a).The presence of larger boulders would speed up the accumulation process, due to their higher mass and larger collision cross section.Finally, this model neglects gravity between moonlets, assuming that they undergo perfect mergers whenever they enter their mutual Hill spheres.In this work, we will show that this is often not the case, and that the merger process of moonlets is highly chaotic, leading to a wide range of outcomes.

YORP spin-up and mass shedding
In this study, the considered binary formation scenario is based on the idea that the progenitor body sheds substantial surface materials in a relatively short timescale.However, this mass shedding failure behavior is not unconditional.Prior theoretical work, relying on static stress analyses, suggests that a homogeneous ellipsoidal body would first structurally fail internally at its center instead of at the surface during spin up (Holsapple 2010;Hirabayashi 2015).This internal failure would lead to internal deformation, suppressing surface landslides and mass shedding.Modeling of a rubble-pile body's spin-up using the soft-sphere discrete element method (SSDEM) supports these static analyses, showing that homogeneous cohesionless bodies could indeed fail through internal deformation (Sánchez & Scheeres 2012;Hirabayashi et al. 2015) and homogeneous cohesive bodies could even fail through global tensile disruption (Sánchez & Scheeres 2016;Zhang et al. 2018).
For mass shedding to occur on the surface, the body's interior must be mechanically stronger and/or denser than the exterior layer.Research by Sánchez & Scheeres (2018) using SSDEM simulations on a spherical rubble pile with an internal core occupying 70% of its radius demonstrates that the interior needs to be 3 to 4 times stronger than the shell to prevent internal failure.Similarly, Ferrari & Tanga (2022) use a polyhedron discrete element method to show that a rigid core with a volume fraction exceeding 25% can also facilitate mass shedding, as shown for the case of Didymos.By employing higher particle resolution (∼10 5 particles) and a highly polydisperse particle size distribution, the SSDEM spin-up simulations of Zhang et al. (2017Zhang et al. ( , 2021) ) reveal that the large particle size difference can reduce the internal porosity and increase the mobility of surface regolith in a rubble-pile.The resultant small scale heterogeneity could trigger surface failure via mass shedding, even if mechanical interactions at the particle level are uniform throughout the body.A similar effect is produced by the interactions between non-spherical particles with irregular shape, where geometrical interlocking mechanism adds mechanical strength to the inner structure (Ferrari & Tanga 2020).In all cases, maintaining some internal shear strength and a low surface cohesion are crucial to initiate surface mass shedding (Sánchez & Scheeres 2020;Zhang et al. 2022).
Beyond considering the internal structure, the evolution of a rubble-pile body after reaching its spin limit is also heavily influenced by numerous other factors, such as shape and surface morphology (Hirabayashi & Scheeres 2019;Hirabayashi et al. 2020;Zhang et al. 2022).This complexity implies that the surface failure conditions could vary across different bodies, which determines why certain asteroids evolved into binary systems while others didn't, and could play some role in the apparent deficit of small binaries and pairs among primitive asteroid types relative to silicate-rich bodies (Minker & Carry 2023).In Section 2.1, to validate the assumed binary formation scenario for Didymos, we will present an SSDEM simulation example to showcase the feasibility of Didymos' mass shedding failure behavior under reasonable conditions.
We finally note that the mass shedding doesn't have to occur solely as a result of YORP spin up.Once YORP spins up a body near its critical limit, a natural impact could trigger the mass shedding event, for example.In fact, this study is somewhat agnostic to the exact process causing the mass shedding, we merely require that the mass shedding occur all at once such that enough material is put into orbit to allow a rapid reaccretion process.

METHODOLOGY
We use pkdgrav, a gravitational N -body tree code to model the accumulation of a rubble-pile satellite (Richardson et al. 2000;Stadel 2001).Particle contacts are handled using the soft-sphere discrete element method (SSDEM), in which a Hooke's Law spring provides a normal force between overlapping particles (Schwartz et al. 2012).Parameters such as the spring constant and coefficients of static, twisting, and rolling friction and cohesion can be adjusted to represent desired material properties.We refer the reader to Zhang et al. (2017) for a detailed description of the model.In short, k N and k T are the two spring constants that control the particle's stiffness in the normal and tangential directions respectively, ϵ N and ϵ T are the coefficients of restitution in the normal and tangential direction for controlling energy dissipation, µ S , µ R , and µ T are the coefficients of static, rolling, and twisting friction.Finally, the shape parameter β is used to approximate the non-spherical nature of real particles by statistically defining their contact area, enabling the calculation of the associated rotational resistance.
The normal spring constant (k N ) and timestep are selected such that typical particle overlaps don't exceed ∼1% of a particles radius and that ∼30 timesteps take place over the course of a collision (Schwartz et al. 2012).This ensures that particle contacts are resolved properly and that particles do not undergo a nonphysical level of interpenetration.Based on our typical particle mass, sizes, and collision speeds this corresponds to a k N ∼ 4 × 10 6 kg s −2 and a punishingly small timestep of ∼0.15 s.Following common practice, we set k T = 2 7 k N (Walton 1995;Schwartz et al. 2012).In all simulations presented here, we set the restitution coefficients to ϵ N = ϵ T = 0.55, a typical value for rocky material (Chau et al. 2002).µ R and µ T are set to 1.05 and 1.3, respectively, to represent the rough surfaces of medium hardness rocks (Jiang et al. 2015), leaving µ S and β as the only two parameters that are varied to achieve different friction angles.These values and their corresponding friction angle are provided in Table 1.

Satellite formation via mass shedding
Here, we show an example simulation to test Didymos's structural failure behaviors at its spin limit.Didymos is modeled as a granular aggregate composed of 87,635 spheres with radii ranging from about 4 to 16 m.The particle size distribution follows a cumulative power-law distribution with an exponent of −2.5, aligning with the boulder size frequency distribution with similar radii found on Dimorphos, as observed in the images taken by the camera onboard the DART spacecraft (Pajola et al. 2023a).The granular aggregate was configured based on the pre-impact Didymos shape model (Barnouin et al. 2023) to ensure accurate representation of Didymos.To allow the initiation of surface mass shedding and maintain the stability of Didymos at its current spin period of 2.26 h, the shape parameter β is adopted to be 0.8, and µ S is taken to be 1.0, representing a material internal friction angle of ∼40 • (Zhang et al. 2022).This friction angle is within the typical range of compacted dry sand, i.e., 33 • -43 • (Bareither et al. 2008).To resolve the quasi-static mechanical contact between particles, we adopted a smaller timestep of ∼0.02 s and a larger k N ∼8 × 10 7 kg s −2 .The bulk density of the body is set to 2.7 g cm −3 , consistent with Didymos' latest bulk density estimate constrained by the updated shape model and the binary separation, i.e., 2.760 ± .130g cm −3 (Naidu et al. 2023) 2 .Our numerical investigation found that cohesion is no longer required to maintain the bulk structural stability of Didymos at its current spin state for a bulk density of ≳2.7 g cm −3 , and, therefore, the interparticle cohesive strength was set to zero 3 .All other parameters are the same as those introduced in the previous section.
Didymos is quasi-statically spun-up (as a proxy for YORP) until structural failure is detected, which we define by the body's longest dimension changing by more than 1% relative to the starting shape.Then the spin-up procedure is halted and the granular system evolves purely under its own self gravity.The results show that the primary structurally failed at a spin period of ∼2.2596 h, where it then shed surface particles from mid-to-low latitudes, putting ∼3% of its total mass onto low inclination orbits.Much of this material rapidly clumps into moonlets, and a Dimorphosmass satellite is formed within days following a series of moonlet mergers.As a result of the conservation of angular momentum, the primary's spin period drops to ∼2.5 h by the end of the simulation due to mass shedding and reshaping.Material that falls back onto Didymos preferentially lands on the equator, which contributes to forming Didymos's equatorial ridge, similar to the process demonstrated by Hyodo & Sugiura (2022).The present-day shape of Didymos is the subject of other ongoing and future studies (Barnouin et al. 2023, Zhang et al., in prep).
Snapshots of this simulation are shown in Fig. 1 along with a time-series plot of the satellite's shape, mass, orbit, and attitude in Fig. 2. In this simulation, we form an approximately Dimorphos-mass satellite in a matter of only a few days.The satellite is also relatively oblate compared to the measured shapes of other binary systems, having axis ratios of a/b∼1.15 and b/c∼1.5, based on the satellites dynamically equivalent equal-volume ellipsoid (DEEVE) 4 .The satellite is initially in synchronous rotation with a libration amplitude of ∼45 • .However, this tidally-locked state is broken when the satellite has a close encounter with a smaller moonlet, sending it inwards where it undergoes a partial tidal disruption followed by an immediate merger with an another moonlet.This sequence of events breaks the satellite's synchronous rotation state and happens to lead to a relatively oblate shape.We encourage the reader to view the provided movie of this simulation.
This simulation robustly demonstrates that a Dimorphos-like satellite can form rapidly from a single rotational disruption event.In addition, the satellite's series of close encounters and collisions with other moonlets demonstrates that this process is highly chaotic.An infinitesimal change to the initial conditions of this simulation could lead to a very different outcome in terms of the satellite's physical and dynamical properties.However, due to the computational expense of these high resolution, full spin-up-to-satellite-formation simulations, they are impractical for longer simulations as well as studying outcomes statistically, which is the focus of the rest of this study.Before proceeding to the rest of this study, we note that the simulation presented in this subsection is a new result in its own right, although its primary purpose is to motivate the initial conditions used in the rest of this study.This simulation also highlights some key differences between our study and recently published papers on binary formation, which we address here.
2 Due to uncertainties in the size and volume of Dimorphos, which has some degeneracy in the body separation, there is significant uncertainty in Didymos's bulk density.The formal uncertainty of ±.130 g cm −3 is likely an underestimate and a realistic uncertainty might be larger. 3We also performed simulations using a friction angle of 35 • (µ S = 0.6, β = 0.5), and found that a bulk cohesive strength of about 8 Pa is needed for the structural stability.The mass-shedding structural failure behavior at faster spin is similar to the case presented here. 4For a body of mass m and principal moments of inertia A, B, C, its corresponding dynamically equivalent equal-volume ellipsoid axis lengths a, b, c are given by the following relations:  The satellite is initially tidally locked to the primary (starting at ∼1.3 d), but then undergoes a partial tidal disruption followed by a merger with a moonlet, causing it to break from synchronous rotation (at ∼2.5 d).The satellite has an oblate shape, similar to Dimorphos.However, the satellite's shape, mass, orbit and rotational state will continue to change as it continues to accrete more material.
Comparison with Sugiura et al. (2021) and Hyodo & Sugiura (2022) -Recently, Sugiura et al. (2021) studied the shape evolution of a rubble pile under YORP spinup using smoothed-particle hydrodynamics (SPH) simulations.In contrast to DEM simulations where particles represent discrete objects (i.e., rocks, boulders, etc.), SPH is used to simulate continuum mechanics where the particles sample local quantities such as density, internal energy, and pressure.Sugiura et al. (2021) find that spin-up can result in a "top-shaped" body for friction angles exceeding 70 • .However, it is challenging to compare this study to the result obtained due to differences in initial conditions, material models, and spin-up procedure.First, the work of Sugiura et al. (2021) starts from spherical shapes, whereas this simulation starts from a Didymos-like shape.Second, these simulations consider friction angles exceeding those of terrestrial and lunar granular material (typically ∼ 35−45 • , e.g., Bareither et al. 2008;Al-Hashemi & Al-Amoudi 2018;Mitchell et al. 1972) as well as friction angle estimates of recently visited asteroid surfaces (also on the order of ∼35 • , e.g., Fujiwara et al. 2006;Watanabe et al. 2019;Barnouin et al. 2019;Barnouin et al. 2023;Robin et al. 2023).This is because Sugiura et al. ( 2021) consider a single "effective friction angle" that accounts for both friction and cohesion.Cohesion is defined as the shear strength at zero pressure, where the shear strength of a granular material can be written as Y = tan(ϕ)p + c, where ϕ is the friction angle, p is the confining pressure, and c is the cohesion.Sugiura et al. (2021) instead adopt a material mod4el of the form Y = tan(ϕ)p where ϕ is now an "effective friction angle" that encompasses the effects due to cohesion, which justifies their choice to explore values of ϕ that are significantly higher than real granular materials.With this material model, we can see that as the confining pressure goes to zero (i.e., near the surface), the shear strength goes to zero, which effectively results in the material having zero cohesion.Therefore, we suspect that the failure mechanisms observed by Sugiura et al. (2021) are a result of the rubble pile implicitly having a cohesionless upper layer with a stronger internal core.Ryugu and Bennu seem to have sub-Pa levels of cohesion at their surfaces and potentially some strength at depth, so this implicit assumption by Sugiura et al. (2021) isn't unreasonable (Arakawa et al. 2020;Jutzi et al. 2022;Walsh et al. 2022;Barnouin et al. 2019).However, recent work by Zhang et al. (2022) demonstrates that rubble pile failure mechanisms are highly sensitive to small changes in cohesion at a fixed friction angle, and to small changes in friction at constant cohesion, indicating that friction and cohesion cannot be combined into a single parameter and should be treated separately.
Finally, the spin-up procedure used here differs from Sugiura et al. (2021) in a critical way.Both methods apply a similar angular acceleration to the rubble piles (∼10 −10 rad s −2 ), which ensures that the artificially induced Euler acceleration is always negligible compared to the centrifugal acceleration.However, Sugiura et al. (2021) spin up the rubble pile until 1% of its mass is put into orbit, whereas this study halts the spin-up process at the moment of failure to ensure a realistic mass shedding process.Although the total angular momentum added during the mass-shedding process by Sugiura et al. ( 2021) is small, this can make a critical difference for a rubble pile on the edge of stability.As a result, Sugiura et al. (2021) typically find that ∼10% of the primary's mass ends up going into orbit following the spin-up process.Hyodo & Sugiura (2022) then track the formation of moons after such a mass-shedding event, where they find similar formation timescales to what is presented in this study.However, due to their simulations starting with 10 − 20% of the primary's mass in a debris disk, their simulations result in a system of satellites with a combined mass of ∼5% of the primary's mass.Although many binary asteroids have similarly high mass ratios (on the order of 5 − 10%), this model does not produce binary systems with smaller satellites, such as the Didymos system, where the mass ratio is ∼1% (Pravec et al. 2016(Pravec et al. , 2019)).
Comparison with Madeira et al. (2023) -In comparison to the simplified 1D model of Madeira et al. (2023), there are several key differences to highlight based on this simulation.Of course, this model has the advantage of modeling the system for a much longer timescales that what is possible with a direct N -body approach.In addition, their model can provide useful insight into the formation of binaries.However, the result presented here (as well as the result of Hyodo & Sugiura (2022)) demonstrates that the mass-shedding and accretion timescales are relatively short, meaning that it may not be necessary to simplify the problem to one dimension.In addition, the insight of the simplified model of Madeira et al. (2023) is only useful to the extent to which it reflects the true nature of the problem.Madeira et al. (2023) suppose an idealized disk in which only a single satellite accretes from the Roche limit at a time.On the contrary, this simulation demonstrates that following a realistic YORP-driven mass-shedding event, multiple moonlets can form simultaneously before undergoing a chaotic series of close encounters and mergers until one satellite remains, which cannot be captured by a 1D model.It is not clear that mass shedding would actually lead to a disk matching the assumptions of Madeira et al. (2023).
Recently, Madeira & Charnoz (2024) extended their study, claiming that their model explains the oblate shape of Dimorphos, the contact binary shape Dinkinesh's recently discovered satellite, Selam, and the prolate shapes of other binary asteroids, where the key difference between these three outcomes solely comes down to the mass-shedding timescale.Some caution should be applied here until the model can address some key issues.First, the mass-shedding timescales explored by Madeira & Charnoz (2024) are chosen arbitrarily and need some physical justification, as there is currently no connection between the initial conditions of their model and YORP processes or the geophysical properties of the primary body.In addition, these 1D disk-moon models do not directly account for the gravitational interactions and subsequent mergers between moonlets.Rather, mergers are assumed to occur anytime two moonlets enter one another's Hill sphere.The shape outcome of these mergers is not demonstrated either; rather a post-merger shape is assumed based on simulations of mergers of Saturn's inner moons (Leleu et al. 2018).However, the geophysics and dynamics of mergers leading to the formation of Saturn's inner moons are very different than in the case of two merging satellites of an S-type asteroid.For example, the moons Pan, Atlas, and Prometheus have densities below 1 g cm −3 and orbit within or near Saturn's Roche limit while Dimorphos is thought to have a bulk density on the order of ∼2.4 g cm −3 and any mergers leading up to its formation are thought to occur beyond the Roche limit, according to Madeira et al. (2023).Any model that can explain the shapes of Dimorphos, Selam, other binary satellites needs to self-consistently address the spin-up, mass-shedding, accretion, and (potentially) merger processes.For example, the shape of a post-merger satellite will depend on parameters such as its density, friction angle, and cohesion.At the same time, these parameters are intimately connected with the primary's failure mechanisms and the mass shedding process, which will then determine the initial conditions of any orbiting debris.This is a challenging problem, and the work presented here only attempts to address a small part of this process.

Simulation set up
The simulation shown in Figs. 1 and 2 are computational expensive requiring ∼105 particles, are limited to short integration times and are impractical for determining statistical outcomes.Therefore, we use the above result, along with known properties of the Didymos system, to inform a simplified initial condition for a large suite of simulations.This allows us to determine the properties of the resulting satellite from a statistical point of view.It is important to note that the primary's shape, density, internal structure, and material properties will all play some role in determining the initial conditions of the shed material (e.g., Sánchez & Scheeres 2020;Walsh et al. 2012;Zhang et al. 2022;McMahon et al. 2020).Therefore, the simplified initial conditions presented here are intentionally generic and may not be perfectly representative of a fully realistic scenario.However, we will show this generic initial condition produces many of the properties of both Dimorphos and other binary systems.
The primary was treated as a uniform oblate spheroid with semiaxes a = b > c, with an equatorial radius of a = b = 400 m and a/c = 1.5, and a bulk density of 2.4 g cm −3 , similar to the best estimate for Didymos's shape and density at the time this study began (Daly et al. 2023b).In order to reduce the total number of particles in the simulation, the primary is modeled as a hollow shell of particles locked together into a rigid aggregate.In order for the primary to have moments of inertia as if it were a solid body (which is important for both gravity and collisions), we include a single point-mass particle at the center, and then adjust the mass of the central particle and the shell particles to achieve moments of inertia of an equivalent uniform oblate spheroid 5 .In other words, the point-mass and surrounding spheroidal shell collectively behave as if they are a single, solid body with the correct moments of inertia.As shown in Fig. 3, the radius of the shell particles are over-inflated to prevent small particles from falling through the cracks.This approach enables a physically realistic solid-body primary without requiring an excessive number of particles in the interior.The spin period of the primary is initially set to 2.5 h in all simulations, approximately matching the rotation period of the primary from Section 2.1 following its rotational disruption.
The orbiting material is generated following a power-law size-frequency distribution (SFD), with a power-law index of -3 and particle radii ranging between 3 and 10 m.This SFD is informed by (but does not match exactly) the observed boulder SFD on Dimorphos (Pajola et al. 2022;Pajola et al. 2023b), as well as the boulder SFD seen on other rubble pile asteroids (e.g., Dellagiustina et al. 2019;Burke et al. 2021;Michikami et al. 2019;Michikami & Hagermann 2021).Including the smallest boulders observed on Dimorphos is impractical due to computational constraints; however, the smallest particles in the simulation are 6 m in diameter, approximately the size of the two largest boulders Bodhran Saxum (∼6.1 m) and Atabaque Saxum (∼6.5 m) near the DART impact site (Daly et al. 2023b).Therefore, we consider this particle resolution sufficient for modeling the formation and shape of a Dimorphos-like satellite, as the particle sizes in these simulations approach the actual sizes of individual large boulders on Dimorphos's surface.Cohesion between boulders is ignored in these simulations, as cohesion forces on rubble piles likely arise from Van der Waals forces between fine-grained material less than ≲10 µm (e.g., Sánchez & Scheeres 2014).In a mass shedding scenario, much of this fine-grained material (if present) would be released when the mass is first shed, and when moonlets undergo collisions and tidal disruptions.For a Dimorphos-like system, solar radiation pressure would rapidly remove these fine grains (e.g., Ferrari et al. 2022).Therefore, we might expect the cohesion of the secondary to be no higher than the primary's cohesion (if present).This may explain why Didymos's surface requires a small amount of cohesion to maintain its surface stability while Dimorphos does not (Barnouin et al. 2023).
A disk of orbiting debris is placed in a circular region extending out to 1.5 times the equatorial radius of the primary, which approximately corresponds to the effective Roche limit for a body with a friction angle of 35 • (Holsapple & Michel 2006).The disk also has a finite thickness of 40 m (roughly two diameters of the largest disk particle), to allow enough space to initialize the required number of particles.The example simulation presented in Section 1.2 demonstrates that mass shedding does not lead to a stable disk but instead material clumps almost immediately to form moonlets.To approximate this effect here, particles are purposefully not given any initial velocity dispersion to ensure a disk that is gravitationally unstable (i.e., a Toomre Q < 1) so material will begin clumping immediately.All disk particles have a density of 3.5 g cm −3 , to match the grain density of L and LL chondrites (Flynn et al. 2018), which are the best meteoritic analogues for Didymos (de León et al. 2006;Dunn et al. 2013).
Each disk particle is placed on a circular orbit, accounting for the mass of the primary and its oblateness, as well as the mass of all other enclosed disk particles.Each disk particle's mean motion (i.e., its initial orbital angular velocity), can be written as a function of its radial distance r i , where M A , R eq A , and J 2 are the primary's respective mass, equatorial radius, and oblateness, M (< r i ) is the mass of all the disk particles enclosed within particle i's position, and G is the gravitational constant.Owing to the finite thickness of the disk, particles are initialized with inclinations on the order of 2 − 3 • .Particles are generated in a symmetric disk to simplify the process of generating initial conditions, which is not a realistic starting condition.However, we emphasize that, owing to the disk being gravitationally unstable, collisions and close encounters quickly excite the eccentricities and inclinations of the disk particles, which we deomnstrate below rapidly leads to particle orbits representative of the more realistic starting conditions following a mass-shedding event, such as that shown in Fig. 1.
Any particles that reach a distance of 40 km (100R eq A ) are automatically deleted from the simulation.This boundary was set purely as a precaution to prevent a single particle from being ejected from the system, which could significantly slow down pkdgrav's tree due to the single extremely distant particle, and corresponds to ∼2/3 of the system's Hill sphere if it were located at 1 AU.We found that a small fraction of particles end up reaching this boundary and are deleted.Typically only 1-2% of the disk's initial mass is ejected and we verified that the vast majority of the particles that hit this boundary and were deleted were either on escape trajectories (e orb > 1), or had an apoapse distance well outside the Hill sphere and would have not returned to the binary system.Therefore, deleting these particles has a negligible effect on the binary's formation.In a some cases, a large aggregate is ejected from the system forming an asteroid pair and is discussed in Section 3.7.1.
A core motivation of this study was to demonstrate that a Dimorphos-sized satellite could plausibly form via a single mass shedding event.Therefore, most simulations simply vary the initial mass in the disk (i.e., number of particles), however we also vary the friction angle of the material between 29 • and 40 • .For each set of parameters, we randomize the initial locations of the particles to understand the chaotic nature of the satellites formation.For each disk mass (M disk ) and friction angle (ϕ), we run a given number of 'clone' simulations, which are listed in Fig. 1.In total, we run 120 simulations.96 of the simulations have ϕ = 35 • and have an initial disk mass ranging between 0.02M A and 0.04M A .These 96 simulations constitute the bulk of the results shown in the following sections.With the remaining 24 simulations, M disk is kept fixed at 0.03M A and friction angles of 29 • , 32 • , and 40 • are tested.All simulations were limited to ∼100 d (5.8 × 10 7 steps) due to the high computational cost imposed by the small timestep.Each simulation is assigned a number between 1 and 120, and the result of each simulation is tabulated in Table 3 in the Appendix.Note-In all simulations, the primary consists of 2,371 particles, while the number of debris particles varies between ∼4, 200 and ∼8, 400, depending on the total mass in the disk.See Table 3 to see a full listing of each simulation and its result.3. RESULTS

An example case
First, we show a representative simulation to demonstrate the model.In Fig. 4 we show 8 snapshots of a simulation in which the M disk = 0.02M A and ϕ = 35 • .Almost immediately, particles begin clumping and forming short-lived spiral arms due to Keplerian shear, as a result of the disk being gravitationally unstable (e.g.Kokubo et al. 2000).A couple days later, several moonlets form and undergo a chaotic series of close encounters and mergers until a single large satellite remains.Interestingly, the satellite is born in a tidally-locked configuration with the primary.The immediate tidal locking of the secondary is likely due to several factors, including its low eccentricity and formation near the Roche limit.In addition, all the mergers in this simulation occur rather gently without significantly perturbing the secondary's spin state.As we will see in the next several subsections, the chaotic nature of the satellite's formation leads to a wide range out outcomes in terms of the satellite's physical and dynamical properties.Therefore, this simulation is not necessarily typical of all cases.Starting from the top, we plot the satellites DEEVE axis ratios, mass ratio, orbital distance, and Euler angles, which describe its orientation in a frame rotating with the orbit.In only a handful of days, the satellite reaches a near-circular orbit at ∼3RA with a mass just under 0.01MA.The rotation state of the satellite is excited, but tidally locked, given by all three Euler angles being large but still well under 90 • .
In order to compare whether the simplified initial condition of the disk is a reasonable representation of the more realistic initial conditions following mass shedding, as demonstrated in Section 2.1, we compare the distributions of semimajor axis, eccentricity, and inclinations of orbiting debris at early times.Fig. 6 compares the orbital elements of orbiting material between the more realistic mass-shedding example from Fig. 1 and the simplified case shown in Fig. 4. Owing to the gravitationally unstable disk and the role of collisions, the orbits of disk particles are rapidly excited to wider orbits, with higher eccentricities and inclinations, despite initially starting with circular, nearly co-planar orbits.At least qualitatively, the distribution of orbital elements several days into the simulation reasonably reflects the orbits of post-mass-shedding debris (i.e., Fig. 1), although the more realistic initial conditions have a slightly wider distribution in both inclination and eccentricity.1 and the simplified initial conditions from the simulation shown in Fig. 4. The orbital elements are plotted when moonlets first begin forming.For the mass-shedding simulation, this corresponds to ∼0.5 d after mass is initially shed, while for the simplified disk, this occurs around ∼4 d into the simulation.Although the distribution of orbital elements between the two cases does not match perfectly, due to the gravitationally unstable initial conditions of the simplified disk, we arrive at a similar range in semimajor axis, eccentricities, and inclinations seen in the more realistic simulation.

Satellite mass and density
In Fig. 7, we plot the secondary-to-primary mass ratio for all simulations with ϕ = 35 • in order to understand how the mass of the initial disk determines the resulting satellite mass.Assuming the primary and secondary have the same bulk density (which is approximately true in these simulations), we also provide the size ratio DB DA ≈ MB MA 1/3 on the second y-axis.In the 96 simulations shown here, the satellite tends to reach its final mass in the first few tens of days apart from a select few special cases that undergo late mergers or disruptions.For context, the Dimorphosto-Didymos size ratio is ∼0.2 which corresponds to a mass ratio ∼0.01 if the bodies have equal densities (Daly et al. 2023b).Therefore, we find that an initial disk mass of only ∼0.02 − 0.03M A is capable of producing a Dimorphos-mass satellite in only a matter of days.Although their study focuses on a regime where significantly more mass and angular momentum is put into orbit, we find that our formation timescale is broadly consistent with that found by (Hyodo & Sugiura 2022).Our result is orders of magnitude lower (both in time and mass) than the calculation of Madeira et al. (2023) which requires 25% of Didymos's mass to be shed into a ring that will then take years to form a Dimorphos-mass satellite.This disagreement likely stems from the different approach to the problem, namely in the initial conditions.Our model starts with a disk of particles initialized on much wider orbits, given the existing body of literature on spin-up and mass shedding (e.g., Yu et al. 2019;Zhang et al. 2018;Zhang et al. 2022;Hyodo & Sugiura 2022) and based on the spin-up example provided in Section 2.1, whereas the Madeira et al. (2023) model supposes that the orbiting debris starts in a narrowly confined region at the surface of the primary and slowly spreads outwards, which substantially increases the timescale for the satellite's formation.
We plot a histogram of the bulk density of each satellite at the end of the simulation in Fig. 8.The volume of the rubble pile is computed by the concave hull (or "alpha shape"; (Edelsbrunner 1995)) of the set of points defined by the edges of the outermost spheres which make up the satellite.This method of calculating bulk density provides high accuracy by providing a "tighter fit" to the true shape of the rubble pile than other volume estimates such as the volume of the convex hull or the dynamically equivalent ellipsoid.Since some of the large void spaces could be filled with smaller boulders that are below the resolution of the simulations, the measured bulk density here is only notional and should be thought of as a soft lower limit.0.0 0.17 Figure 7.The secondary-to-primary mass ratio over time, for all simulations with ϕ = 35 • .The accretion of the satellite is highly efficient, occurring in only a matter of days.The second y-axis shows the secondary-to-primary size ratio, assuming that the two bodies have the same bulk density (which is approximately true in these simulations).The spread in final mass among disks with the same initial mass is due to the chaotic formation history of each system.Note: most discontinuities in this plot are real and are due to mergers or tidal disruptions of the satellites.However, there are a couple discontinuities that result from the nearest-neighbor search algorithm misidentifying the largest orbiting fragment.
2.12 2.14 2.16 2.18 2.20 2.22 2.24 ρ bulk [g cm The density of all individual spheres is 3.5 g cm −3 , so these rubble piles have packing fractions of ∼0.62.Realistically, it is possible for the "real" bulk density to be slightly higher than what is measured here, as small particles which are below the resolution limit of these simulations would fill in some void space.

Satellite orbit and rotational state
The majority of the simulations end with a single satellite on a near-circular orbit having a semimajor axis between ∼2.5 and ∼4 primary radii.Fig. 9 shows the final semimajor axis and eccentricity for the satellite in the 96 simulations where ϕ = 35 • .Generally, we find that more massive disks tend to produce a satellite on a wider, more eccentric orbit, simply as a result of the disk having a larger initial angular momentum.This is demonstrated in Fig. 9 where we show the mean a orb and e orb along with their standard deviations for each of the 3 different disk mass cases.
For each satellite, we determine its rotation state based on the 1-2-3 Euler angle set, consisting of the satellites roll, pitch, and yaw angles in a frame rotating with its orbit (see Agrusa et al. (2021)).When the roll and pitch angles are small, the yaw angle can be thought of as the classic planar libration angle, i.e., the angle between the secondary's long axis and the line-of-centers.However, many of the satellites have undergone several mergers or close gravitational encounters, and are in excited, non-planar rotation states, necessitating the use of the Euler angle convention.In Fig. 10, we plot the maximum roll and yaw angle for each satellite over the final 10 days of the simulation.We see three distinct regions in this plot.If the yaw angle stays below ≲60 • , we consider the satellite to be in synchronous rotation, since its long axis stays pointed in the direction of the primary.Out of the synchronous rotators, a minority are in "pure" (albeit highly excited) synchronous rotation where the roll and pitch angles are also ≲60 • .Most of the synchronous rotators, however, are in a rolling state about their long axis.In this so-called "barrel instability", the satellite remains on-average tidally locked to the primary although it continues to roll about its long axis ( Ćuk et al. 2021).This non-principal axis rotation state within the 1:1 spin-orbit resonance could be long-lived as this mode would dissipate inefficiently by tides.Since binary YORP (BYORP) requires a synchronous secondary, this would also substantially weaken the BYORP effect or prevent it entirely ( Ćuk & Burns 2005;Quillen et al. 2022).Finally, we also see many satellites in an end-over-end tumbling state, where their roll and yaw angles both reach 90 • .Generally, the tumblers have a higher eccentricity than the synchronous rotators as indicated in the plot, however, the onset of chaotic tumbling also depends on the satellite's inertia ratios (Wisdom et al. 1984;Agrusa et al. 2021) as well as its collision history.
Of course, it is no surprise that these satellites are never perfectly tidally-locked, as these are short-term simulations in which the satellite has undergone many collisions and mergers so its rotation state is naturally excited.However, synchronization is the fastest-evolving tidal effect (Goldreich & Sari 2009), and we would therefore expect the satellite's free libration to damp on relatively short timescales, potentially within hundreds of years since they are already in synchronous rotation (Meyer et al. 2023).This may provide a simple explanation for why we observe so many synchronous satellites (Pravec et al. 2016), without needing to invoke the more complicated dynamical processes of rotational fission to achieve this equilibrium (Jacobson & Scheeres 2011).Instead, if the satellite forms via accumulation of shed material, then it has a reasonable chance to form in or near a tidally-locked state.
We are not aware of any studies that demonstrate a binary asteroid can immediately form in a synchronous rotation state following rotational disruption, although this seems like a relatively natural outcome.Formation with synchronous rotation has been found in other studies on gravitational accumulation near the Roche limit, such as the accretion of moonlets from Saturn's rings (Karjalainen & Salo 2004) and from circumplanetary disks (Hyodo et al. 2015).Although the variance is quite large, a larger disk mass typically leads to a higher eccentricity and larger semimajor axis, due to the disk having a higher angular momentum and more collisions that can drive up the largest satellite's eccentricity. .The maximum roll and yaw angles for each satellite, based on the final 10 days of the simulation.There are three distinct post-formation rotation states for the satellites, which are annotated on the plot to guide the reader's eye.We consider the satellite to be in synchronous rotation if its long axis remains aligned towards the secondary (i.e., yaw ≲ 60 • ).The "barrel instability" is a subset of the synchronous rotators where the secondary's long-axis stays largely aligned with the primary, yet it is able to roll about its long axis like a barrel.Finally, a significant fraction of satellites form in a tumbling state where all three Euler angles (roll, pitch, and yaw) can reach 90 • .

Satellite shape
While there is significant literature on binary asteroid formation, there are no studies to our knowledge that directly model the expected shapes of the satellite.Many studies consider the shape of the secondary, but do not model that shape being formed directly (e.g.Jacobson & Scheeres 2011;Davis & Scheeres 2020).
It is important to be clear with the definition of the satellite's shape.In this study, we define the shape of the satellite by its three principal axis a, b and c, which correspond to the body's three principal moments A, B, and C. We measure its axis lengths in two different ways.The first is simply the physical extents of the body along these axes.The second measure is the axis lengths of the dynamically equivalent equal-volume ellipsoid (DEEVE).This is a uniform density ellipsoid having the same mass and moments of inertia of the rubble pile.If the body has an approximately ellipsoidal shape, then these two measures of its shape will match closely.Measuring the body's axis ratios by its DEEVE can be useful as they don't fluctuate significantly as a result of the motion of a single particle on the surface, which is common as the satellite is forming.Therefore, we use the DEEVE axis ratios in any time series plots as they are much less noisy.However, the physical extent of the body tends to better represent the "true" shape of the body and is most analogous to real life observations.In our simulations, these two measures tend to differ significantly for highly irregular shapes (high a/b and/or b/c).This is demonstrated in Fig. 11, where the DEEVE axis ratio tends to be quite larger than the physical extent axis ratio for large axis ratios.Figure 11.The DEEVE axis ratios of each simulated satellite compared to the body's axis ratios defined by its physical extents.Points lying under the diagonal line have a DEEVE axis ratio (either a/b or b/c) that is greater than its physical extents and vice versa.Points lying on the line indicate that the DEEVE shape is a good approximation of the physical shape.For reference, we include the axis ratios of Dimorphos defined by its physical extents and best-fit ellipsoid, demonstrating that Dimorphos's shape is exceptionally well-fit to an ellipsoid.Dimorphos's a/b and b/c axis lengths, including 1σ uncertainties, are indicated in blue and black, respectively (Daly et al. 2023a) Due to the discrepancy between the DEEVE-derived and extend-derived shapes highlighted above, we compare both quantities to known secondary shapes for thoroughness.In Fig. 12 we show both the DEEVE-and extent-derived semi-axis ratios a/b and b/c at the end of each simulation with ϕ = 35 • .We also include the shapes of the satellites of 66391 Moshup (formerly 1999 KW 4 ), 2000 DP 107 , and 2001 SN 263 , which are the only publicly available radar-derived shapes for the satellites of other binary asteroids (Ostro et al. 2006;Naidu et al. 2015;Becker et al. 2015), as well as the axis ratios of Dimorphos (Daly et al. 2023b,a).The updated shape of Dimorphos provided by (Daly et al. 2023a) differs slightly from the initial assessment in (Daly et al. 2023b), but the difference is small enough that we only plot the latest values to avoid confusion.For each real asteroid system, we include 1σ uncertainties in a/b and b/c, assuming that the reported uncertainties in a, b and c from each respective paper are uncorrelated, which may not be true.Therefore, the uncertainties should only be used to guide the reader's eye.The satellites of Moshup and DP 107 are the best comparisons with Dimorphos, as they are both S-type binaries, whereas SN 263 is a C-type triple (with large uncertainties in the satellite shapes).Fig. 12 demonstrates that, generally speaking, these simulations produce satellites that are more elongated (high a/b) and more flattened (high b/c or a/c) than the radar-derived secondaries.However, there are many cases that produce shapes similar to radar-observed secondaries, given their uncertainties.Although no satellites are produced within the 1σ uncertainty region of Dimorphos's shape, several simulations do come close.Dimorphos's b/c ratio lies approximately in the middle of the simulated b/c range.Although there is a preference for elongated satellites, several simulations result in a low elongation (a/b ≲ 1.1), like Dimorphos.These simulations demonstrate that more elongated shapes are preferred, although immediately forming a Dimorphos-like shape by mass shedding is not implausible.Due to the satellite's accretion near the Roche limit, this strong preference for more elongated shapes comes as no surprise and has been seen in analogous studies (e.g.Porco et al. 2007;Hyodo et al. 2015).It is important to note that the simulations here are run for only 100 days and there could be longer term processes that will modify the satellite's shapes (impacts, landslides, etc.), so it is important to interpret any comparison between the simulated and observed shapes with caution.
We also compare our results to the lightcurve-derived shapes from Pravec et al. (2016Pravec et al. ( , 2019)).Fig. 13 shows a histogram of the simulated a/b ratio compared to the dataset of Pravec et al. (2019) along with some new (not yet published) secondaries (Pravec, personal communication).Here, we only include lightcurve secondaries in which the primary is less than 20 km in diameter and the secondary-to-primary size ratio is less than 0.6, to ensure that we are only comparing with binaries likely formed by YORP.
We find a remarkable agreement between the shapes of lightcurve secondaries and those simulated in this work.Both data sets show that most satellites have a/b < 1.6, aside from a few outliers.The sharp drop-off in satellites high a/b has been attributed to chaotic dynamics of satellites with a/b ⪆ √ 2 that ultimately destroys (or alters) them ( Ćuk & Nesvorný 2010;Pravec et al. 2016).While this is a real dynamical effect, we suggest that it may also be that secondaries simply don't form with extremely high elongations to begin with.
Interestingly, there are several secondaries with low elongations (a/b ≲ 1.1) measured by Pravec.Naively, it appears as if the simulations do a good job at matching this low elongation population.However, it is important to note that lightcurves are heavily biased against detecting satellites with a a/b ≲ 1.05 (Pravec, personal communication).This fact, combined with Dimorphos apparently having a/b < 1.1 suggests the existence of a significant population of secondaries with low a/b.If the simulations presented here reasonably capture the formation of these secondaries, we would expect the simulations to produce an overabundance of low elongation secondaries compared to the observed population.This suggests that either (1) there are other longer-term processes at play that could reshape some satellites over time or (2), there is an effect not captured by our model resulting from our assumptions or initial conditions.
Figure 13.Lightcurve-derived shapes are only sensitive to a/b, so we construction histograms of a/b from the simulations (defined by their physical extents) to compare with a/b of satellites derived from lightcurves (Pravec et al. 2016(Pravec et al. , 2019)).The y axis is not normalized since we are comparing a similar number of real systems (87) and simulations (96).

Effect of friction angle
We conducted an additional set of simulations where the disk mass was held constant at M disk = 0.03M A and the friction angle was varied between 29 • and 40 • .For each friction angle (29 • , 32 • , and 40 • ), we conducted 8 simulations with randomized initial conditions.With this smaller sample size, we are dealing with small number statistics but can still draw a few basic conclusions.
As the friction angle is decreased, there is significantly less variation between simulations in general, including the satellite's formation distance, eccentricity, shape, and rotation state.In Fig. 14a we show the satellite's final semimajor axis and eccentricity, along with an average and standard deviation for each value of ϕ.Although this is small number statistics, there is a clear trend with friction; satellites with a lower friction angle tend to form on a closer, more circular orbit.This is because the violent processes such as collisions and close gravitational encounters that lead to highly excited orbits tend to disrupt lower friction moonlets.This means that the moonlets, which eventually merge to form the final satellite, are necessarily on more circular orbits.Therefore, the final satellite will tend to have a less excited orbit.
The axis ratios (based on the physical extents) are given in Fig. 14b.In general, there is a smaller variation in the secondary's shape as the friction is reduced.In addition, a lower friction tends to result in a less flattened object (i.e., more prolate).As the friction angle is lowered, the body is forced to reshape itself until it is structurally stable, and bodies with higher frictions are able to maintain a wider range of shapes.In the limit that the friction angle were reduced to 0 • , the satellite effectively behaves like a fluid and would take on the shape of a Roche ellipsoid, at its given orbital distance (Holsapple & Michel 2006).This is not demonstrated here, simply because a simulation with ϕ = 0 • would take an excessive amount of time to form a satellite.
A histogram of the secondary's bulk density as a function of friction angle is shown in Fig. 14c, where there is a clear trend.Bodies with a lower effective friction are able to achieve a more efficient packing arrangement, thus having less void space and a higher bulk density.At a grain density of 3.5 g cm −3 , this figure should be thought of as a lower limit for the secondary's bulk density.In reality, some void space could be filled with boulders that are below the resolution of these simulations, so the "true" bulk density might be be higher.The axis ratios of the satellite, based on its physical extents, as a function of friction angle.(c) The bulk density of each satellite as a function of its friction angle ϕ.Generally, bodies with a lower effective friction are able to achieve a more efficient packing and thus have a higher bulk density.

Matching the Shape of Dimorphos
Although these simulations demonstrate that immediately forming an oblate satellite is relatively rare, we show two cases with a ϕ = 35 • that achieve shapes similar to Dimorphos.Although neither cases are a perfect match, they do suggest that an oblate spheroid-like shape can plausibly be formed from a single mass shedding event.In Fig. 15, we show the two examples, both of which have an initial disk mass of M disk = 0.04M A .In the first example (Fig. 15a), the satellite grows rapidly over the first several days until it undergoes a close encounter with another satellite (not shown), which causes the satellite to move inwards and undergo a partial tidal disruption event with the primary.This torques the remaining mass onto a much wider, eccentric orbit and also causes the satellite to rotate asynchronously yet maintain principal axis rotation (as demonstrated by the roll and pitch angles remaining small and the yaw angle circulating).As a result of the asynchronous rotation, the satellite begins accreting material in an azimuthally symmetric manner, as it has a different orientation at each periapse passage (where there is material available to accrete).The satellite also undergoes a merger with the last remaining moonlet at ∼20 d, which abruptly This process suggests that an oblate shape like that of Dimorphos could plausibly be acquired if the satellite accretes material after some dynamical process triggers an asynchronous rotation state.
In the second example (Fig. 15b), the satellite achieves its oblate shape simply by undergoing two near-head-on mergers with other moonlets at early times (around ∼3.5 d and ∼7 d, respectively).These mergers serendipitously lead to a more oblate shape due to their favorable geometry.In addition, the satellite's relatively circular orbit, which is well outside the Roche limit, prevents its shape from undergoing significant further modification, despite its tumbling rotation state.Mergers among similar-sized moonlets have been proposed to explain the oblate shapes of some of Saturn's small moons (Leleu et al. 2018).
In Fig. 16, we show renderings of these two examples in each body's respective body-fixed frame, using the same lighting conditions and viewing geometry as the DART spacecraft on its approach to Dimorphos.Although these two examples demonstrate that an oblate satellite can plausibly form by asynchronous rotation during accretion or by lucky mergers, these events are relatively rare.Out of over 100 simulations, Dimorphos-like shapes are only formed a few percent of the time.Therefore, it is feasible that Dimorphos obtained its shape through some other unexplored mechanism.

Unique outcomes
Due to the chaotic orbital and collisional evolution of the moonlets, some simulations resulted in outcomes that were not expected at the conception of this study.This includes the creation of asteroid pairs and triple systems.To be clear, we do not claim that all or even a majority of pairs or triples originate from a single mass shedding event, rather that this formation mechanism is plausible for some systems.

Asteroid Pairs
An asteroid pair is where two unbound asteroids have similar enough orbital elements (and colors/spectra; see Pravec et al. (2019), and references therein) such that they can be traced back to a common origin.An asteroid cluster is where multiple unbound bodies can be traced back to the same system.Many pairs and clusters are thought to form either through collisions, or through rotational fission (Vokrouhlický & Nesvorný 2008;Pravec et al. 2010Pravec et al. , 2018)).Because the dynamics of a newly formed binary system are chaotic and have a positive free energy for mass ratios below ∼0.2, the secondary can be ejected purely by the internal dynamics of the system (Scheeres 2007b;Jacobson & Scheeres 2011).This phenomenon does an excellent job at reproducing the observed correlation between the primary's rotation period and the mass ratio of the asteroid pair (Pravec et al. 2010(Pravec et al. , 2019)).There is a distinct correlation between asteroid pair mass ratio and the primary's spin rate; pairs with a large mass ratio tend to have slower primary rotation periods, presumably because more angular momentum was transferred from the rapidly rotating primary to the secondary in order eject it from the system.
Due to the azimuthal symmetry of the primary in these simulations, there is very weak coupling between the primary's rotation and the satellites orbit, preventing these simulations from exploring this correlation.However, our simulations demonstrate that pairs/clusters can also form via three-body interactions between the moonlets during their formation.In all cases where we form a successful pair or cluster, there is always a secondary still bound to the primary.In other words, in any case where we form a pair, it is always a "paired binary".
As an example, we show a time-series plot of the semimajor axis and eccentricity of the two largest satellites from Disk 040 (M disk = 0.03M A , ϕ = 35 • ) in Fig. 17.In this case, the smaller satellite (M 2 ) has a series of close encounters with the larger satellite (M 1 ), which increase its semimajor axis and eccentricity until its eccentricity eventually exceeds 1, placing it on an escape trajectory.The particle then reaches the simulation boundary and is removed from the simulation.In Fig. 18, we show snapshots of this simulation at the moment when the M 2 is ejected.The two satellites have a close encounter at 48.3 d, where M 1 torques M 2 onto a hyperbolic orbit, ejecting it from the system.Meanwhile, as a consequence of angular momentum conservation, M 1 's eccentricity is raised (i.e., its periapse distance decreases) causing it to have a close approach with the primary at 48.9 d.M 1 then undergoes a partial tidal disruption, loosing some mass, reshaping, and torquing it back onto a higher orbit safe from further tidal disruption.This process also breaks M 1 from its previous synchronous rotation state.Figure 17.An example of a paired binary being formed as the larger satellite ejects the smaller satellite from the system.After an extremely close encounter, the smaller satellite is ejected from the system (e orb > 1) before it hits the simulation boundary at 100R eq A where it is then removed.This particular simulation is Disk 040 (M disk /MA = 0.03, ϕ = 35 • ) and a movie of this simulation is provided.
Figure 18.Simulation snapshots corresponding to Fig. 17 of the final close encounter that ejects the smaller satellite from the system.When this occurs, larger satellite's periapse distance is lowered and its eccentricity is raised.This leads to a partial tidal disruption where the remaining satellite loses some mass and breaks from synchronous rotation.A movie of this simulation is provided.
In Fig. 19a, we plot the size of each ejected fragment in relation to both the primary and the secondary, from all simulations where ϕ = 35 • .It is important to note that this figure excludes fragments consisting of one or two particles, of which numerous are ejected from every simulation.Here, we only show ejected fragments that contain 3 or more particles.Also, this figure only includes fragments that have been ejected within the ∼100 d pkdgrav integrations.However, as we will see in the next section, there are several systems which still contain two satellites after 100 d that will be ejected at a later time and would also count as a successful pair.In general, we find that a larger initial mass shedding event tends to eject larger fragments, as there is more material and angular momentum available.For context, the smallest size ratio observed for an asteroid pair is roughly ∼0.1, so this is a size range that would be potentially discoverable (Pravec et al. 2019).We also show the spin rate of each ejected fragment (Fig. 19b), which demonstrates that most fragments are spinning relatively fast, due to the high angular momentum nature of ejecting the fragment.The observed population of asteroid pairs shows a much broader distribution of spin rates of ejected secondaries, ranging between ∼1 − 9 d −1 (Pravec et al. 2019).In addition, many of the pairs formed in these simulations are in non-principal axis (NPA) rotation, whereas all observed secondaries of asteroid pairs for which there is sufficient lightcurve coverage appear to be in principal axis rotation.This suggests that these fragments undergo significant spin evolution after being ejected (through damping and/or YORP) or other mechanisms besides three body interactions are responsible for the production of most asteroid pairs.The axis ratios are also shown in Fig. 19c and are consistent with the axis ratios of observed pairs inferred from their lightcurve amplitudes (Pravec et al. 2019).Finally, we also show the eccentricity of the fragment before it is deleted from the simulation.The vast majority are on hyperbolic (escape) trajectories, having e orb > 1.However, some of these fragments are still bound to the system.If these fragments were kept in the simulation, they would likely have been ejected or collided during a future close encounter.If these are near-Earth asteroids, these satellites would be lost anyways, having an apoapse outside of the system's Hill sphere.However, if these are main-belt asteroids with larger Hill spheres, there is a chance that these objects deemed "pairs" could remain bound with wide eccentric orbits.
For example, asteroid (3749) Balam is a main-belt triple with a fast-rotating (2.8 h spin period), ∼4 km primary with a relatively close inner satellite, separated by ∼13 km on a near-circular orbit (Pravec et al. 2019).The outer satellite is on a wide orbit with an eccentricity between 0.35 and 0.77, (Merline et al. 2002;Vachier et al. 2012).The inner and outer satellites have respective satellite-to-primary size ratios of ∼0.46 and ∼0.24 (Pravec et al. 2019), making them relatively large.In addition, Balam is a member of an asteroid pair (Vokrouhlický 2009), with backwards numerical integrations indicating that asteroid 2009 BR 60 separated from Balam less than one million years ago (Vokrouhlický 2009;Pravec et al. 2019), making the Balam system very young.Due to its high eccentricity, it has been suggested that Balam's outer satellite could have been formed as the result of a collision (Durda et al. 2004) or by capture (Marchis et al. 2008), while the inner satellite is consistent with YORP spin-up and rotational disruption due to the primary's fast rotation (Polishook et al. 2011).Since the simulations in this study considered much more mild mass shedding events and removed these eccentric satellites, there are no direct analogs to the Balam system in our study.However, our results suggest that a plausible explanation for the Balam system could be formation via a single, much more massive rotational disruption event.In such a scenario, an inner satellite, an outer eccentric satellite, and an unbound satellite could all be produced at the same time, although the details of such a scenario require further study.

Triples
In some cases, we find that 3-body interactions lead to the formation of systems with two satellites on non-crossing orbits.However, these simulations are run for less than 100 days, making it impossible to determine their stability with pkdgrav alone.In rare cases where two satellites are still remaining at the end of the pkdgrav simulation, we perform a brief test of the system's stability.We hand off the masses, positions, and velocities of the primary and its satellites to the REBOUND N -body code (Rein & Liu 2012;Rein & Spiegel 2015;Rein & Tamayo 2017).These simulations are relatively simple; the two satellites are treated as spheres but we include the effect of the primary's J 2 using the REBOUNDx package (Tamayo et al. 2020).As these simulations do not include higher-order effects such as spin-orbit coupling of the two satellites, tidal dissipation, BYORP, etc., we only integrate the system for 10 3 years and then check if both satellites still remain.In reality, higher-order effects may strongly affect the system's evolution and stability.Additionally, the primary is largely azimuthally symmetric, and a more realistic primary with a nonnegligible C 22 may make these close-in triple systems less stable.Here, we merely demonstrate that the formation of a triple system from a single rotational disruption event is plausible.In total, we find 11 triple systems after the 100 d pkdgrav simulations, of which 7 survive after further integration with REBOUND.
We find that many of the stable triples achieved their stability due to capture into mean motion resonances (MMRs) during the pkdgrav stage.While moonlets are growing, their orbits are constantly perturbed by collisions, mergers, and close encounters, so capture into a MMR is effectively random.Out of the 7 stable triples, 5 are in MMRs, while two are not in resonance.
In three of the resonant cases, the two satellites form in a co-orbital configuration (i.e., the 1:1 MMR) by co-accretion.When this occurs, the smaller satellite occupies the L 4 or L 5 Lagrange point of the larger body.Snapshots of the three cases are shown in Fig. 20.These systems satisfy the condition for stability, in which M1+M2 M +M1+M2 ≲ 0.04, where M is the primary's mass and M 1 and M 2 are the two satellite masses (Murray & Dermott 2000;Laughlin & Chambers 2002;Figure 20.Three cases where the co-orbital configuration of two satellites remains stable after a 10 3 yr integration with REBOUND.The top panels show a view of the pkdgrav simulation after ∼100 d prior to handing off to REBOUND, looking down from the primary's spin pole.We encourage the reader to view the provided movies of these simulations.The corresponding bottom panels show the position of the two satellites in the rotating frame of the larger satellite at each output over the 10 3 yr REBOUND integration, demonstrating that the smaller satellite resides in a Lagrange point (either L4 or L5).In the left case (Disk 003), the smaller satellite librates around the trailing Lagrange point (L5) with an amplitude of ∼7.2 • , while in the middle and right cases (Disk 010 and Disk 108), the smaller satellite librates around the leading Lagrange point (L4) with amplitudes of ∼11.4 • and ∼18.7 • , respectively.Movies of these threes simulations are provided.Sicardy 2010).We also numerically test the system's stability using REBOUND and find that the smaller satellite remains in stable libration about the Lagrange point for at least 10 3 yr.We restrict the REBOUND simulations to only 10 3 yr, because we neglect tides and non-gravitational forces which can tighten or loosen the tadpole orbits, depending on the direction of migration (e.g., Fleming & Hamilton 2000).We find that co-orbital systems are only formed when the initial disk mass is small (M disk /M A ≲ 0.03), in agreement with the study of Hyodo et al. (2015) on the formation of multiple satellites with N -body simulations in circumplanetary disk.
We note that the idealized initial conditions of these simulations may increase the likelihood of forming co-orbital satellites and that the primary's azimuthal symmetry likely contributes to their long term stability.Here, we merely point out that forming a co-orbital system via a single mass shedding event is possible, albeit unlikely.A future detection of a triple asteroid with co-orbital satellites, however unlikely, could be explained by co-accretion following a single mass shedding event.
The remaining two resonant cases are capture into the 2:1 and 5:2 MMRs.This occurs after a satellite is first torqued onto a wider orbit (≳4R A ), while there is still material near the Roche limit available to accrete.This allows a second body to grow closer in, which can capture into a nearby MMR with the outer satellite while it is accreting.Once the initial capture occurs, the resonance can further stabilize due to the smooth growth of the inner satellite.
In Fig. 21 we show an example of capture into the 2:1 MMR, where the inner satellite is the most massive.In Fig. 21a, we show a top-down view of the triple system at the end of the 100 d pkdgrav integration along with a time-series plot in Fig. 21b.Early on, both M 1 and M 2 grow rapidly until they temporarily merge into a single body which is then quickly tidally disrupted.The tidal disruption torques one large fragment (now M 2 ) onto a wide orbit.Then, M 2 continues to grow and migrate outwards through a series of mergers and close encounters with other moonlets.Meanwhile a new inner satellite forms (now M 1 ) from one of the disrupted remnants and migrates inwards.This process places M 1 and M 2 close enough to the 2:1 MMR that they capture.The resonance is then further stabilized by the smooth growth of M 1 , the inner satellite.We can confirm the existence of the resonance due to the libration of the resonance argument ϕ res = 2λ 2 − λ 1 − ϖ 2 , where λ i are the respective body's mean longitudes and ϖ 2 is M 2 's longitude of pericenter.We attribute the drift in the resonance argument to the growth in the satellite masses and collisions, both of which make the process non-adiabatic.We confirmed that the resonance argument continues to librate when handed off to REBOUND.We also find one example of capture into the 5:2 MMR, shown in Fig. 22. Capture into such a high-order resonance came as a surprise and may have been enabled by the high eccentricity of the outer satellite.If a triple asteroid were to be observed in such a resonance, then this process should be studied further.
We find two other long-term stable triples (disk 075, 085) with a more massive outer satellite having period ratios close to 2:1 and 4:1 respectively, although we confirm that are not in resonance.However, these systems have wellseparated satellites and thus remain long-term stable in the REBOUND simulations.
Comparing our simulated triples to the known population is a problem of small number statistics.There are very few known triple asteroids, none of which are known to be in mean motion resonances.Out of the near-Earth asteroids, there are three confirmed triples: 1994 CC, 2001 SN 263 , and 3122 Florence (Brozovic et al. 2009;Nolan et al. 2008).Both CC and SN 263 are well-characterized, not near MMRs, and also have some mutual inclination, which as been attributed a previous close planetary encounter (Brozovic et al. 2009;Fang & Margot 2012;Becker et al. 2015).However, no detailed dynamical studies have been published for 3122 Florence.Preliminary estimates have put the orbit periods of the two satellites at 22 h and 7 h (Brozovic et al. 2018), which is not far from the 3:1 MMR, although a more comprehensive analysis of the satellite orbits would be needed to confirm this.Despite our simulations occasionally resulting in a stable triple in a mean-motion resonance, a single mass shedding event is certainly not the only way to achieve such a configuration.For example, a second mass-shedding event occurring after the first satellite has migrated outwards, followed by convergent migration (due to tides and/or BYORP) is an equally plausible mechanism to form a MMR triple.

Simulation outcomes
In Fig. 23, we show the simulation outcomes for all simulations with ϕ = 35 • , summarizing the resulting rotation state of the secondary and whether the system formed a pair or triple.We consider the satellite to be in synchronous rotation if its maximum yaw angle over the final 10 days of the simulation is < 60 • .If a synchronous rotator has a roll angle > 60 • over the final 10 d, we also consider it to be in the "barrel instability" (see Fig. 10).If a fragment consisting of 3 or more particles is ejected from the system, we count the simulation as successfully forming a pair.Triples are only counted if the system remains stable after further integration with REBOUND.Each category is not exclusive, in the sense that a triple could have also formed a pair, for example.
We find that a lower disk mass leads to a higher likelihood of the satellite forming with synchronous rotation, likely due to fewer violent collisions and mergers.For similar reasons, we find that more massive initial disks are more likely to form asteroid pairs as there is more mass (and angular momentum) that can launch fragments onto escape trajectories.Given that the final outcomes are highly dependent on the initial angular momentum in the disk, we caution against over-interpreting Fig. 23.A higher fidelity study with more realistic initial conditions would provide better estimates for the fraction of triples, pairs, and synchronous rotators that could plausibly form from a single mass shedding event.

CONCLUSIONS
We demonstrate that satellites can rapidly form by accretion following a single mass shedding event.The formation time is relatively fast, with the satellite reaching its terminal mass in a matter of days.Satellites tend to form on near-circular orbits within ∼4 primary radii, although 3-body interactions with other moonlets can sometimes put the final satellite onto wider and more eccentric orbits.Due to the satellite's accretion near the Roche limit, the satellite tends to be prolate in shape and is oftentimes in synchronous rotation.Many of these synchronous rotators, however, are in a rolling state about their long axis while remaining tidally locked, which would terminate or significantly weaken any BYORP evolution until this mode is dissipated ( Ćuk et al. 2021;Quillen et al. 2022).
The simulated secondary shapes are broadly consistent with the radar-derived shapes of other asteroid satellites, given their uncertainties (Ostro et al. 2006;Naidu et al. 2015;Becker et al. 2015).In addition, the distribution of secondary elongations (a/b) agrees well with the lightcurve-derived data set from Pravec et al. (2019).Our simulations rarely produce satellites with a/b ≳ 1.5 in strong agreement with the observed population of secondaries.Despite the strong preference to form elongated satellites, oblate shapes like Dimorphos are occasionally formed in our simulations.However, the fact that lightcurves are heavily biased against measuring satellite shapes with low a/b, combined with Dimorphos's low a/b suggests the existence of a significant population of oblate secondaries (unless Dimorphos is an outlier).This implies that there may be a longer term, yet-unidentified pathway that can reshape satellites to become oblate.
We tested four different friction angles (29 • 32 • , 35 • , 40 • ), finding only subtle differences.Generally, higher friction satellites have more irregular shapes.In addition, they are more resistant to disruption during close gravitational encounters and collisions, allowing them to form on wider, more excited orbits, whereas a lower friction satellite would disrupt and reaccrete during this process.The most noticeable difference is that higher friction tends to form satellites with lower bulk densities, due to their ability to maintain larger interior void spaces.
The accretion process is highly chaotic, resulting in a wide range of outcomes.In some cases, moonlets are ejected from the system via three-body interactions, effectively forming asteroid pairs.We note that this process does not explain the population of asteroid pairs in general, as there are many other processes that can form pairs. Here, we merely demonstrate that 3-body interactions may play some role in the production of asteroid pairs.
In a few percent of cases, we form stable triple systems, including satellites in mean motion resonances and coorbital satellites.Although unlikely, a future discover of a system with co-orbital satellites would be strong evidence for formation from a single mass shedding event, as co-accretion would naturally explain their configuration.
There are some caveats and limitations of our model.These simulations are only run for 100 days due to their computational cost and therefore cannot capture longer term effects.In addition, the primary is azimuthally symmetric which may affect the likelihood of forming pairs and triples.Finally, this study was originally focused solely on forming a Dimorphos-mass satellite, therefore we only explored a narrow range of disk masses.Although the initial conditions of these simulations were intentionally simplified, we would expect the general results and trends to hold.
There are a wide range of parameters available to explore in future work, and we list a few avenues to explore.With stronger spin-orbit coupling, an asymmetric primary may change the fraction of systems that form asteroid pairs, and would also allow the observed correlation between pair mass ratio and primary spin rate to be explored.In addition, this may affect the likelihood of forming a triple, as well as their long-term stability.In addition, a larger disk mass should be explored to probe the formation of systems with larger mass ratios such as Balam, and would also strongly affect the production of unbound fragments.Rather than a single mass shedding event, its plausible that these systems form from multiple events.If this is the case, simulated additional mass shedding events after the satellite is initially formed may change the distribution of secondary shapes.
The Hera mission will investigate the Didymos system for 6 months in early 2027 and will provide crucial information for binary formation models (Michel et al. 2022).In particular, it will measure the mass of Dimorphos and its internal properties, allowing us to check some of the predictions from our modeling.Through the landing and possible bouncing of its CubeSats on Dimorphos' surface, we will be able to access some of the mechanical parameters used in these simulations.However, since the shape of Dimorphos may have changed as a result of the DART impact (Raducan et al. 2023), we may not be able to put constraints on the formation models leading to either oblate or prolate secondaries, although all efforts will be done to relate the measured shape by Hera to the original one.
Adriano Campo Bagatin, Kate Minker, and John Wimarsson for helpful discussions and feedback.We are grateful to David Minton and the anonymous referee, whose feedback substantially improved the quality of this manuscript.

B. SIMULATION RESULTS
Table 3 shows the results of each simulation in this study.Note-Column M disk /MA is the initial disk-to-primary mass ratio, ϕ is the material's friction angle.The rest of the columns are the simulation outcome: the resulting secondary-to-primary mass ratio (MB/MA), the diameter ratio DB/DA, the secondary's axis ratios a/b and b/c (both DEEVE and extent-derived), the secondary's semimajor axis normalized to the primary's volume-equivalent radius (a orb /RA), the secondary's eccentricity (e orb ), orbit period (P orb ).θmax is the secondary's maximum libration angle over the final 10 days of the simulation.If this number is close to 90 • , then the secondary in an asynchronous or tumbling rotation state.N frag is the number of fragments ejected from the system.In cases where there is a second satellite, we include a second row in the table containing the parameters of the satellite.The final column specifies whether that satellite was stable after the handoff to REBOUND.If the satellite was unstable, it was either ejected from the system, or collided with the primary or secondary.

C. FINAL RENDERINGS
Here we show renderings of a subset of the simulations to give the reader a sense of the variance in satellite shape between similar simulations and as a function of the friction angle.All renderings are cases where M disk = 0.03M A and we show 8 cases for each value of friction (ϕ).Each frame is rendered from a view looking down from the body's principal rotation axis.ϕ = 29 • :

Figure 1 .
Figure 1.An example of a satellite forming after mass shedding from a Didymos-shaped primary.Particles are colored in white being shed from the primary.The initial spiral arms are caused by mass shedding occurring at localized regions rather than across the entire body.By ∼0.5 d, this asymmetry largely smooths out in azimuth before moonlets begin forming.A movie of this simulation is provided.

Figure 2 .
Figure 2. A time series plot showing the shape (a/b and b/c, based on the satellite's DEEVE), mass ratio (MB/MA), body separation (r orb /RA), and the attitude of the satellite formed via mass shedding.Within only several days, a Dimorphos-mass satellite is formed.The satellite is initially tidally locked to the primary (starting at ∼1.3 d), but then undergoes a partial tidal disruption followed by a merger with a moonlet, causing it to break from synchronous rotation (at ∼2.5 d).The satellite has an oblate shape, similar to Dimorphos.However, the satellite's shape, mass, orbit and rotational state will continue to change as it continues to accrete more material.

Figure 3 .
Figure 3.A view of the initial conditions from the disk simulations.The primary (gold colored particles) behaves as a single rigid body.

Figure 4 .
Figure 4.An example of a simulation where the secondary is formed in synchronous rotation.This is Disk 008 in Table 3, and has an initial disk mass of M disk = 0.02MA and a friction angle of ϕ = 35 • .A movie of this simulation is available.

Figure 5 .
Figure 5. Time-series plots of the example simulation from Fig. 4 (Disk 008).Starting from the top, we plot the satellites DEEVE axis ratios, mass ratio, orbital distance, and Euler angles, which describe its orientation in a frame rotating with the orbit.In only a handful of days, the satellite reaches a near-circular orbit at ∼3RA with a mass just under 0.01MA.The rotation state of the satellite is excited, but tidally locked, given by all three Euler angles being large but still well under 90 • .

Figure 6 .
Figure6.A comparison of orbiting particles between the more realistic, full spin-up simulation shown in Fig.1and the simplified initial conditions from the simulation shown in Fig.4.The orbital elements are plotted when moonlets first begin forming.For the mass-shedding simulation, this corresponds to ∼0.5 d after mass is initially shed, while for the simplified disk, this occurs around ∼4 d into the simulation.Although the distribution of orbital elements between the two cases does not match perfectly, due to the gravitationally unstable initial conditions of the simplified disk, we arrive at a similar range in semimajor axis, eccentricities, and inclinations seen in the more realistic simulation.

Figure 8 .
Figure8.The bulk density of all satellites with ϕ = 35 • at the end of the simulations.The density of all individual spheres is 3.5 g cm −3 , so these rubble piles have packing fractions of ∼0.62.Realistically, it is possible for the "real" bulk density to be slightly higher than what is measured here, as small particles which are below the resolution limit of these simulations would fill in some void space.

Figure 9 .
Figure9.The semimajor axis (in terms of the primary's radius) and eccentricity of the newly formed satellite, over all simulations with ϕ = 35 • .The colors indicate the starting disk mass, and the three points with error bars show the mean and 1σ standard deviation for the three different disk masses.Although the variance is quite large, a larger disk mass typically leads to a higher eccentricity and larger semimajor axis, due to the disk having a higher angular momentum and more collisions that can drive up the largest satellite's eccentricity.
Figure10.The maximum roll and yaw angles for each satellite, based on the final 10 days of the simulation.There are three distinct post-formation rotation states for the satellites, which are annotated on the plot to guide the reader's eye.We consider the satellite to be in synchronous rotation if its long axis remains aligned towards the secondary (i.e., yaw ≲ 60 • ).The "barrel instability" is a subset of the synchronous rotators where the secondary's long-axis stays largely aligned with the primary, yet it is able to roll about its long axis like a barrel.Finally, a significant fraction of satellites form in a tumbling state where all three Euler angles (roll, pitch, and yaw) can reach 90 • .

Figure 12 .
Figure 12.The satellite's shape, parameterized by its (a) DEEVE or (b) physical axis ratios a/b and b/c for the 96 simulations with ϕ = 35 • .The colored dots indicate the simulation results with varying initial disk masses.The shape and 1σ uncertainty of Dimorphos, along with the radar-derived shapes of the satellites of Moshup (1999 KW4), DP107, and SN263 are shown for comparison(Daly et al. 2023a;Ostro et al. 2006;Naidu et al. 2015;Becker et al. 2015).

Figure 14 .
Figure 14.Plots showing the effect of friction on the resulting satellite.(a) The final semimajor axis and eccentricity from each simulation.(b)The axis ratios of the satellite, based on its physical extents, as a function of friction angle.(c) The bulk density of each satellite as a function of its friction angle ϕ.Generally, bodies with a lower effective friction are able to achieve a more efficient packing and thus have a higher bulk density.

Figure 15 .
Figure15.Two examples that result in a Dimorphos-like shape (i.e., low a/b).(a) This satellite achieves an oblate shape due to a partial tidal disruption torquing the remaining mass onto a wide, eccentric orbit.Then, due to the body's asynchronous rotation, it accretes material in an azimuthally symmetric manner, as it sweeps through periapse with a different orientation every orbit.(b) This satellite achieves an oblate shape due to serendipitous mergers having a favorable geometry and is able to maintain its shape through the rest of the simulation.Movies of these simulations are provided Figure16.A rendering of two satellites that have similar shapes to Dimorphos, with the same lighting conditions and viewing geometry as DART's approach.

Figure 19 .
Figure19.Simulation results for all ejected fragments with ϕ = 35 • .These plots ignore all fragments containing less than 3 particles, of which there are many.The different colors show the initial mass of the disk.Only a minority of simulations eject fragments at all, but in some cases, multiple fragments are ejected in a single simulation.(a) The size of each ejected fragment as a function of both the fragment-to-primary and fragment-to-secondary diameter ratio.Unsurprisingly, a more massive initial disk ejects larger fragments.(b) A histogram of the spin rate of each ejected fragment.(c) The shapes of each ejected fragment, defined by their physical extents.(d) The eccentricity each fragment relative to the primary.Most fragments are on escape trajectories (e orb > 1) while the others have their apoapse outside of the systems Hill sphere if it were located at 1 au.
(a) Final frame from pkdgrav simulation (b) Time series showing capture into 2:1 MMR.

Figure 21 .
Figure 21.An example of a stable triple system formed from a single mass shedding event.(a) A snapshot of the system at the end of the simulation.(b) A time-series plot showing the satellite mass, semimajor axis, eccentricity, and resonance argument, ϕres.Capture into the 2:1 MMR enables the system to remain stable.
Figure 22.The two satellites capture into the 5:2 MMR around 15 d.Due to M2's high eccentricity, the two satellites have a close encounter, sending M1 inwards where it undergoes a tidal disruption around 30 d, which destroys both M1 and the resonance.Around 60 d, a new inner satellite begins growing and quickly captures again into the 5:2 resonance again with the more massive outer satellite, this time remaining stable.

Figure 23 .
Figure23.A summary of simulation outcomes showing the percentage of cases that have a given outcome, considering only the 96 simulations where ϕ = 35 • .A satellite is counted as a synchronous rotator if its yaw angle remains < 60 • over the final 10 d of simulation.Similarly, a satellite is considered to be in the barrel instability if its roll angle is > 60 • over the final 10 d.These thresholds are set arbitrarily and should only be considered notional.

Facility:
Simulations and analysis were performed on the ASTRA cluster administered by the Center for Theory and Computation, part of the Department of Astronomy at the University of Maryland and on Mésocentre SIGAMM hosted at the Observatoire de la Côte d'Azur.APPENDIX A. NOTATION For reference, Table2lists the notation used in this paper.

Table 1 .
Number of simulations for each combination of M disk and ϕ, along with the static friction coefficient (µS) and shape parameter (β) used to achieve the given friction angle (ϕ).

Table 2 .
Agrusa et al. (2021)ength is measured either based on the body's physical extents, or by the axis lengths of the corresponding dynamically equivalent equal-volume ellipsoid.dKeplerianelement, based on the instantaneous body positions and velocities.eTheangle between the secondary's long axis and the line connecting the body centers f The 1-2-3 Euler angle set coordinated in the secondary's rotating orbital frame (seeAgrusa et al. (2021)).
a M1 and M2 are used when there is more than one satellite.b Based on the body's volume-equivalent sphere.c

Table 3 .
Tabulated simulation outcomes for each simulation.Table 3 continued on next page

Table 3
(continued)Table 3 continued on next page