High-fidelity Modeling of Rotationally Fissioned Asteroids

Binary asteroids represent an important aspect of the dynamical evolution of small bodies and may provide insight into the evolutionary history of these populations as a whole. Many past studies have focused on Yarkovsky–O’Keefe–Radzievskii–Paddack-driven spinup and disruption as a pathway for binary formation in the inner solar system. While these studies have shown the likelihood that such a process occurred, they are generally limited by assumptions and simplifications in their dynamics models. In this study we apply a high-fidelity and computationally efficient model of binary asteroid dynamics in order to understand the potential effect of higher-order gravity terms and nonplanar dynamics on binary fission. We apply this dynamics model to 66391 Moshup (1999 KW4), 8567 (1996 HW1), and 185851 (2000 DP107) as a representative set of binary and contact binary systems to understand the implications for their fission and formation. Our analysis supports the importance of secondary fission for stable low-mass binary formation as initially suggested by Jacobson and Scheeres. Additionally, we find that the inclusion of higher-fidelity dynamics distorts the dynamical structure of the system, creating a pathway for the secondary to re-collide with the primary. The increased complexity from the inclusion of nonplanar dynamics also suggests more excited spin states of the asteroids during disruptive events such as secondary escape and fission.


Introduction
The solar system's asteroid populations contain a dynamically rich variety of systems. It is estimated that 16% of observed asteroids are binary asteroids, in addition to roughly 70 candidate asteroid pairs . These dynamic objects are widely believed to be produced by rotational fission of rubble pile asteroids driven by the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect (Walsh et al. 2008;Jacobson & Scheeres 2011). A number of past studies have supported the plausibility of the YORP effect as the driving evolutionary process that led the asteroid population to its current dynamical richness. In general the proposed evolutionary pathway begins with YORP-driven spinup of a rubble pile asteroid. Eventually, the asteroid spin rate exceeds some critical value, leading to mass-shedding or outright fission (Scheeres 2007;Walsh et al. 2008). After fission the expelled mass coalesces into a secondary asteroid, which either settles into a stable orbit about the primary asteroid, re-collides with the parent body, or escapes from the system and forms an asteroid pair.
This pathway was initially proposed as an explanation for observations of a wide range of asteroids by Pravec et al., who show an upper spin limit of about 2.2 hr throughout the observed population (Pravec et al. 1998;Pravec & Harris 2000). In response, theorists began to explore YORP spinup as an explanation. Work by Scheeres (2007) explored YORP spinup as a driver for contact binary fission, while Walsh et al. (2008) studied YORP-driven rubble pile mass-shedding as a formation process for binary asteroids. Scheeres focused on contact binary fission, assuming that the asteroid was initially two distinct bodies resting on one another. On the other hand, Walsh et al. treated the asteroids as a collection of small spherical masses that would be shed from a single asteroid and then coalesce into the secondary asteroid. Both of these works showed the potential of YORP spinup to explain both the observed spin barrier of near-Earth asteroids (NEAs), as well as the dynamical complexity of observed NEA systems. However, the binary systems formed by both theories showed significant energetic excitement compared with the energetically relaxed binaries commonly observed among NEAs. This required a powerful source of energy dissipation to match these predictions to the observed population. Work by Jacobson & Scheeres (2011 explored the potential sources of this energy dissipation. They modeled the postfission dynamics for the contact binary case with a simplified gravity model, the effects of solar gravity, tidal torques, and the cohesive strength of the secondary asteroid. Their work identified secondary fission, or structural failure of the secondary, driven by mutual momentum transfer between the asteroids as a powerful source of energy dissipation for systems with a low mass ratio (0.2). For systems with a larger mass ratio the study suggested a slower relaxation into the doubly synchronous equilibrium, eventually evolving into a contact binary or asteroid pair. With these two processes their work was able to describe a complete evolutionary map for NEA binaries, asteroid pairs, and other NEA systems (JS2011). However, their analysis ignored both nonplanar effects and higher-order gravity terms. Work by Asphaug & Scheeres (1999) had employed a fictitiously disrupted model of 4769 Castalia under nonplanar dynamics to theorize the potential of nonplanar interactions as a source of non-principal axis rotation in asteroid populations. Scheeres et al. (2000) also suggested nonplanar interactions between two asteroids and between asteroids and planets as a potential source of the observed nonprincipal axis rotation in the NEA population. Later work by Boldrin et al. (2016) did begin to investigate the effects of nonplanar and second-order dynamics in binary fission, but they limited their investigation to the formation of contact binaries. Their analysis hypothesized nonplanar binary fission as a source of non-principal axis rotation found among NEAs.
Since the JS2011 analysis, the computational efficiency of binary asteroid dynamics models has improved significantly. The standard dynamics model for binary asteroids is the full two-body problem (F2BP). The F2BP treats both bodies as an arbitrary mass distribution. The bodies influence each other's motion through the mutual gravity potential, a function of their relative position and orientation that is used to model both gravity accelerations and torques. This leads to a coupling between the rotational and translational motions of each body, giving the dynamics nine degrees of freedom for the nonplanar problem (Maciejewski 1995). Recent work by Hou et al. (2017) reformulated the mutual gravity potential as a recursive summation wherein the mass distributions are described using the inertia integrals. The authors have developed an opensource simulation tool for binary asteroids building upon this model. A description of both the model and the dynamics tool can be found in the appendices and the cited repository (Davis 2019). We apply this dynamical model to evaluate the conclusions of past work under a more rigorous dynamics model.
Of particular interest for this study are the expected formation rates and conditions for the most common multibody categories of asteroid system. The categories of interest are asteroid pairs, contact binaries, low-mass-ratio binaries, and high-mass-ratio binaries. In this nomenclature asteroid pairs are defined as a collection of two or more bodies that can be shown to have similar orbital elements and to have had small relative velocities and separations in the past. It has been suggested that asteroid pairs form after the fission of a parent body similar to observations of P/2013 R3 (Jewitt et al. 2017). Contact binaries are defined as bi-lobed asteroids that appear to be two distinct bodies resting on one another. Examples are 25143 Itokawa and 8567 (1996 HW1) (Fujiwara et al. 2006;Howell et al. 2010). The binary asteroid population of the NEAs is commonly separated into low-mass-ratio and high-mass-ratio binaries. Low-mass-ratio binaries are defined to have a mass ratio below 0.2, but tend toward a mass ratio of 0.05 with a topshaped primary, similar to 66391 Moshup and 65803 Didymos (Ostro et al. 2006;Naidu et al. 2020). The majority of lowmass-ratio binaries also exhibit a singly synchronous arrangement where the smaller secondary is tidally locked and the primary has a nearly constant spin rate, faster than the secondary's orbital rate. High-mass-ratio binaries, on the other hand, tend toward bodies of similar size, such as 1994 CJ1 and 2017 YE5 (Taylor et al. 2014(Taylor et al. , 2018. These systems are generally in the doubly synchronous equilibrium where both bodies are tidally locked. In addition to these multibody asteroid categories, we are also interested in the source of the equatorial craters seen on 185851 (2000 DP107) and 341843 (2008 EV5) and the potential for these to be the result of a fission event (Tardivel et al. 2018). While both JS2011 and Walsh et al. provided constraints on the formation of these various objects, the aim of this work is to further constrain and explore the feasibility of these formation processes.
Our study focuses on the evolution of binaries from the moment of fission until the secondary collides with the primary, escapes, or experiences secondary fission. To model these dynamics we apply the F2BP with the inclusions of Hill relative solar gravity effects and mutual tidal torques. To model structural failure of the secondary a simple cohesion model of a spherical body is used to compute a spin rate at which failure occurs. This model is used to evaluate the evolution of three NEA binaries, which sample the the different classifications of this subpopulation. The binaries we model are Moshup, 1996 HW1, and 2000 DP107 (Figure 1, Table 1). These example systems represent low-mass-ratio binaries (Moshup and DP107) and high-mass-ratio binaries (HW1). We also use DP107 to test the feasibility of the proposal by Tardivel et al. that NEA equatorial craters may be markers of binary fission. For all three systems, the initial fission geometry (Figure 1) is analogous to that of a fissioning contact binary. The table collects the estimated physical parameters of these binaries based on current radar observations. The mass parameters in the table are computed assuming a constant density of the asteroid radar shape models. The densities for Moshup and DP107 have been estimated as 2.0 and 1.38 g cm −3 (Ostro et al. 2006;. HW1 does not have an estimated density so it is assumed to be 2.0 g cm −3 based on convention. The mass fraction used in the table is defined as n = + M M M 2 1 2 ( ).

Dynamics Models
The goal of our dynamics model is to capture the higherorder effects of the mutual gravity interactions as well as other significant perturbations. We use the model of the F2BP of Hou et al. (2017) to capture the mutual gravity interactions of the bodies. This model captures nonplanar effects and gravity terms up the fourth order for both bodies (Hou et al. 2017). It also assumes the bodies to be rigid, unlike previously mentioned work by Walsh et al. and work by others since (Walsh et al. 2008;Ferrari et al. 2017). The details of the mutual gravity potential used in this study are presented in Appendix A. Because our study is concerned with the dynamical evolution as opposed to the fission mechanisms at work, we do not include the effects of non-rigidity beyond tidal torques. In order to account for the mutual tidal torques between the bodies we apply the model described by Murray & Dermott (1999;see also Goldreich & Gold 1963). This model treats the tidal torque as the effect of a tidal bulge on a spherical body driven by an orbiting point mass perturber. While this model ignores the effects of shape on the tidal torques, it is able to approximate the dynamical effect to first order without requiring specific knowledge of cohesive or structural properties. To be comparable to past work we use the value of the standard tidal Love number, k=10 −5 , and quality factor, Q=100 (Goldreich & Gold 1963;Jacobson & Scheeres 2011). Our solar gravity model builds on the relative Hill acceleration models described in JS2011 as well as Fahnestock & Scheeres (2006). The form of the model we use assumes the binary's barycenter to be on a fixed circular orbit about the Sun and computes a relative solar gravity acceleration between the asteroids.
In addition to the integrated dynamics, our post-processing analysis evaluates the dynamics for secondary fission. To identify secondary fission, we compute a critical spin rate for cohesive failure and check the secondary's spin rate throughout a simulation. The critical spin rate is computed as a function of elongation, α, density, ρ, strength, σ, angle of friction, f, and the gravity constant, G (Holsapple 2001(Holsapple , 2004: To determine the angle of friction we refer to the values studied by Sánchez & Scheeres to analyze cohesive failure, ranging from 45°to 90° (Sánchez & Scheeres 2014. Using values from either end of this range, our analysis does not see differences in the critical spin rate beyond fractional percentages for our simple model.

Simulation Conditions
The goal of the simulation for all three systems studied is to evaluate their evolution and understand the influence that fission conditions have on their evolved state. Each of the three targets is selected to explore the evolution of different system types: Moshup represents current low-mass-ratio binaries, 1996 HW1 represents contact binaries and high-mass-ratio binaries, and 2000 DP107 further explores the regime of low mass ratio while probing equatorial crater formation of the secondary. For each system we identify the conditions for its fission based on the body shapes and equilibria. A Monte Carlo analysis is then performed with initial conditions perturbed from the equilibrium and simulated for one year. For simplicity it is assumed that all fissions occur in the ecliptic plane. Each simulation is then post-processed to evaluate the fate of the fissioned binary. The fates fall into four categories: collision, escape, secondary fission, and chaotic orbital capture. Collision is assumed to occur when the system geometry aligns such that the best-fit ellipsoids for each of the two bodies intersect. Escape occurs when the secondary crosses the sphere of influence of the primary's gravity. Secondary fission occurs when the secondary spin rate exceeds the critical spin rate as defined in Equation (1). The numerical integration is performed with an adaptive Dormand-Prince RK7(8), while collision detection, escape, and secondary fission are evaluated in post-processing (Prince & Dormand 1981).

Fission Conditions
We define the simulations to begin from the moment of fission. Within the framework of the F2BP a binary system must gain enough energy from YORP or other processes to pass through the inner unstable equilibrium of the F2BP. The inner unstable equilibrium acts as an energetic barrier between collision and chaotic orbital capture (Scheeres 2007(Scheeres , 2009). This barrier is identified by fixing the angular momentum of a binary system and varying the relative separation. When the bodies are separated by a few percent of the primary's radius a local maximum occurs (Figure 2). This maximum acts as the dynamical barrier to fission. In our Monte Carlo analysis we use perturbations from this ideal fission condition to account for the potential increase in the fission spin rate due to weak cohesive forces.
Using the methodology for computing F2BP equilibria described in our previous work, we evaluate the inner unstable equilibrium at the fourth order (Davis & Scheeres 2020). The angular momentum value for the calculation is selected as the circular orbit rate for near-contact geometry and depends on the system of interest. Our evaluation of the equilibrium at the fourth order finds that the symmetric alignment of the principal axes present for the second-order models breaks down. In Figure 3 we illustrate the relative geometry of the unstable equilibrium for the second-and fourth-order models of the Moshup system.  Note.a orbit identifies the semimajor axis of the binary orbit. a, b, and c are the best-fit ellipsoidal semi-axes for each body, respectively the long, intermediate, and short axes. The top value for each semi-axis give the primary's best-fit ellipsoidal dimension, while the bottom value corresponds to the secondary.
This change in alignment, combined with the higher dimensionality of the nonplanar problem, warps the zerovelocity curves at the equilibrium such that the equilibrium no longer acts as a boundary to collision. The change in alignment is caused by the higher-order gravity terms associated with the asymmetric shape torquing the secondary away from axial alignment. The asymmetry has no effect when the mutual gravity is truncated at the second order because second-order terms can always be rotated into symmetric principal alignment.
The Monte Carlo conditions for each system are seeded at this equilibrium and perturbed uniformly in their relative position, velocity, fission spin rate, and relative orientation. Perturbation of the relative position and velocity are directed only along the relative separation vector to avoid aliasing with the effects of the spin rate and orientation perturbations. The perturbations of spin rate maintain the direction of the system angular velocity at fission, but vary the magnitude. Relative orientation is perturbed by rotating the orientation of the secondary in all three dimensions relative to the equilibrium state. For perturbations of the position, velocity, and spin rate the values are only increased to avoid reducing the energy of the system below the minimum for mutual orbit. To maintain consistency across the three sample systems we compute the perturbation range as a relative percentage of the computed equilibria.

Moshup (1999 KW4)
Asteroid Moshup, previously known by its provisional name 1999 KW4, is an asymmetric NEA binary that has been extensively studied since its initial discovery. Radar measurements and shape modeling were performed by Ostro et al. and adopted for this study. The system characteristics are reported in Table 1. We select the Monte Carlo perturbation to be relatively small, but account for the potential effects of cohesion. The selected initial conditions and perturbations are provided in Table 2.  (2009) that depicts the energy of a second-order planar ellipsoid-ellipsoid binary as a function of the primary's radius for a fixed angular momentum. The local maximum at a radius value near 1 is the unstable inner equilibrium. For a conservative system a binary initialized under these conditions must remain between the red zero-velocity curve and the blue energy line. Thus the binary can either collide left of the equilibrium or be chaotically captured right of the equilibrium with the potential to further separate and escape (Scheeres 2009). . Relative geometry of inner unstable equilibrium evaluated at order 2 (above) and order 4 (below). The view for both is isometric, looking 45°d ownward and perpendicular to the relative separation vector. The relative twist caused by the higher-order terms results in a breakdown of the equilibrium as a dynamical barrier. We refer to the larger, gray body as the primary and the smaller as the secondary.
Of the 150 simulated cases the majority either collide or experience secondary fission, with a subset of secondaries escaping and a handful remaining in a chaotic captured orbit. Table 3 provides the statistical breakdown of these results along with the median time for the different evolutionary events to occur.
A few immediate comparisons can be made between these results and those of JS2011. Most significant is the occurrence of re-collision after the moment of fission. Re-collisions appear to be limited to within the first two days after fission, with the majority occurring within the first half day. In the case of secondary fission we can also see that our rate of secondary fission, 36.7%, is consistent with the rate of 40%±4% seen in JS2011. This comparison is not one-to-one, however, because our model for secondary fission accounts for cohesion. JS2011 has no cohesion, but assumes full disruption of the secondary. Finally, we see that the rate of chaotic capture is almost identical to the 2% seen in JS2011. The clearest point of disagreement is the median time for each event, with escape and secondary fission taking roughly 5-10 times as long to occur. This slowing of the evolution and disruption of the system results from the slower angular momentum transfer from the primary to the secondary under the higher-fidelity dynamics. JS2011 in fact predicts this.
Looking more closely at the re-collision events, we can better understand the nature and potential result of the collision. The two key factors of interest here are the location of the impact on the primary and the surface velocity of the impact. We define the surface velocity as the relative velocity between the surfaces of the primary and secondary at the moment of collision, as opposed to the relative velocity between the centers of mass. In Figure 4 the locations of collision on the Moshup primary are identified along with the magnitude of the surface velocity at impact.
It is clear that the collisions bias toward the location of maximum extent. This occurs because that location extends roughly 20 m further from the center of mass than the rest of the equatorial bulge. Thus secondaries orbiting just above the surface of the primary will hit this location first. The surface velocity for these collisions ranges from about 20-90 mm s −1 . This would be a slow grazing collision, unlikely to cause significant fracturing of either body. For context a fall from 0.1 mm on Earth will have an impact velocity of about 50 mm s −1 . We further probe the impact conditions by mapping the directional distribution of the surface velocity at impact ( Figure 5).
Note. We sample uniformly between the upper and lower perturbation bounds. The relative Euler angles define the relative orientation between the primary and secondary.  The directional distribution is described using the instantaneous radial, cross-track, and normal directions at the moment of impact. The radial direction is defined from the center of mass of the primary toward that of the secondary. The crosstrack direction is defined along the direction of the orbit. The normal direction is the orbit normal. The range of the radial impact velocities, the color in the figure, is relatively smaller than in the other directions and almost entirely negative, or toward the primary. While positive radial velocity at impact may not be intuitive, it results from combining the rotation of each body and their relative velocity into the surface velocity. One can imagine that the secondary may be moving away from the primary but rotating from an orientation with the longest axis pointing along the orbit to one with it pointing toward the primary. This would swing the long axis of the secondary into the surface of the primary. The wide range of cross-track and normal velocities for the impact also implies that these impacts are not only slow but likely would be grazing or rolling collisions. Likely this means that simulating the results of any potential re-formation of the system would require sophisticated models of the cohesive forces and continuum mechanics at play.
The most interesting aspect of secondary fission in our model is the out-of-plane behavior of the secondary at the moment of secondary fission. In past studies the dynamics have been assumed to be planar, such that all debris or tertiary bodies remain within the plane. Our model shows that secondary fission is likely to be a much more complex process. To illustrate this we define the dynamic inertia, I D , and effective spin of the secondary, ω l , at the moment of secondary fission as where H and T are respectively the angular momentum and kinetic energy of the secondary. Figure 6 illustrates the effective spin and dynamic inertia for the fissioning secondaries of the Moshup system superimposed over the secondary's moments of inertia. The spin state of the primary is not shown because it shows virtually no change in its rotational state. The dynamic inertia of the bodies is clustered about the intermediate axis, implying non-principal axis rotation or separatrix motion. This means that any fissioned debris is likely to either enter an excited out-of-plane orbit or collide with the primary randomly across its surface.
In comparing the JS2011 secondary fission behavior to this analysis several differing assumptions must be pointed out and compared. The JS2011 fission model assumes the secondary to be a contact binary made up of two similarly sized spheres without cohesion. In their model the fission of the secondary can be seen as a balance between the spin of the secondary, the mutual gravity between the two spheres making up the secondary, and the effects of the primary's gravity on each lobe of the secondary (Jacobson & Scheeres 2011). Our model assumes the material expelled from the secondary to be much smaller than the secondary, thus the gravity of the secondary alone dominates and cohesive forces have a significant effect. While we do not presume which model is the more accurate, there is value in comparing the results of the two models. To implement the JS2011 secondary fission model we define the mass fraction of the components making up the secondary where M a is the mass of the sphere initially closer to the primary and M b is the mass of the sphere initially further from the primary. The centers of each sphere are oriented such that the bi-lobed shape approximates the secondary's attitude. The bi-lobed secondary is only used within the evaluation of secondary fission and does not affect any other aspect of the dynamics model or later analysis. In Table 4 we present a comparison of the Moshup results processed with our cohesive secondary fission model as well as the JS2011 bi-lobed secondary fission model given a range of values for ν s .
The first comparison to be made is between the cohesive case and the isolated bi-lobed case. The isolated case is created Figure 5. Instantaneous surface velocities of colliding secondaries at the moment of impact for Moshup simulations. The velocities are provided in the frame of the secondary's orbit where the radial direction is from the center of mass of the primary toward that of the secondary, the cross-track direction is along the orbit in the direction of the orbital motion, and the normal direction is along the orbit normal. The surface velocity is the relative velocity between the surfaces of the primary and secondary at the moment of impact, accounting for relative spins and velocity.
for more direct comparison between the two models. It assumes the secondary's fissioning lobe to be sufficiently small, ν s =1.5×10 −5 or about 5 m in radius, and the primary's gravity to be ignorable. The two cases have fairly close median times for secondary fission and see nearly identical secondary fates. Once the gravitational influence of the primary is considered in the lobed model, the results reflect the importance of the secondary mass ratio. Comparing the two ν s =1.5×10 −5 cases, we see that the primary's gravity slows the secondary fission process. Further, if we compare to the ν s =(1-1.5)×10 −5 case the significance of the arrangement of the secondary's lobes is also apparent. In general, if the smaller lobe begins closer to the primary (ν s >0.5) then the secondary fission process is slower than when the larger lobe is closer to the primary (ν s <0.5). As the lobes grow closer in size the secondary fission process speeds up. This is because, as the smaller lobe grows relative to the larger lobe, the separation between their centers grows. Under the inverse square form of the lobes' mutual gravity this reduces the force that each lobe exerts on the other. The growth in mass of the secondary also increases the gravitation from the primary. The structure of the secondary, whether it be mostly one coherent mass or an amalgamation of several large masses, will have significant influence on the secondary fission process. For both models it is clear that in most cases secondary fission is a feasible evolutionary event for a low-mass-ratio binary.
Because the secondary fission analysis is done as a postprocessing step, the choice of model affects only the time of fission and not the system dynamics. As an example, Figure 7 identifies the complex spin state of the secondary at fission for the bi-lobed case when ν s =0.5. This case has a very short median time of secondary fission at roughly 1 day, yet still sees significant excitement of the spin state toward complex or nonprincipal axis rotation.
The cases of escaping secondaries can give us insight into a potential formation process for asteroid pairs. Here the relative velocity at escape and the spin states of each body at escape are of interest. Because the escape condition is the moment when the secondary crosses the primary's sphere of influence, we can provide a histogram of relative velocities when a potential asteroid pair could have formed (Figure 8).
We see that the relative velocities are fairly slow, of the order of 5-10 cm s −1 , and fairly well distributed across this range. This distribution provides a point of comparison for studies of asteroid pairs. Once again plotting the dynamic inertia and The x-axis provides both the effective spin rate and effective spin period, with the rate above and the period below in parentheses. The dynamic inertia is normalized by the mass-equivalent sphere inertia and the effective spin period is normalized by 2.33 hr, the disruption period of a 2 g cm −3 object.  Note.The median time for secondary fission is provided in parentheses in the last column. The cohesive model is denoted by COH and is described by Equation (1).
The bi-lobed sphere model is denoted by BL and the value of ν s . The entry labeled IBL identifies a secondary fission model that uses the bi-lobed framework but ignores the gravitational influence of the primary.
effective spin of the secondary, we can probe the behavior at the moment of escape (Figure 9). These spin states imply a fairly excited state for the ejected secondary that would likely have a significant effect on its continued evolution as a single and small complex rotating asteroid. An asteroid pair formed from binary fission would thus be expected to exhibit tumbling or non-principal axis rotation. The same cannot be said of the primary, whose spin state is largely unaffected by the ejection of the secondary, beyond the minor effects of the initial fission.
The small handful of captured secondaries remaining in orbit after a year of integration are likely statistically insignificant, but can still provide useful insight. Probing their secondary spin states after a year of integration reveals very excited systems likely not far from a disruptive event. Because these systems are highly excited and few in number, their spin states are widely distributed but without a clear pattern and are not shown in a figure.

1996 HW1
HW1 is used to represent two categories of NEA binary: fissioning contact binaries and high-mass-ratio binaries. The current geometry of HW1 is that of a contact binary, but we use a fictitiously fissioned model for this study (Figure 1). The current observed characteristics of HW1 are presented in Table 1 (Magri et al. 2011).
The fictitiously fissioned shape models were developed in past work by Hirabayashi & Scheeres (2019), which derived the minimum strength of the asteroid for its current geometry. While these shape models both have sharp angles at the neck where the fission is assumed to occur, the approximation of the mass distribution at the fourth order smooths this discontinuity because it does not converge identically to the shape. Table 5 identifies the uniform distribution of initial conditions used for the Monte Carlo simulations of the fission.
For the simulations of HW1, the larger secondary in such a system significantly slows the transfer of angular momentum from the primary to the secondary. This results in a much different evolutionary pathway. The statistical results from these simulations are presented in Table 6.
The most significant difference from the low-mass-ratio Moshup system is the lack of secondary fission and rarity of escape. We also see that the rate of collisions decreases slightly and the rate of chaotic capture increases significantly. While the rate of collision is likely dependent on the dynamical state, the rate of capture results from the slower angular momentum transfer, which prevents disruption of the system by escape or secondary fission. This agrees with the results of JS2011, whose longer integrations of high-mass-ratio binary evolution indicate a slower relaxation toward the outer stable doubly synchronous equilibrium of the order of 10 4 -10 6 yr.
While the collision rate for HW1 is relatively close to that for the Moshup system, the behavior is much different. In   In this case we see that all of the collisions occur immediately after the fission, rather than more slowly falling inwards as was the case for Moshup. This implies that the topology of the inner unstable equilibrium for HW1 is actually shifted further away from the body by higher-order effects. Likely, the initial location for many of these simulations is within this expanded collision barrier. Figure 11 shows an example of the evolution of the relative separation for a collision case.
For this case we see that the system is initialized just above the collision separation. At roughly 0.02 days the asteroids experience their initial collision. Within the post-processing this would be considered the time of collision. The evolution of the dynamics beyond this point shows a quick descent as opposed to a slow orbital decay over several periods.
Turning next to the distribution of the surface velocity at impact, we see similarly slow velocities ( Figure 12). The radial velocity range is much larger in magnitude than that of Moshup and does not come nearly as close to 0. This results from the more immediate collisions seen for this case, likely due to the warped topology of the barrier to collision.
Although there are only two escaping cases and they are likely outliers, it is still of interest to investigate their states at the moment of escape. We provide the complex spin state of the primary (Figure 13). We see that the spin state of the primary can be much more excited at the moment of escape than for the case of Moshup. The excited primary spin state results partially from its elongated shape being more susceptible to gravity torques from the secondary but is also due to the larger mass of the secondary. This may imply that for high-mass-ratio asteroid pairs the spin states will be tumbling or will exhibit non-principal axis rotation.
In addition to the particular details of the asteroid dynamics at the time of escape, HW1 provides a good example to explore the process of angular momentum transfer that excites binary orbits from captured to escape trajectories. This process has been studied extensively by Scheeres, Harris, Walsh, and The spin states are widely distributed, indicating a generally excited spin state of the secondary as it leaves the system. The x-axis provides both the effective spin rate and effective spin period, with the rate above and the period below in parentheses. The dynamic inertia is normalized by the mass-equivalent sphere inertia and the effective spin period is normalized by 2.33 hr, the disruption period of a 2 g cm −3 object. Note.We sample uniformly between the upper and lower perturbation bounds. The relative Euler angles define the relative orientation between the primary and secondary away from the current contact geometry. For the radial velocity and system spin rate two values for each entry are provided, where the first corresponds to the cohesion-less fission and the second corresponds to the fission with 4 hr spin period. The essential idea of the process is that as the secondary moves toward periapsis its orbital motion couples with the primary's spin, and the primary's shape is able to gravitationally torque the orbit of the secondary, increasing the angular momentum and energy stored in the mutual orbit. In Figure 14 we identify this interaction in a simulation of HW1 that occurs just before the secondary escapes the system. As the secondary approaches periapsis, in frames 1 and 2, it is moving with a lower angular velocity than the primary's spin. Its periapsis occurs nearly in exact alignment with the primary's major axis, maximizing the gravitational torque it experiences. After periapsis passage the secondary has been orbitally torqued such that its angular velocity surpasses that of the primary. Shortly afterwards the secondary is moving along a trajectory to escape the system.
The abundance of chaotically captured secondaries at the end of a year of integration is in stark contrast to the low-mass-ratio Moshup results. As previously described, the higher mass ratio of HW1 prevents the large flow of angular momentum from the primary to the secondary. Instead the bodies exchange angular momentum between each other and the orbital state more freely. While this prevents either body from reaching a critical spin rate for fission regardless of the fission model used, it does lead to a much more excited spin state of each body and a rapidly changing mutual orbit. Figure 15 provides the complex spin state for the chaotically captured primaries after a year of integration.
The excitement of the primary toward the state of nonprincipal axis rotation is extreme when compared to the lowmass-ratio cases. The secondary spin states show similar levels of excitement and are thus excluded. Because no secondary fission occurs for these cases, even when the simulations are extended for an additional year of integration, slower methods of energy dissipation dominate. Over these longer periods of time, the planar second-order dynamics modeled in JS2011 likely dominate and the systems would slowly relax toward the doubly synchronous state, as described in JS2011.
We further explore the effect of the initial fission spin rate and potential for cohesive fission to change the end state of the binary's evolution. To do this we run an additional set of simulations for HW1 with a fission spin rate of 4 hr, as opposed to the previous 4.74 hr. Apart from the spin rate, all other states are initialized with the same values as presented in Table 5 above. We do not recompute the unstable equilibrium in order to isolate the potential effects of cohesion under the same fission geometry. Table 7 provides the statistics for this case of faster spinning fission.
The increase in the initial energy and angular momentum at fission results in a much faster and nearly inevitable disruption of the system, with 93% of secondaries escaping in a median time of 6.6 days. All other systems collide with the primary. This suggests that, at least for high-mass-ratio binaries, the Figure 10. Location of impacts for collision cases shown in the body-fixed frame of the HW1 primary. The points are colored by the magnitude of the surface velocity at impact. All secondaries are initialized as roughly resting on the flattened surface of the primary shape model. Figure 11. Evolution of the relative separation in a colliding case for HW1 simulations. The collision separation (red) represents the minimum distance between the two asteroid barycenters before they collide given their instantaneous relative attitude and position throughout an integration. The relative separation (blue) shows the actual separation between the two asteroid barycenters throughout the integration. Once the relative separation falls below the collision separation the asteroids are considered to have collided. Here collision occurs just before 0.02 days of integration. initial fission must be relatively gentle and not dominated by cohesive forces.

2000 DP107
Previous studies of rotational asteroid fission and binary pairs predict equatorial craters to be a result of the inset mass ejection fission model and identify potential craters on binary pair shape models (Tardivel et al. 2018;Polishook & Aharonson 2020). Both NASA's OSIRIS-REx mission and JAXA's Hayabusa 2 identify this type of topography on their target asteroids, potentially suggesting that the fission process is more common in the NEA population than previously thought Watanabe et al. 2019). Radar shape models of DP107 identify an equatorial crater on the primary similar in size to the secondary asteroid . The current state based on these measurements is provided in Table 1. DP107 thus provides a convenient example to test the feasibility of the equatorial crater as a topographic marker for this model of fission by testing formation rates of binaries and asteroid pairs. To approximate this type of excavating fission we seed the unstable equilibrium solver with the secondary set within the current crater (Figure 1). Because the secondary is placed within the crater and not principally aligned, the equilibrium solver identifies a facsimile to the unstable equilibrium that is essentially an orbit with a temporary doubly synchronous state. The distribution of initial conditions sampled for the Monte Carlo analysis is provided in Table 8. DP107 was studied both to probe the low-mass-ratio binaries and to contrast the inset mass ejection formation with the Figure 12. Instantaneous surface velocities of colliding secondaries at the moment of impact for HW1 simulations. The velocities are provided in the secondary's orbit frame where the radial direction is from the center of mass of the primary toward that of the secondary, the cross-track direction is along the orbit in the direction of the orbital motion, and the normal direction is along the orbit normal. The surface velocity is the relative velocity between the surfaces of the primary and secondary at the moment of impact, accounting for relative spins and velocity. Figure 13. Complex spin state of the primary at the moment of secondary escape for HW1 simulations. The horizontal bars indicate the values of the secondary's moments of inertia. The x-axis provides both the effective spin rate and effective spin period, with the rate above and the period below in parentheses. The dynamic inertia is normalized by the mass-equivalent sphere inertia and the effective spin period is normalized by 2.33 hr, the disruption period of a 2 g cm −3 object. contact or resting model. Because the crater on the DP107 primary is not principally aligned, this example cannot be initiated identically at the unstable equilibrium, which complicates the comparison to Moshup. The statistics from the DP107 simulations are presented in Table 9.
Unfortunately, for the case of DP107 we see that all secondaries collide for the nominal Monte Carlo perturbations within a small fraction of a day. The likely cause of this is the placement of the secondary within the recess of the crater. At the point of fission the secondary does not have the outward velocity or energy to quickly move out of the crater without hitting the crater's edge, nor has enough time elapsed for the rotational state of the primary or secondary to become excited.
Looking at the impact locations (Figure 16), we can see that this explanation describes a large fraction of the collisions, but for a subset of cases the secondary is able to leave the confines of the crater, but it re-collides shortly after. These cases likely fission onto a highly elliptical or parabolic orbit and begin their descent back toward the primary before their orbits have been sufficiently circularized.
The distribution of the surface velocity directions (Figure 17) appears to further support this explanation. First, it shows that a majority of the collisions occur at low, and sometimes outward, Figure 14. Example of angular momentum transfer in the HW1 system. The six frames show the secondary moving through periapsis and gaining angular momentum from the primary's spin and shape before moving toward an escaping trajectory. The system is shown in the primary's principal frame with its angular momentum oriented out of the page. Figure 15. Complex spin state of the primary in chaotically captured cases after one year of integration for HW1 simulations. The horizontal bars indicate the values of the secondary's moments of inertia. The spin states are broadly distributed throughout the dynamical space, with a slight skew toward the intermediate moment of inertia. The x-axis provides both the effective spin rate and effective spin period, with the rate above and the period below in parentheses. The dynamic inertia is normalized by the mass-equivalent sphere inertia and the effective spin period is normalized by 2.33 hr, the disruption period of a 2 g cm −3 object. radial velocities. These points correspond to the secondaries that fail to escape the crater, colliding within the crater surface or edge. These cases impact mostly with cross-track and normal velocity on their way out of the crater. Comparing the impacts with larger radial velocity, we see that the velocity magnitude in the other directions is much smaller, implying a more direct collision than what was seen for Moshup or HW1. An argument can be made for the DP107 cases that the secondary may have sufficient energy to fission and launch one portion of the body into the orbit under the bi-lobed secondary fission model. This argument can be taken to the extreme by applying the model with equally sized lobes, equivalent to a mass fraction of 0.5. Unfortunately, when this is explored (Table 10) the majority of the re-collisions occur so quickly that only a single secondary experiences fission before returning to the surface of the primary.
Further exploring the results of the bi-lobed secondary fission model for DP107, we can identify the bounds of the mass fraction of the secondary that enable secondary fission to occur before re-collision. The identified condition is n   0.107 0.816, s telling us that, in addition to the rarity of secondary fission occurring quickly enough, it requires that the secondary be composed of elements with relatively large mass. Based on the results of this initial set of simulations it is clear that for the DP107 secondary to successfully fission from the equatorial crater, the initial conditions must be more energetic. This would allow the secondary to escape the crater and enter a captured orbit. We use two approaches to try to explore the conditions for DP107 formation. The first is to assume that cohesion plays a more significant role in the fission, as was attempted with HW1. In this case the Monte Carlo conditions and perturbation bounds are maintained, but the system's initial spin period is sped up from 5.08 to 3.85 hr. The second approach is to explore the dynamical space more widely by increasing the perturbation bounds while maintaining the same initial conditions. The perturbation bounds for this wider Monte Carlo analysis are provided in Table 11.
The results of these Monte Carlo sets are best compared using both the cohesive and bi-lobed secondary fission models. Note.We sample uniformly between the upper and lower perturbation bounds. The relative Euler angles define the relative orientation between the primary and secondary away from the current contact geometry.  Figure 16. Location of impacts for collision cases shown in the body-fixed frame of the DP107 primary. The points are colored by the magnitude of the surface velocity at impact. All secondaries are initialized within the crater seen on the right side of the equator in this view.
In Table 12 the statistical results for both Monte Carlo sets post-processed with both secondary fission models are provided. The broad conclusion to be drawn from these results is that in both cases the majority of secondaries re-collide or escape unless the secondary is composed of elements with fairly large mass. In the case that the secondary does fission it occurs very quickly, potentially in less than an hour. This suggests that formation of binaries with equatorial craters such as DP107, assuming the secondary fission pathway described in JS2011, is only likely given a particular case of a rubble pile asteroid with sufficiently large boulders or other massive elements. Looking more closely at the case of a 3.85 hr initial spin period, the dynamics appear to dictate that similar numbers of secondaries either collide or escape, ignoring secondary fission of these cases. There is a delicate balance based on the initial spin period that leads to this; if the initial spin period is reduced to 4 hr the secondaries all re-collide, and if the spin period is increased to 3.75 hr the secondaries all escape. This presents an interesting comparison with the 2.2 hr spin barrier traditionally observed in the NEA population, suggesting that systems fissioning near this barrier would be more likely to form asteroid pairs than binary systems. In the case of wider Monte Carlo perturbation bounds the vast majority of secondaries recollide, with a small subset having the potential to escape due to the relatively large range of velocity perturbations. One possible explanation for the difficulty in identifying a formation pathway for DP107 is the initialization of the secondary near the crater as opposed to it being closer to the unstable equilibrium. This would reinforce the importance of this state as a boundary or gateway for binary formation. As such it may have implications for the mass distribution of binaries with equatorial craters, assuming the crater was the original site of the secondary before fission.

Discussion
Many past studies of binary formation and evolution have relied on a range of assumptions in order to broadly study the system dynamics or particular processes. Within this work we are able to apply computationally efficient high-fidelity dynamics models in order to lift these assumptions and more rigorously test the theories and conclusions presented in past work. In broad strokes, we reconfirm many past findings. We find that low-mass-ratio binaries such as Moshup and DP107 tend toward disruption and require a fast acting means of energy dissipation, such as secondary fission, in order to relax into their current observed orbital states. Likewise, the slow angular momentum transfer available to high-mass-ratio binaries prevents excitement to secondary fission, instead suggesting a slower evolutionary pathway as described in JS2011. We additionally find explicit examples of the angular momentum transfer from an elongated primary into the mutual orbit of a binary leading to escape in the HW1 simulations. Similar behaviors occur for Moshup and DP107 although they are not as easily illustrated (Scheeres 2001). Asteroid pairs formed from escaping secondaries by this process are also expected to have tumbling spin states, at least after their initial separation. In studying the formation of binaries via inset mass ejection we identify the difficulty of this process without the inclusion of secondary fission. Likely this can be argued to be a result of the dynamical structure of the F2BP away from the unstable equilibrium. In general, we show that the inclusion of higher-order gravity terms and nonplanar dynamics slows the evolutionary processes at play but does not significantly alter the formation of binaries. The more complex dynamics do topologically morph the zero-velocity curves near Figure 17. Instantaneous surface velocities of colliding secondaries at the moment of impact. The velocities are provided in the secondary's orbit frame where the radial direction is from the center of mass of the primary to that of the secondary, the cross-track direction is along the orbit in the direction of the orbital motion, and the normal direction is along the orbit normal. The surface velocity is the relative velocity between the surfaces of the primary and secondary at the moment of impact, accounting for relative spins and velocity. the unstable equilibrium, introducing the possibility of lowvelocity re-collision with the primary. Given the relatively high rate of re-collision we see, it would be expected that binaries or contact binaries may have small craters or impacts from these events. There is also a possibility that under precise conditions a secondary could bounce after re-collision, dissipating energy and being captured into orbit. This could be of particular interest for the case of DP107 and equatorial crater formation. A better understanding of these interactions would be necessary to identify the exact observable markers. In detailed comparison with JS2011 we identify slower and more fraught formation processes. For a low-mass-ratio binary, such as Moshup, we see that both escape and secondary fission occur on a much slower timescale, albeit at the same rate. As a result of the re-collision events possible in our model, we also see a significantly decreased rate of escaping secondaries. The rates of chaotic capture and secondary fission remain roughly the same. This suggests potentially fewer asteroid pairs resulting from this type of binary fission than may have been expected in JS2011. For the high-mass binaries, such as HW1, JS2011 shows that most secondaries remain chaotically captured and relax slowly toward the doubly synchronous equilibrium. Our study of HW1 supports this with the caveat of potential re-collision. While our study of DP107 deviates from JS2011 it does support the need for secondary fission as a means to dissipate the energy of the secondary in order to form a stable low-mass-ratio binary.
In our study of secondary fission two models were tested against one another, one assuming only cohesive forces while the other focused on the gravitational effects of the primary on the components making up the secondary. In the case of the Moshup system we show that, regardless of the secondary fission model employed and the size of the fissioning particles, between one third and two thirds of systems experience secondary fission. While the ejection of meter-sized objects from the secondary is not likely to alter its orbit as significantly as larger masses, realistic secondaries are likely to be more complex than the bi-lobed model used. This suggests that secondary fission remains a likely pathway for the evolution of low-mass-ratio binaries, while also indicating that a more complex and rigorous study be performed. The difficulty of stabilizing a binary such as DP107 before re-collision without secondary fission further supports this argument, although the system is a unique subset of binary fission. In comparison with previous work by Walsh et al. (2012), where the primary and secondary asteroid are made up of objects near 150 m, our secondary fission study of Moshup (Table 4) seems to suggest many of these binaries are likely to experience secondary fission.

Conclusions
Building on past work by Scheeres, Jacobson, Walsh, and others, we are able to apply improved simulation tools for binary asteroid to further analyze the formation and evolution of these systems. Moshup, 1996 HW1, and 2000 DP107 are used as example systems to capture the particular evolutionary pathways identified in past studies. For each of these systems a Monte Carlo analysis of their conditions at formation is performed to understand the dynamics at play. Our analysis identifies many new questions to be addressed. Re-collision brings up two topics to be pursued. First, could these recollision events leave characteristic craters that may be identifiable through close observation and what constraints could they provide? Second, it highlights the need for a more thorough dynamical analysis of the structure at the unstable equilibrium. What are the exact sources of the structural warping that allows for re-collision and can the warped dynamical structure be mapped in a generic way to provide broader insight? Study of re-collisions that model the continuum dynamics and energy dissipation due to these events could be used to explore a bouncing secondary as a pathway to binary formation. While our comparison of secondary fission does begin to shed more light onto these Note.We sample uniformly between the upper and lower perturbation bounds. The relative Euler angles define the relative orientation between the primary and secondary away from the current contact geometry. Note.The median time for secondary fission is provided in parentheses in the last column. The cohesive model is denoted by COH and is described in Equation (1). The bi-lobed sphere model is denoted by BL and the value of ν s . processes, long-term simulation of fissioned particles or lobes is important, particularly study of the kind of evolutionary pathways that exist for debris from the secondary and how quickly the models of secondary fission dissipate energy. That analysis will require the implementation of a broader dynamics model of the full N-body problem. Finally, we highlight the tendency toward non-principal axis rotation of fissioned asteroids. Further study would have implications for asteroid pairs and asteroids with unexplained complex rotation states.
This material is based upon work supported by the National Science Foundation Graduate Research Fellowship Program under grant No. DGE 1650115. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

Appendix A The Hou Mutual Gravity Potential
The general-order reformulation of the mutual gravity potential by Hou et al. begins from the double integral description of the mutual gravity potential:   The variables R and e x , e y , and e z represent the magnitude and unit direction of the relative separation. The variables T A and ¢ T B are the mass-normalized inertia integral sets for the primary and secondary bodies, where the prime denotes that the inertia integrals of the secondary are rotated into the frame of the primary.

A.1. Inertia Integrals
Central to this reformulation of the mutual potential is the use of inertia integrals to describe the mass distribution. This aspect of the reformulation is accomplished by the application of a Legendre polynomial expansion to describe the mass distributions, where the Legendre coefficients are referred to as inertia integrals: In this way the inertia integrals can be considered analogous in use to spherical harmonics (Tricarico 2008). The mathematical form of the inertia integrals is similar to that of the moments and products of inertia for a rigid body wherein each term represents the mass distribution about some axis; however, the inertia integrals are expanded to order N whereas moments of inertia are linear combinations of second-order inertia integrals.
Here we provide the mass-and length-normalized form of an inertia integral along with the zeroth-and second-order coefficients in terms of the normalized moments and product of inertia:  T  I  T  T  I  T  T   I  T  T  I  T   I  T  I  It is of note that the zeroth-order inertia integral is equal to the mass for the non-normalized form, thus it is equal to one in the normalized form.

Appendix B General Use Binary Asteroid Simulator
As a part of the analysis for this paper we developed a tool for dynamical propagation of binary asteroids of arbitrary shape and expansion order using the Hou mutual gravity potential described in Appendix A (Hou et al. 2017). We have provided the software tool for free use at https://github.com/ alex-b-davis/gubas. The tool, referred to as the General Use Binary Asteroid Simulator, is intended to provide the planetary science community with an easily used, fast, and high-fidelity simulation tool for the numerical integration of binary asteroid dynamics. It does not include the tools for the fundamental frequency analysis performed in this paper.
The software was designed and implemented to be highly modular, and thus to enable a wide set of uses and allow for easy integration into larger tool sets. For this reason the architecture was centered around a C++ executable wrapped in a Python shell. The C++ executable performs the numerical integration and calculation of the inertia integrals while the Python wrapper pre-processes user input from a configuration file to initialize the executable and post-processes the results. This approach allows the user to easily modify the Python shell script to fit their needs. In the standard architecture all interactions are handled through the configuration file and a single command line call to initialize the process. While the software and a detailed user guide can be found by following the github link, Figure B1 shows a basic flowchart of the software process.