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

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 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 are highly sensitive to small changes in the primary's shape, meaning that each landslide, mass shedding event, and natural impact 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. ¹⁶ 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. (2008, 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 Sánchez & Scheeres 2012, 2016; Murdoch et al. 2015). 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 that 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 post-fission 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 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 reimpacts the primary, but it can also escape the system or reimpact 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, 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).

Jacobson & Scheeres (2011) suffer 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, 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), 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. 2024; Raducan et al. 2024). 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, 2011; Salmon et al. 2010). In this model, material migrates outward 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 from 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., Zhang et al. 2017, 2021; Yu et al. 2018). 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. 2024). 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.

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–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⁵ particles) and a highly polydisperse particle size distribution, the SSDEM spin-up simulations of Zhang et al. (2017, 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 nonspherical particles with irregular shape, where a 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 did not, 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's mass-shedding failure behavior under reasonable conditions.

We finally note that the mass shedding does not 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.

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–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 (SFD) with similar radii found on Dimorphos, as observed in the images taken by the camera on board the DART spacecraft (Pajola et al. 2024). The granular aggregate was configured based on the preimpact Didymos shape model (Barnouin et al. 2024) 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 hr, 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 time step of ∼0.02 s and a larger k_N ∼ 8 × 10⁷ kg s⁻². The bulk density of the body is set to 2.7 g cm⁻³, consistent with Didymos's latest bulk density estimate constrained by the updated shape model and the binary separation, i.e., 2.760 ± 0.130 g cm⁻³ (Naidu et al. 2024). ¹⁷ 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⁻³; therefore, the interparticle cohesive strength was set to 0. ¹⁸ 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 hr, 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 Dimorphos-mass 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 hr 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. 2024; Y. Zhang et al. 2024, in preparation).

Snapshots of this simulation are shown in Figure 1 along with a time-series plot of the satellite's shape, mass, orbit, and attitude in Figure 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 satellite's dynamically equivalent equal-volume ellipsoid (DEEVE). ¹⁹ 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 inward, 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 in Zenodo at DOI: 10.5281/zenodo.8387043.

Zoom In Zoom Out Reset image size

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.

Comparison with Sugiura & Kobayashi (2021) and Hyodo & Sugiura (2022). Recently, Sugiura & Kobayashi (2021) studied the shape evolution of a rubble pile under YORP spin-up using smoothed particle hydrodynamics (SPH) simulations. In contrast to discrete element method 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 & Kobayashi (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 & Kobayashi (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., Mitchell et al. 1972; Bareither et al. 2008; Al-Hashemi & Al-Amoudi 2018) as well as friction angle estimates of recently visited asteroid surfaces (also on the order of ∼35°; e.g., Fujiwara et al. 2006; Barnouin et al. 2019; Watanabe et al. 2019; Barnouin et al. 2024; Robin et al. 2024). This is because Sugiura & Kobayashi (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 (\phi )p\,+\,c$, where ϕ is the friction angle, p is the confining pressure, and c is the cohesion. Sugiura & Kobayashi (2021) instead adopt a material model of the form $Y=\tan (\phi )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 & Kobayashi (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 & Kobayashi (2021) is not unreasonable (Barnouin et al. 2019; Arakawa et al. 2020; Jutzi et al. 2022; Walsh et al. 2022). 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 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 & Kobayashi (2021) in a critical way. Both methods apply a similar angular acceleration to the rubble piles (∼10⁻¹⁰ rad s⁻²), which ensures that the artificially induced Euler acceleration is always negligible compared to the centrifugal acceleration. However, Sugiura & Kobayashi (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 & Kobayashi (2021) is small, this can make a critical difference for a rubble pile on the edge of stability. As a result, Sugiura & Kobayashi (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, 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 timescale 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 of 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 postmerger 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⁻³ 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⁻³, 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, and other binary satellites needs to self-consistently address the spin-up, mass-shedding, accretion, and (potentially) merger processes. For example, the shape of a postmerger 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.

The simulation shown in Figures 1 and 2 are computational expensive, requiring ∼10⁵ particles, and are therefore limited to short integration times and 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., Walsh et al. 2012; McMahon et al. 2020; Sánchez & Scheeres 2020; Zhang et al. 2022). 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 that 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, an equatorial radius of a = b = 400 m and a/c = 1.5, and a bulk density of 2.4 g cm⁻³, similar to the best estimate for Didymos's shape and density at the time this study began (Daly et al. 2023). 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. ²⁰ 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 Figure 3, the radii of the shell particles are overinflated 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 hr in all simulations, approximately matching the rotation period of the primary from Section 2.1 following its rotational disruption.

Zoom In Zoom Out Reset image size

The orbiting material is generated following a power-law 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 exactly match) the observed boulder SFD on Dimorphos (Pajola et al. 2022, 2023), as well as the boulder SFD seen on other rubble-pile asteroids (e.g., Dellagiustina et al. 2019; Michikami et al. 2019; Burke et al. 2021; 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. 2023). 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 in 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. 2024).

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⁻³, to match the grain density of L and LL chondrites (Flynn et al. 2018), which are the best meteoritic analogs 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 ,

$$\begin{eqnarray*}\begin{array}{rcl}{n}_{i}^{2} & = & \displaystyle \frac{{{GM}}_{A}}{{r}_{i}^{3}}\left(1+\displaystyle \frac{3}{2}{J}_{2}{\left(\displaystyle \frac{{R}_{A}^{\mathrm{eq}}}{r}\right)}^{2}\right)\\ & + & \displaystyle \frac{{M}_{\mathrm{disk}}(\lt {r}_{i})}{{r}_{i}^{3}},\end{array}\end{eqnarray*}$$

where M_A , R^eq _A , and J₂ 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 demonstrate below rapidly leads to particle orbits representative of the more realistic starting conditions following a mass-shedding event, such as that shown in Figure 1.

Any particles that reach a distance of 40 km ($100{R}_{A}^{\mathrm{eq}}$) 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 ∼two-thirds 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 not have returned to the binary system. Therefore, deleting these particles has a negligible effect on the binary's formation. In 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 satellite's formation. For each disk mass (M_disk) and friction angle (ϕ), we run a given number of "clone" simulations, which are listed in Figure 1. In total, we run 120 simulations. Ninety-six of the simulations have ϕ = 35° and 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 days (5.8 × 10⁷ steps) due to the high computational cost imposed by the small time step. Each simulation is assigned a number between 1 and 120, and the result of each simulation is tabulated in Table 3 in the Appendix.

First, we show a representative simulation to demonstrate the model. In Figure 4, we show eight snapshots of a simulation in which M_disk = 0.02M_A and ϕ = 35°. Time series plots from this simulation are shown in in Figure 5. 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 of outcomes in terms of the satellite's physical and dynamical properties. Therefore, this simulation is not necessarily typical of all cases.

Zoom In Zoom Out Reset image size

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. Figure 6 compares the orbital elements of orbiting material between the more realistic mass-shedding example from Figure 1 and the simplified case shown in Figure 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 coplanar 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., Figure 1), although the more realistic initial conditions have a slightly wider distribution in both inclination and eccentricity.

Zoom In Zoom Out Reset image size

In Figure 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 $\tfrac{{D}_{B}}{{D}_{A}}\approx {\left(\tfrac{{M}_{B}}{{M}_{A}}\right)}^{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 Dimorphos-to-Didymos size ratio is ∼0.2, which corresponds to a mass ratio of ∼0.01 if the bodies have equal densities (Daly et al. 2023). 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 are 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., Zhang et al. 2018; Yu et al. 2019; Hyodo & Sugiura 2022; Zhang et al. 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 outward, which substantially increases the timescale for the satellite's formation.

Zoom In Zoom Out Reset image size

We plot a histogram of the bulk density of each satellite at the end of the simulation in Figure 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 that 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.

Zoom In Zoom Out Reset image size

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. Figure 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 Figure 9, where we show the mean a_orb and e_orb along with their standard deviations for each of the three different disk mass cases.

Zoom In Zoom Out Reset image size

For each satellite, we determine its rotation state based on the 1-2-3 Euler angle set, consisting of the satellite's 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, nonplanar rotation states, necessitating the use of the Euler angle convention. In Figure 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.

Zoom In Zoom Out Reset image size

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 that 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 of gravitational accumulation near the Roche limit, such as the accretion of moonlets from Saturn's rings (Karjalainen & Salo 2004) and circumplanetary disks (Hyodo et al. 2015).

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 axes, 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 extent of the body along these axes. The second measure is the axis lengths of the 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 do not 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 Figure 11, where the DEEVE axis ratio tends to be much larger than the physical extent axis ratio for large axis ratios.

Zoom In Zoom Out Reset image size

Due to the discrepancy between the DEEVE- and extent-derived shapes highlighted above, we compare both quantities to known secondary shapes for thoroughness. In Figure 12, we show both the DEEVE- and extent-derived semiaxis 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₄), 2000 DP₁₀₇, and 2001 SN₂₆₃, which are the only publicly available radar-derived shapes for the satellites of other binary asteroids (Ostro et al. 2006; Becker et al. 2015; Naidu et al. 2015), as well as the axis ratios of Dimorphos (Daly et al. 2023, 2024). The updated shape of Dimorphos provided by Daly et al. (2024) differs slightly from the initial assessment in Daly et al. (2023), 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₁₀₇ are the best comparisons with Dimorphos, as they are both S-type binaries, whereas SN₂₆₃ is a C-type triple (with large uncertainties in the satellite shapes). Figure 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.

Zoom In Zoom Out Reset image size

We also compare our results to the lightcurve-derived shapes from Pravec et al. (2016, 2019). Figure 13 shows a histogram of the simulated a/b ratio compared to the data set 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.

Zoom In Zoom Out Reset image size

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 with high a/b has been attributed to the chaotic dynamics of satellites with $a/b\gtrapprox \sqrt{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 do not 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/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.

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 eight 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 Figure 14(a), 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.

Zoom In Zoom Out Reset image size

The axis ratios (based on the physical extents) are given in Figure 14(b). In general, there is a smaller variation in the secondary's shape as the friction is reduced. In addition, 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 friction are able to maintain a wider range of shapes. In the limit that the friction angle was reduced to 0°, the satellite effectively behaves like a fluid and takes 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 Figure 14(c), 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⁻³, 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 higher.

Although these simulations demonstrate that immediately forming an oblate satellite is relatively rare, we show two cases with ϕ = 35° that achieve shapes similar to Dimorphos. Although neither case is a perfect match, they do suggest that an oblate spheroid-like shape can plausibly be formed from a single mass-shedding event. In Figure 15, we show the two examples, both of which have an initial disk mass of M_disk = 0.04M_A . In the first example (Figure 15(a)), 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 inward 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 days, which abruptly increases b/c. After the merger, the satellite continues accreting on all sides, leading to a gradual reduction in a/b and increase in b/c, until there is almost no material left to accrete. 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.

Zoom In Zoom Out Reset image size

In the second example (Figure 15(b)), the satellite achieves its oblate shape simply by undergoing two near-head-on mergers with other moonlets at early times (around ∼3.5 and ∼7 days, 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 Figure 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.

Zoom In Zoom Out Reset image size

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.

3.7.1. 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. 2010, 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, 2019). There is a distinct correlation between the 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 satellite's 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 Figure 17. In this case, the smaller satellite (M₂) has a series of close encounters with the larger satellite (M₁), 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 Figure 18, we show snapshots of this simulation at the moment when M₂ is ejected. The two satellites have a close encounter at 48.3 days, where M₁ torques M₂ onto a hyperbolic orbit, ejecting it from the system. Meanwhile, as a consequence of angular momentum conservation, M₁'s eccentricity is raised (i.e., its periapse distance decreases), causing it to have a close approach with the primary at 48.9 days. M₁ then undergoes a partial tidal disruption, losing some mass, reshaping, and torquing it back onto a higher orbit safe from further tidal disruption. This process also breaks M₁ from its previous synchronous rotation state.

Zoom In Zoom Out Reset image size

Download figure:

Video Standard image High-resolution image

In Figure 19(a), 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 three or more particles. Also, this figure only includes fragments that have been ejected within the ∼100 days of pkdgrav integrations. However, as we will see in the next section, there are several systems that still contain two satellites after 100 days 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 (Figure 19(b)), 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 and 9 day⁻¹ (Pravec et al. 2019). In addition, many of the pairs formed in these simulations are in non-principal-axis 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 Figure 19(c) 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 NEAs, these satellites would be lost anyway, 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.

Zoom In Zoom Out Reset image size

For example, asteroid (3749) Balam is a main-belt triple with a fast-rotating (2.8 hr 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 backward numerical integrations indicating that asteroid 2009 BR₆₀ separated from Balam less than 1 million yr 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 milder 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.

3.7.2. Triples

In some cases, we find that three-body interactions lead to the formation of systems with two satellites on noncrossing 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₂ 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³ yr 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₂₂ 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 days of pkdgrav simulations, of which seven 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 an MMR is effectively random. Out of the seven stable triples, five are in MMRs, while two are not in resonance.

In three of the resonant cases, the two satellites form in a coorbital configuration (i.e., the 1:1 MMR) by coaccretion. When this occurs, the smaller satellite occupies the L₄ or L₅ Lagrange point of the larger body. Snapshots of the three cases are shown in Figure 20. These systems satisfy the condition for stability, in which $\tfrac{{M}_{1}+{M}_{2}}{M+{M}_{1}+{M}_{2}}\lesssim 0.04$, where M is the primary's mass and M₁ and M₂ are the two satellite masses (Murray & Dermott 2000; Laughlin & Chambers 2002; 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³ yr. We restrict the REBOUND simulations to only 10³ yr, because we neglect tides and nongravitational forces that can tighten or loosen the tadpole orbits, depending on the direction of migration (e.g., Fleming & Hamilton 2000). We find that coorbital 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 a circumplanetary disk.

Zoom In Zoom Out Reset image size

We note that the idealized initial conditions of these simulations may increase the likelihood of forming coorbital satellites and that the primary's azimuthal symmetry likely contributes to their long-term stability. Here, we merely point out that forming a coorbital system via a single mass-shedding event is possible, albeit unlikely. A future detection of a triple asteroid with coorbital satellites, however unlikely, could be explained by coaccretion 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 Figure 21, we show an example of capture into the 2:1 MMR, where the inner satellite is the most massive. In Figure 21(a), we show a top-down view of the triple system at the end of the 100 days of pkdgrav integration along with a time-series plot in Figure 21(b). Early on, both M₁ and M₂ 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₂) onto a wide orbit. Then, M₂ continues to grow and migrate outward through a series of mergers and close encounters with other moonlets. Meanwhile, a new inner satellite forms (now M₁) from one of the disrupted remnants and migrates inward. This process places M₁ and M₂ close enough to the 2:1 MMR that they capture. The resonance is then further stabilized by the smooth growth of M₁, the inner satellite. We can confirm the existence of the resonance due to the libration of the resonance argument ${\phi }_{\mathrm{res}}=2{\lambda }_{2}-{\lambda }_{1}-{\varpi }_{2}$, where λ_i are the respective body's mean longitudes and ϖ₂ is M₂'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 nonadiabatic. We confirmed that the resonance argument continues to librate when handed off to REBOUND.

Download figure:

Video Standard image High-resolution image

We also find one example of capture into the 5:2 MMR, shown in Figure 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.

Download figure:

Video Standard image High-resolution image

We find two other long-term stable triples (disks 075 and 085) with a more massive outer satellite having period ratios close to 2:1 and 4:1, respectively, although we confirm that they are not in resonance. However, these systems have well-separated 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 MMRs. Out of the NEAs, there are three confirmed triples: 1994 CC, 2001 SN₂₆₃, and 3122 Florence (Nolan et al. 2008; Brozovic et al. 2009). Both CC and SN₂₆₃ are well characterized, not near MMRs, and also have some mutual inclination, which has 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 and 7 hr (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 an MMR, 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 outward, followed by convergent migration (due to tides and/or BYORP), is an equally plausible mechanism to form a MMR triple.

3.7.3. Simulation Outcomes

In Figure 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 days, we also consider it to be in the "barrel instability" (see Figure 10). If a fragment consisting of three 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.

Zoom In Zoom Out Reset image size

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 overinterpreting Figure 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.