© A. Ho et al. 2022

1 Introduction

Since the first binary asteroid system, (243) Ida and its moon Dactyl, were discovered by the Galileo spacecraft (Chapman et al. 1995), many more have been identified among near-Earth objects, in the main belt, and in the Kuiper belt (see, e.g., Margot et al. 2015, and references therein). Roughly 27 000 near-Earth asteroids (NEAs) are currently known, the majority of them with diameters less than 1 km (Harris & Chodas 2021). NEAs are thought to originate from the main belt and, due to resonances with Jupiter, to have migrated into Earth-crossing orbits with perihelion distances of < 1.3 AU (Morbidelli et al. 2002). It is estimated that roughly 16% of near-Earth objects are binaries (Margot et al. 2002).

It is believed that smaller binary systems among asteroids are formed through rotational fission (Margot et al. 2002; Pravec & Harris 2007). Small asteroids, typically with diameters ~0.1–10 km (Walsh 2018), are “rubble piles”, porous collections of irregularly shaped boulders and finer grains held together by gravity and possibly weak cohesion forces (Hirabayashi et al. 2015; Li & Scheeres 2021). In the rotational fission model a rubble pile asteroid is spun up by the Yarkovsky-O‘Keefe-Radzievskii-Paddack (YORP) effect (Rubincam 2000). Once the asteroid reaches a critical spin rate, it will start to shed some of its mass (Scheeres 2007; Walsh et al. 2008). This model also matches the observations of rapidly rotating primaries of asteroid pairs (Pravec et al. 2010, 2019).

Other binary creation processes have also been proposed, such as binary creation by collisions and even creation via tidal disruptions from nearby planets (see, e.g., Margot et al. 2002; Merline et al. 2002; Richardson & Walsh 2006). The first mechanism is likely to describe formation of binaries of large asteroid systems (see, e.g., Walsh & Jacobson 2015). However, it is believed that creation of binaries among the NEA population is highly unlikely through these mechanisms.

Various works have studied the dynamics of an asteroid binary system during and after the fission process. Walsh et al. (2008) modeled asteroids as rubble piles consisting of numerous self-gravitating spheres. In their model, the YORP spin-up would eject some of these spheres, and they found that the formation of a satellite was more efficient for a spherical and oblate primary. The work of Scheeres (2007) considered a slightly different scenario in which the asteroids are initially resting on each other; these are known as contact binaries. Scheeres studied the limits in which fission would take place, considering an ellipsoid-sphere model and extended this to a two-ellipsoid model to study the stability of the binary system post-fission (Scheeres 2009). However, the systems predicted by these theories are highly energetically excited. In order to stabilize the systems and prevent the secondary from escaping, a form of energy dissipation mechanism is necessary.

Work by Jacobson & Scheeres (2011) studied the creation of various NEA binary systems, including doubly synchronous binaries, high-e binaries, ternary systems and contact binaries. They introduced a new binary process, secondary fission, as a mechanism to decrease the energy level of the system. This was extended by Boldrin et al. (2016) to include non-planar effects, and they found that secondary fission can take place at higher mass ratios, compared to Jacobson & Scheeres (2011), as a nonplanar configuration allows for higher energy levels. They also found that the secondary acquired nonprincipal axis rotations as a consequence of the nonplanar effects.

Davis & Scheeres (2020) further studied post-fission dynamics by including higher order gravity terms in addition to nonplanar effects, and also included tidal torques. Davis and Scheeres compared their results with Jacobson & Scheeres (2011), and found that the formation processes remain unaltered, but that the process itself is slower. Additionally, due to the possibility of re-collision in their model, they found that the rate of escaping secondaries is lower.

In this paper, we study the dynamical evolution of asteroid binary systems immediately after fission occurs. Our work is similar to the work of Boldrin et al. (2016), who assumed rotational fission of a contact binary where the secondary was given different initial tilt angles about its body-fixed y-axis. We investigate the outcome of the rotational fission for a number of different mass ratios and configurations of the contact binary. Whereas the Boldrin et al. study was restricted to systems with mass ratios q ≤ 0.3 where the density and shape of the secondary was identical to that of the primary, we include the whole range of mass ratios from 0.01 to 1, and also allow for different densities and shapes of the secondary. Our work applies the method we recently developed that computes the forces and mutual torques between two bodies without using approximations (Wold & Conway 2021; Ho et al. 2021). When expanding the mutual potential, for instance through spherical harmonics, higher order terms have a more significant role in the dynamics of the system when the bodies are closer. Furthermore, Hou et al. (2017) showed that higher order terms are required when the bodies are also more elongated. This means that using an exact method may provide more accurate results of the dynamics of asteroid binaries or pairs after the initial separation.

The structure of this paper is as follows. In Sect. 2, we present the mathematical framework and initial conditions used for our models. In Sect. 3, we describe the models used and present the results of the simulations. Finally, we summarize and discuss our results in Sect. 4.

2 Dynamical model

The model consists of two triaxial ellipsoids. Initially they are attached to each other as a contact binary. We assume that the contact binary undergoes rotational fission, a process where the two components separate when a certain limiting rotational speed is reached (Bottke et al. 2002; Scheeres 2007; Walsh et al. 2008). The initial setup is shown in Fig. 1, and is similar to that used by Boldrin et al. (2016), with the secondary centered on the long semiaxis of the primary and rotated an angle θ₀ about its body-fixed y-axis.

To compute the force and torque on body i in the gravitational field of body j, we apply the surface integral equations described by Conway (2016):

$$F_{\text{G}}=G{\rho }_i{\displaystyle \underset{S_i}{\iint }{\text{Φ}}_j\left(r\right)n\text{d}S},$$(1)

$$M=-G{\rho }_i{\displaystyle \underset{S_i}{\iint }{\text{Φ}}_j\left(r\right)n\times r\text{d}S},$$(2)

where the mutual potential between the two bodies is written as:

$$U=\frac{G{\rho }_i}{3}{\displaystyle \underset{S_i}{\iint }\left(r{\text{Φ}}_j\left(r\right)-\frac{1}{2}{\left|r\right|}^2{g}_j\left(r\right)\right)}·n\text{d}S.$$(3)

In these formulae ρi is the density of body i (assumed to be constant throughout the body), and Φ_j(r) and g_j(r) = ∇Φ_j are the scalar potential and gravitational field of body j at a position r on the surface of body i. The vector normal to the surface of body i at position r is n, and dS is the surface element at that position. The gravitational constant is denoted as G.

It is customary to use second- or fourth-order approximations of the mutual gravitational potential for two-body interactions of nonspherical bodies, and from that compute force and torque (Fahnestock & Scheeres 2008; Boldrin et al. 2016; Hou et al. 2017; Davis & Scheeres 2020). The mutual potential is thus expressed as a sum of several terms, which in fact suffers from a truncation error. However, our approach uses exact expressions in the form of surface integrals, and will therefore not suffer from truncation errors. For ellipsoids the potential of the bodies (Φ) can be expressed using well-known analytical expressions (MacMillan 1930). The surface integration scheme thus becomes a surface integration over an ellipsoid surface (see Wold & Conway 2021, for a more detailed outline of the surface integration).

We propagate the binary after rotational fission by solving the rotational and translational equations of motion in an inertial frame of reference, formulating it as a standard initial value problem. The rotational motion of the bodies is solved in the body-fixed reference frames using Euler parameters (e₀, e₁, e₂, e₃) in order to avoid the singularities related to Euler angles. For the integration of the equations of motion, we use the ninth-order Runge-Kutta method by Verner (2010). While it is convenient to use an adaptive time stepper, we use the solver with a fixed time step of Δt = 19 minutes in order to compare the time evolutions between various simulations. Furthermore, we do not make use of an adaptive time stepper because our simulations are relatively short. The end results did not have significant changes when the time step was smaller, nor did an adaptive time stepper affect the outcome.

Fig. 1 Contact binary before separation (top) and after separation (bottom). Top: Configuration of the contact binary the moment of fission where the gravitational force FG and the centrifugal force FC are equal. The cross indicates the center of mass of the system. The long and short semiaxes (a and c) are aligned with the body-fixed x- and z-axes of the respective bodies. Bottom: Angular velocities of the bodies after the contact binary separates.

Fig. 2 How the secondary moves closer to the primary when the angle θ0 of the secondary is increased.

2.1 Rotational fission

Throughout the rest of this paper, all variables with subscript p and s correspond to variables describing the primary and secondary, respectively.

Initially, before separation, the contact binary rotates about an axis passing through the center of mass of the system and perpendicular to the xy-plane of the primary, as shown in Fig. 1. When the rotational speed reaches a certain limit ω₀, the centrifugal force on the secondary matches the gravitational attraction between the primary and secondary, and the contact binary fissions.

The initial angular velocity ω₀, which we use to start our simulations, is therefore the limit for rotational fission given by:

$${\omega }_0=\beta \sqrt{\frac{F_{\text{G}}}{m_{\text{s}}r}},$$(4)

where m_s is the mass of the secondary and r is the distance between the centroid of the secondary and the center of mass of the system (see Fig. 1). We found during our simulations that it was necessary to assume a value of ω₀ that is slightly higher than the theoretical limit, hence we multiplied the theoretical limit by the factor β = 1.01. The β-factor can be interpreted as some cohesion between the two components, and small amounts of cohesion may allow rubble pile asteroids to rotate faster than the theoretical limit (Holsapple 2007; Sanchez & Scheeres 2014).

2.2 Initial conditions

As the system is not affected by external forces or torques, linear and angular momentum is conserved. Furthermore, no energy is added or removed at the instant of fission. Thus, immediately after fission the primary and the secondary both experience the same angular velocity ω₀. Therefore, the initial translational velocities of these two objects right after fission can be found as:

$$v_{0\text{,p}}={\omega }_0\times \left(r_{\text{0,p}}-r_{\text{cm}}\right),$$(5)

$$v_{0\text{,s}}={\omega }_0\times \left(r_{\text{0,s}}-r_{\text{cm}}\right),$$(6)

where r_0,p and r_0,s are respectively the initial positions of the primary and secondary in the inertial frame, r_cm is the position of the center of mass of the system, and ω₀ = [0,0, ω₀] is the initial angular velocity vector in the center of mass system. After the bodies have separated, the angular velocities of the bodies in the inertial frame are equal to those of the contact binary before separation, as shown below the dashed line in Fig. 1. The angular velocities in the body-fixed frames are determined as:

$${\widehat{\omega }}_0=R^{\text{T}}{\omega }_0,$$(7)

where R^T is the transpose of the rotation matrix at the time of separation.

The configuration is varied by changing simultaneously the initial angle θ₀ of the secondary and the centroid-to-centroid distance between the primary and secondary, under the condition that the separation between the two surfaces at their closest point is kept at ∆r = 1 cm. When θ₀ = 0°, the initial positions of the primary and the secondary are:

$$r_{0\text{,p}}=\left[0,0,0\right],$$(8)

$$r_{0\text{,s}}=\left[a_{\text{s}}+a_{\text{p}}+\text{Δ}r\text{,}0,0\right],$$(9)

where a_p and a_s are respectively the long semiaxes of the primary and secondary and ∆r = 1 cm is the separation between the surfaces.

When θ₀ increases from 0 to 90°, the surface-to-surface distance increases. In order to keep this distance at 1 cm, the secondary’s centroid has to be moved closer to the primary’s centroid, as illustrated in Fig. 2. In this manner, we ensure that the initial separation between the surfaces is always 1 cm. In practice, when θ₀ changes, the initial position of the secondary r_0,s is recalculated by a separate algorithm.

By keeping the initial distance between the surfaces to 1 cm regardless of the value of θ₀, the limiting value of ω₀ for the initial fission will increase. This is a consequence of r becoming smaller in Eq. (4). The variation in ω₀ with θ₀ is shown in Fig. 3. The top panel shows that ω₀ increases as a function of θ₀ when ∆r = 1 cm (blue crosses). However, when the centroid-to-centroid distance is kept constant, which leads to an increasing gap between the surfaces, the value of ω₀ decreases slightly as a function of θ0 (red crosses). Our model therefore takes into account that the limiting rotational speed for fission changes as the tilt angle of the secondary changes. The bottom panel shows the relative difference between these two cases for three different mass ratios of the primary to the secondary. The relative difference amounts to ≈ 15–20% when θ₀ approaches 90°. We also note that the relative difference grows larger as the mass ratio increases.

Throughout all the simulations the shape and density of the primary are fixed. The semiaxes are (a_p, b_p, c_p) = (1.0,0.7,0.65) km, equal to the numbers used by Boldrin et al. (2016), and the density is ρ_p = 2.0 g cm⁻³, which is a commonly used density to model rubble pile asteroids (Pravec et al. 2010; Jacobson & Scheeres 2011; Boldrin et al. 2016). Some observed asteroids also have densities close to this value, for example 25143 Itokawa (Fujiwara et al. 2006; Kanamaru et al. 2019), as do some primaries of asteroid binaries, such as (66391) 1994 KW4 (Moshup) (Ostro et al. 2006; Scheirich et al. 2021) and (88710) 2001 SL9 (Scheirich et al. 2021).

For each configuration defined by sets of θ₀ (and consequently r_0,s), our aim is to study how the dynamics of the binary system evolve while varying the mass ratio q = m_s/m_v. We ran simulations for 37 different mass ratios q = 0.01,0.05,0.10, q ∈ [0.15,0.30] in increments of 0.01 and q ∈ [0.32,1.00] in increments of 0.04. For each mass ratio q, we considered 45 different initial angles θ₀ of the secondary in the range θ₀ ∈ [0.001°, 90°]. All simulations were run with a time span of 4800 h (200 days), unless they were terminated earlier due to collision (or impact) between the two bodies.

Fig. 3 How angular velocity changes with the initial tilt angle θ0 of the secondary. Top panel: ω0 is affected when there is contact between the surfaces (blue) and when the surface-to-surface distance increases (red) as the secondary is tilted (see also Fig. 2). The mass ratio is q = 0.1 for the data in the top panel. Bottom panel: relative difference in ω0 between the two approaches shown in Fig. 2, but for three different mass ratios q. The mass ratio is defined as q = ms/mv.

3 Results

We examined the dynamics as a function of q and the initial tilt angle θ₀. The mass ratio can be written as:

$$q=\frac{{\rho }_{\text{s}}}{{\rho }_{\text{p}}}\frac{a_{\text{s}}{b}_{\text{s}}{c}_{\text{s}}}{a_{\text{p}}{b}_{\text{p}}{c}_{\text{p}}}.$$(10)

Because we keep the shape and density of the primary fixed, varying the mass ratio of the system mainly affects the mass and volume of the secondary. Moreover, increasing the mass ratio also changes the total energy of the system, as shown in Fig. 4. The total energy is the sum of kinetic and potential energy. The systems where the total energy is negative are bound; in systems where the total energy is positive the two components can undergo mutual escape.

First, we consider models with three fixed values of ρ_s/ρ_p while keeping the ratio of the secondary’s semiaxes equal to that of the primary. In the next three models the secondary can take different geometrical shapes, but now the density is kept constant and equal to that of the primary.

In order to determine whether the secondary has escaped or exists in an unstable orbit, we utilize its orbital eccentricity e. The eccentricity is an osculating Keplerian element, and therefore changes with time. The secondary is considered to have escaped when e ≥ 1 for at least 50 time steps. This is to ensure that cases where e ≥ 1 for only a shorter period of time are not classified as already escaped. Increasing this limit to more than 50 time steps did not change the outcome. If, however, the eccentricity is less than unity at the end of the simulation, and the total energy of the system is positive, we classify it as residing in an unstable orbit. The secondary in systems with negative total energies is classified as being in a stable orbit. If the ellipsoid surfaces intersect at any time during the simulation, we consider it a collision and end the simulation.

These definitions share some similarities with the definitions provided by Scheeres (2002). For instance, the outcome “eventual escape” outlined by Scheeres, where there are multiple periapsis passages that will eventually terminate, is similar to our definition of an unstable case scenario. The “nonim-pacting and nonescaping” outcome is equivalent to our stable orbit outcome. However, we do not classify immediate escape scenarios, nor we do distinguish between different reimpact events.

Fig. 4 Total energy of the system as a function of the mass ratio q, for each model, using θ0 = 0.001°. The gray dotted line shows the zero energy line.

3.1 Varied densities, models D1-D3

The first set of models considered involves varying the density of the secondary while keeping the semiaxis ratio equal to that of the primary, that is a_s/b_s = a_p/b_p and b_s/c_s = b_p/c_p. In this case the semiaxes of the secondary can be derived from Eq. (10), with a_s written as:

$$a_{\text{s}}={\left(\frac{{\rho }_{\text{p}}}{{\rho }_{\text{s}}}q\right)}^{1/3}{a}_{\text{p}}.$$(11)

The equations of b_s and c_s take similar forms.

We examine models with three different density ratios:

• Model D1:ρ_s/ρ_p = 0.5;

• Model D2:ρ_s/ρ_p = 1.0;

• Model D3:ρ_s/ρ_p = 2.0.

As the density of the primary is fixed at ρ_p = 2.0 g cm⁻³, models D1, D2, and D3 have secondaries with densities of 1.0 g cm⁻³, 2.0 g cm⁻³, and 4.0 g cm⁻³, respectively. Model D2 is identical to the model discussed by Boldrin et al. (2016).

Fig. 5 Model outcomes as functions of mass ratio q and initial tilt angle θ0. In the top row, the left, middle, and right panels correspond to models D1, D2, and D3, respectively. In the bottom row, the left, middle, and right panels correspond to models S1, S2, and S3. The dashed lines indicate the value of q where the total energy is zero.

3.2 Varied shapes, models S1–S3

In these models, we investigate cases where we vary the axis ratio of the secondary, but keep the density of the secondary equal to that of the primary. We write the secondary’s semiaxis ratios as:

$$\frac{a_{\text{s}}}{b_{\text{s}}}=f_1,$$(12)

$$\frac{b_{\text{s}}}{c_{\text{s}}}=f_2,$$(13)

We select three combinations of f₁ and f₂:

• Model S₁: f₁ = 1.3, f₂ = 1.03, a secondary that is fairly spherical and almost an oblate spheroid;

• Model S2: f₁ = 1.6, f₂ = 1.2, a cigar-shaped secondary with a ≫ b > c;

• Model S3: f₁ = 2.5, f₂ = 1.2, similar to model S2, but even more elongated.

3.3 Outcome distribution

First we study the outcome of the secondary at the end of the simulations. Figure 5 shows the distribution of the outcomes as functions of mass ratio and initial tilt angle for the six different models. In all six models, most of the simulations (more than 80% of the case results) end up with the two bodies colliding (red area in the figure). These collision events are typically found when θ₀ ≳ 15°.

In general, there are two regions (θ₀ ≲ 15° and θ₀ ≳ 80°) where the components do not impact, but where the secondary either escapes or orbits the primary. Most of these cases are found for configurations with θ₀ < 10–15° over the entire range of q. Those found at higher initial angles mainly take place at low mass ratios, and the number of them residing in this region is low for most models. These two ranges of θ₀ correspond to regions near two configurations (θ₀ = 0° and θ₀ = 90°) where the contact binary is in a relative equilibrium (Scheeres 2009).

The separation between the positive and negative total energy regimes in Fig. 5 occurs between the yellow and green areas. For two spheres this separation occurs at q = 0.2, and for triaxial ellipsoids, as in our case, the separateion fluctuates around this value depending on both shape and configuration (see discussion in Scheeres 2009 and Jacobson & Scheeres 2011). We find that the separation occurs at q = 0.19–0.20, q = 0.22–0.24, and q = 0.26–0.32 for models D1, D2, and D3, respectively. Hence the separation occurs at successively higher mass ratios when the density of the secondary increases. The separation shifts toward slightly higher mass ratios when θ₀ increases, as seen in the top regions of the panels. This occurs because the total energy is raised for these configurations, and also reflects an increased value of ω₀, as illustrated in Fig. 3. A similar trend is seen in the varied shape models, where the separation between positive and negative energy regimes occurs at lower mass ratios when the secondary becomes more elongated.

At low mass ratios where the total system energy is positive, we find a mix of cases where the secondary has escaped and where it is still orbiting the primary in an unstable orbit. With a longer simulation time, we expect to see fewer cases of secondaries in unstable orbits. Boldrin et al. (2016) call these secondaries “escape survivors”, and only find them at q > 0.27 (they only include systems with q < 0.3) after a simulation time of 200 yrs. As our simulations are 200 days long, they represent a snapshot of the situation after a fraction of this time. We therefore have survivor cases also at the lowest mass ratios, as opposed to Boldrin et al. (2016).

Of all the cases at θ₀ ≳ 80° that do not collide, a relatively large fraction have escaped compared to those at lower θ₀, typically making up more than 70% for most models. This indicates that the secondary may escape earlier the more tilted its initial position is. A higher initial angle corresponds to a higher energy configuration of the system, and may therefore be the cause of an earlier escape for the high initial angle systems.

At higher mass ratios q ≳ 0.2 the system total energy is negative, hence all secondaries are gravitationally bound to the primary (unless sufficient energy is added to the system). These are shown in green in Fig. 5. Some also appear when θ₀ ≳ 85°.

3.4 Collisions

The majority of the simulations end up with a collision. The collisions tend to occur for configurations with θ₀ ≳ 15°, but can also take place at lower initial angles when the density ratio ρ_s/ρ_v is lower, when the secondary becomes more elongated, or when the mass ratio increases. Allowing the secondary to become less dense or more elongated also increases the overall number of collision cases.

Although we find that collisions typically happen at θ₀ ≳ 15° (with the exception of model S3, where collisions can happen as low as θ₀ ≳ 8°), Boldrin et al. (2016) report that in their simulations collisions occur for initial tilt angles of θ₀ ≳ 40°. This cannot be due to our study having a shorter simulation time, as we would expect the opposite to happen if that were the case (as we expect more systems to collide over time). The most likely explanation is that the secondary in our study starts out closer to the primary when it is rotated (see Fig. 2). By moving the secondary closer, the probability of collision is also expected to increase. This is especially true when θ₀ is nonzero, as the secondary will fall onto the surface of the primary due to the gravitational torque. This also explains why there are significantly more collisions for the more elongated secondaries, as the gravitational torque is stronger when the secondary becomes more elongated. The simulations that survive at high angles are likely due to higher initial velocities as a result of higher system energies, which thus prevent this type of collision.

There is a sharp horizontal division separating collision and stable cases when θ₀ ~ 10°–15°. However, this is not found at higher angles. This may be because when θ₀ approaches 90 degrees, the secondary approaches an unstable equilibrium, whereas a lower initial angle is closer to a stable equilibrium.

Most of the collisions take place very early in the simulations. More than 95% of the impact events occur within the first five hours. Some of these impacts can occur even within the first two time steps, which make up 82% of the collision outcomes. The collisions that occur between the first and second time step may be considered the “immediate reimpact” events that are mentioned by Scheeres (2002). These early impacts are due to the secondary falling onto the primary as a result of the gravitational torque.

Finally, we study the remaining collision cases that occur later one in the simulation, at t > 5 h. These are shown in Fig. 6, distributed as functions of both q and θ₀. The top panels show that cases that survive longest, in all six models, have intermediate values of the mass ratio, typically between 0.18 and 0.4.

Compared to models D1 and D2, there is a tendency for model D3 to survive longer at both lower and higher mass ratios than this range. For instance, there are a few cases with q ≈ 0.5 and q ≈ 0.1 with a survival time ≳500 h which is not found in models D1 and D2. The collision time of model S1 is, on average, greater than those in models S2 and S3. In model S2, there are only two simulations that experience collision after 500 h, and only one in model S3, that occur when q ≈ 0.25. The bottom panels in Fig. 6 show that nearly all collisions that take place after 5 h have elapsed have secondaries with large initial tilt angles θ₀ > 80°. The one exception is for model D1, where the time before collision is approximately 59 h for a case with θ₀ ≈ 8° and q = 0.3 (corresponding to the “dent” in the green region in the top left panel of Fig. 5). On average, the time before collision, for simulations that last longer than 5h, is 133 h, 166 h, and 143 h for models D1, D2, and D3, respectively, while for the varied shape models the averages are 170 h, 192 h, and 79 h for models S1, S2, and S3.

Fig. 6 Distribution of mass ratio (top) and initial tilt angle (bottom) as functions of time before collision. The left and right panels show the models with varied density and varied shape, respectively. Collisions that occur before 5 h have elapsed are excluded from this plot. For the remaining simulations shown in the figure, the average time it takes to collide is 143 h, 166 h, and 133 h for models D1, D2, and D3, respectively, while the average is 171 h, 193 h and 79 h for models S1, S2, and S3 models, respectively.

Fig. 7 Escape times, averaged over the 45 initial values of θ0, as functions of the mass ratio. The left columns show the varied density models, while the right columns show the varied shape models. The error bars show the standard deviation of the escape times at the given mass ratio.

3.5 Escape cases

The escape cases are mainly found at the low end of the mass ratio spectrum, typically q < 0.2 for most models, as these low mass ratio systems have positive energies. Simulations that result in the secondary escaping make up 1.38%, 2.40%, and 4.86% of the simulations for models D1, D2, and D3, respectively. Thus, it appears that the secondary escapes more easily when the secondary is denser than the primary. Meanwhile, for the varied shape models we find that the escape cases make up 5.23%, 1.20% and 0.48% of the simulations for models S1, S2, and S3, respectively. The lower number of escape cases in models S2 and S3 is likely a consequence of a lower energy configuration in the system, due to the elongated shape of the secondary. However, the torque applied on the secondary, due to the primary, is stronger when it becomes more elongated. It is therefore possible that the small number of escaped secondaries is due to the early collisions.

How long it takes for the secondary to escape varies with both its density and its shape. In Fig. 7, we plot the escape time t_e, averaged over the 45 initial angles, as a function of q. From this figure, we can see that there is a trend that the secondary takes longer to escape as the mass ratio increases, which is similar to the findings of Boldrin et al. (2016). We find that the average escape time is roughly twice as short in model D2 as in the results of Boldrin et al. (2016) at corresponding mass ratios. However, as described in Sect. 2.2, the value of ω₀ becomes larger when the secondary is moved closer, due to an increase in θ₀, and the probability of an early escape increases as the system energy increases. The escape time trends of models D1 and D3 are similar to that of D2, but the escape times are slightly longer when density of the secondary is lower. The average escape times of models S1 and S2 are similar up to q = 0.11. However, the escape time increases significantly with mass ratio in model S3.

For systems where the secondary takes longer to escape, we expect that rotational energy gets transferred to translational energy before the secondary is expelled. At the time of escape (when the eccentricity exceeds 1), the separation between the two bodies is large enough for the rotational and translational motion to be decoupled (Scheeres 2002). Hence, we expect the rotation of the bodies to slow down as time passes in our simulations; after escape we expect the rotation period to stay roughly constant. Because it takes longer for the secondary to escape in systems with higher mass ratios, we expect the rotation of the primary to slow down more in systems with higher mass ratios. We first investigate the rotation of the primary after mutual escape. We calculate the (instantaneous) rotation period of a body as T = 2π/ω, where ω is the magnitude of the angular velocity of the body. The rotation period of the primary, at the time of escape, is displayed in the two top left panels in Fig. 8, showing the rotation period of the primary T_p at the final time step as a function of q. In the figure it can be seen that T_p in all six models, is longer at higher mass ratios after escape of the secondary, indicating a correlation between T_p and q. The Spearman correlation coefficients between T_p and q are shown in Table 1. For all models the correlation coefficients are r_s > 0.9. Furthermore, with the exception of model S3, the p-values are lower than 10⁻⁹. The high p-value in model S3 is likely due to the smaller number of escape scenarios for this model.

We also include the data of asteroid pairs from Pravec et al. (2019) in the figure, for pairs with q < 0.3 (indicated with gray crosses), and most of our results are within the range of the observed data. However, some outliers also exist in the data provided by Pravec et al. (2019), where some asteroid pairs have mass ratios that are too high and some pairs where the primary is rotating too slowly. Pravec et al. believe that these outlier asteroid pairs are not formed by rotational fission.

We also briefly studied the rotation period of the secondary after escape, which is shown in the two bottom left panels in Fig. 8. Unlike the primary, there are no obvious patterns of an increasing rotation period of the secondary when the mass ratio increases. We have also included the rotation period of the secondary of asteroid pairs from Pravec et al. (2019). With the exception of a few outliers in our results, most of the escaped secondaries have rotation periods that are in the range of the data from Pravec et al. (2019).

3.6 Unstable binaries

Some of our simulations, with positive total energy, are still in orbit around the primary after 200 days (the orange regions in Fig. 5). These systems are typically found near the same values of θ₀ as the escape cases, and we refer to them as “unstable”. Of all non-collision systems with positive energy, the unstable scenarios typically make up roughly one-half of them, with the exception of model D3 where the unstable cases make up approximately one-third of the simulations. However, we expect the number of unstable scenarios to decrease and to become either an escape or a collision case if a longer time span is considered.

In the third and fourth columns of Fig. 8, we give the rotation periods of the bodies at the end of the simulations for all unstable cases. The rotation periods of both bodies of these simulations, similar to the scenarios where the secondary has escaped, are also within the range of the observed data from Pravec et al. (2019). The primary is again seen to have longer rotation periods as the mass ratio increases.

Fig. 8 Distribution of rotation periods, as functions of mass ratio, for non-collision systems with positive total energy. The top and bottom rows correspond to the rotation period of the primary and secondary, at t = 200 days, respectively. The first and second columns show the escape cases, while the third and fourth columns show the unstable cases. The D and S models correspond to models D1–D3 and S1–S3, respectively. The gray crosses are data from Pravec et al. (2019), and only asteroid pairs with q < 0.3 are included in the figure.

3.7 Stable binaries

Finally, at mass ratios of q ≳ 0.2 the systems have negative total energy, forming binary systems that are stable against mutual escape. They correspond to the green regions in Fig. 5, and most of them appear at θ₀ ≲ 15°. Although this is called a stable orbit, the secondary may still collide with the primary if a longer time span is considered. Some systems with negative total energies do end up with an impact after 1000 h. The case with the longest time before impact (as seen in Fig. 6) is actually a system with negative total energy. However, we also saw in Fig. 6 that the time before impact is generally shorter at higher mass ratios. It is therefore possible that, for high enough mass ratios, systems that survive longer than ~ 100 h will never collide.

3.8 Rotational motion

In order to examine, the rotational state of the bodies at the end of the simulation, we follow Boldrin et al. (2016) and utilize the dynamic inertia I_D, defined as:

$$I_{\text{D}}=\frac{L^2}{2E_{\text{r}}}$$(14)

(Scheeres et al. 2000), where L is the magnitude of the angular momentum and E_r is the rotational kinetic energy of the body. A body has a uniform rotational motion when I_D = I_Z or I_D = I_x, which corresponds to rotations about the short and long axes, respectively¹. Nonuniform rotation (or tumbling motion) happens when I_x < I_D < I_z∙ This can be categorized as long-axis mode (LAM) when I_x < I_D < I_y, and as short-axis mode (SAM) when I_y < I_D < I_z (Scheeres et al. 2000). Here we only take into consideration the rotational motion in simulations that do not result in the two bodies impacting.

Initially, the primary has uniform rotational motion, where the dynamic inertia is equal to I_z, while the secondary starts off in a tumbling state. For low values of θ₀ the initial dynamic inertia of the secondary is close to I_z, and approaches I_x as θ₀ increases.

At the end of the simulations, we find that both the primary and the secondary, in most cases, are in some state of tumbling. Table 2 summarizes the final rotation state of the two bodies. The primary is mainly found with SAM rotation, which is close to its initial state. For q ≲ 0.25, the primary may be able to retain its uniform rotational motion throughout the whole simulation, and these are mainly found at mass ratios of q ≲ 0.2, as shown in Fig. 9. Most of these situations are found among the escape cases; however, some are also found among the unstable cases. This is a consequence of the secondary being unable to act with a gravitational torque on the primary due to the large separation between the bodies. This is similar to the results of Davis & Scheeres (2020), as they found that the spin state of the primary is, for the most part, unaffected when the secondary escapes. Moreover, simulations where the primary ends with a LAM rotation are more common at high mass ratios.

The secondary is also mostly in a tumbling state. Unlike the primary, LAM rotation is more common for the secondary because most simulations have a secondary with initial LAM rotation. Typically, the initial rotation mode of the secondary is SAM when θ₀ ≲ 27° and LAM otherwise, but it also depends on its shape. For the non-collision cases when θ₀ > 60°, nearly all simulations end with the secondary in a LAM rotation, as shown in Fig. 10, with one exception found in model D3. In some of the simulations (≲5%), the secondary has uniform rotational motion at the end of the simulation, either along the short or the long axis. These are mainly found when q ≤ 0.1, when the secondary has escaped, and when θ₀ = 0.001°. Uniform rotational motion is less common among the stable cases because both the primary and secondary act with torques on each other for a longer time period.

If we isolate the escaped secondaries in this analysis, we find that approximately 35–50% have SAM rotation at the end of the simulation for every model except the S3, where the percentage is 63% instead. Boldrin et al. (2016) found in their study that most escaped secondaries are SAM rotators. Our results are therefore slightly different in that we seem to find fewer with SAM rotation. In particular, we find fewer SAM rotators as the secondary becomes less dense. Davis & Scheeres (2020) also investigated the rotational state of escaped secondaries, and found that every escaped secondary is in tumbling motion.

We also wanted to study how the rotation period of the bodies changes with time when the secondary is still in orbit around the primary. Figure 11 shows the average rotation period of the primary and secondary as functions of time in the top and bottom rows, respectively. The left and right panels correspond to stable and unstable cases, respectively. The averaged data are binned in 48 hour periods.

In the figure, it can be seen that the average rotation period of the primary increases over time, both for the stable and unstable cases. Furthermore, the rotation period of the unstable cases are lower than the stable cases, which is a consequence of the large separations between the bodies, effectively decoupling the translational and rotational motions, similar to the escape cases. The secondary, as shown in the bottom two panels, has rotation periods of typically 10–15 h in the stable systems and, similarly to the primary, rotates slightly faster, typically 8–12 h in the unstable systems. The time evolution of the rotation period of the secondary is far more volatile within the first ~2500h of the simulations, and the figure shows that it experiences frequent speed-ups and slow-downs during this time period. After this, the rotation period of the secondary stabilizes.

We also see that the rotation period of the primary increases with mass ratio for the escape cases. Because the escape times are longer at higher mass ratios, the secondary can act with a gravitational torque for a longer time period.

We also show how the rotation periods change over time for four simulations with different outcomes for model D2. This is illustrated in Fig. 12. As previously mentioned, when the separation between the two bodies becomes large enough, the translational and rotational motion will decouple. As seen in the figure, for the escape and the unstable cases, when the bodies are sufficiently far apart, their rotation periods become approximately constant. For the stable and collision cases, the rotation periods vary far more as the bodies are relatively close to each other.

Fig. 9 Distribution of the rotation state of the primary, at the end of the simulation, as functions of q and θ0. The white regions correspond to simulations that result in collisions. The dashed lines indicate the value of q where the total energy is zero.

Fig. 10 Same as Fig. 9, but for the secondary.

Fig. 11 Average rotation period of the primary (top) and the secondary (bottom), averaged over all the stable (left) and unstable (right) cases, as functions of time. The averages, for t > 0, are binned over 48 hour intervals.

Fig. 12 How the separation between bodies and rotation periods changes over time, for four selected simulations. The top row shows the separation, while the middle and bottom rows show the rotation period of the primary and secondary, respectively. The mass ratios of the collision, stable, escape, and unstable examples are q = 0.29,0.32,0.15, and 0.19 respectively, with corresponding initial rotation periods of T0 = 4.87,5.87,5.42, and 5.56 h.

3.9 Secondary fission

Jacobson & Scheeres (2011) introduced secondary fission as a mechanism to form stable binaries from systems with low mass ratios. During secondary fission, the secondary disrupts or fissions when it is spun up by gravitational torques. Through secondary fission, parts of the energy in the system can be removed if the newly fissioned component escapes or impacts with the primary.

We wanted to investigate whether fission of the secondary can take place in our simulations, and similarly to Boldrin et al. (2016), we applied the rotation limit for surface disruption of the secondary as the critical limit for achieving secondary fission. We define this critical limit T_τ as the rotation rate at which a point mass is lifted off the surface by centrifugal forces. We use Eq. (4) with β = 1.0 to determine this limit. The value of T_r depends on the density, shape, and rotation state of the body. The rotation period required for secondary fission becomes longer when the density becomes smaller or when the body becomes more elongated. Tumbling motion may further increase the spin rate required for fission, and is taken into account during our analysis.

As is evident from the previous section, the average rotation period of the secondary has frequent speed-ups and slow-downs. The secondaries of some systems might obtain rotation periods short enough for secondary fission to occur. Figure 13 shows the percentages of simulations that experience secondary fission as functions of the mass ratio, based on the rotation criterion described above. Secondary fission events are most common when q = 0.01, and decrease as the mass ratio increases. These events may take place up to q = 0.4, with the exception of models S3, where the secondary can still disrupt at mass ratios as high as q = 0.72. The work of Jacobson & Scheeres (2011) and Sharma (2014) also suggests disruption events are more common if the body is more elongated. However, unlike the findings of Jacobson & Scheeres (2011), we find that secondary fission may occur also in systems with positive total energy.

As previously seen, the rotation of the secondary slows down rapidly during the first few hours of the simulation, and spins up again further into the simulation. It is therefore likely that secondary fission events occur early on in the simulation, but they may also take place toward the end of the 200-day simulations. We find that roughly half the secondary fission events may occur before 5 h have elapsed for most models, and for D1 and S3 the percentage is even higher, at 100% and 82%, respectively. Many of these events belong to simulations where the two bodies impact. Thus, for these systems, a ternary (or multiple) can be created early on, and may change the dynamics of the system, possibly preventing the early collisions.

Cases where secondary fission may occur after the initial 5 h are spread out in time. For some models the secondary can disrupt at t > 3000 h as the secondary’s rotation slowly speeds up over time, although the number of these events is small (less than ten in total).

Fig. 13 Percentage of cases that experience secondary fission as a function of mass ratio.

4 Discussion

Table 3 shows a summary of the percentage of each end-case scenario for the models presented. The collision cases make up approximately 80% of the simulations, while the remaining cases are categorized as stable, unstable, or escape. The collisions typically occur when the secondary has a tilt angle in the range θ₀ = 15°–80°. However, for tilt angles smaller or larger than this the system can develop into a stable binary, an unstable binary, or a system with an escaped secondary.

The difference in the end-case distribution does not change significantly when the density of the secondary is changed, but rather when the secondary take different shapes. By allowing the secondary to become more elongated, the number of collisions increases. In the model where the secondary’s shape is close to spheroidal (model S1), ~77% of the simulations end with an impact. This percentage increases to above 90% for the model with the most elongated secondary (model S3).

Most of the collision events take place very early in the simulations. We find that 90% of the collisions occur before 5 h have elapsed. This occurs because we move the secondary closer to the primary when it is rotated with an angle θ₀, such that the surface-to-surface distance is always 1 cm, as described in Sect. 2.2. One consequence of this is that the secondary rotates into the primary early in the simulation, due to the gravitational torque. The gravitational torque is also stronger on the secondary when it is more elongated, and hence the increased fraction of collision events in model S3 compared to S1 and S2. The early impact between the two bodies may help contribute to stabilizing the system. The energy dissipation from these collision events may prevent the secondary from escaping, and thus allow formation of asteroid binaries with low mass ratios. The early collisions we find are similar to the 1996 HW1 simulations and also shorter than the Moshup simulations of Davis & Scheeres (2020), who found that the median collision time is 2.1 h and 0.52 days, respectively.

One of our models is the same as the model used by Boldrin et al. (2016), and when comparing with their work, a larger percentage of our simulations end up with the two bodies impacting. This is another consequence of keeping the surface separation to 1 cm. Furthermore, because the surface-to-surface distance is always 1 cm, we find that collisions can occur at angles as low as θ₀ ~ 8°, while Boldrin et al. (2016) find that collisions do not occur when θ₀ ≲ 40°.

Escape scenarios, which is the likely mechanism behind the formation of some asteroid pairs (Pravec et al. 2010), exist for systems with low mass ratios, and we find that the time it takes for mutual escape to happen is longer when the mass ratio is higher. However, there are cases where the escape time is longer than 1000 h at low mass ratios, but these cases are not frequent. At the lowest mass ratios the escape time tends to be the longest when the secondary has a more elongated shape, as it was seen in model S3, because the energy configuration in S3 is lower than the other models at equal mass ratios. We also found that escape cases were more frequent when the secondary has a higher density, and asteroid pairs with secondaries of higher density may therefore be more frequent in the asteroid pair population.

Because we consider relatively short simulation times, some of the systems will remain as unstable systems throughout the duration of the simulation (200 days). These systems are generally found at intermediate mass ratios, but this will also vary based on the density of the secondary as well as its shape. If a longer time span is considered, such as 200 years as done by Boldrin et al. (2016), these unstable cases either become escape cases, or end up with an impact between the two bodies.

We find that the rotation period of the primary increases with time, hence it loses rotational energy because the rotational energy is converted to translational energy (Scheeres 2002). The rate at which the rotation period increases is slower for the unstable cases compared to the stable cases as the average separation between the bodies is larger for the former case. This is also seen among the escape cases. Higher mass ratios result in both longer escape times and longer rotation periods of the primary. Moreover, changing the shape of the secondary has a larger effect on the rotation of the primary in the stable cases, compared to changing its density. The average rotation period of the primary in model S3 can be nearly twice as long as in model S1.

The average rotation period of the primary is much longer in our simulations compared to that in some of the observed asteroid binaries. For instance, the rotation period of Moshup is estimated to be 2.76 h (Ostro et al. 2006) and 2.26 h for Didymos (Naidu et al. 2020), where the mass ratio of the former system is estimated to be q = 0.057 (Ostro et al. 2006) and q = 0.048 for the latter (Pravec et al. 2006). Observations by Pravec et al. (2016) estimate that the primary bodies have rotation periods of less than ~4.4 h. Meanwhile, the average rotation periods of the primary we find, for the stable cases, are in the range 15–25 h. Although, our simulation time span is very short, adding other physical effects such as tidal torques and the YORP effect may be able to allow the primary to spin up after a longer time period. On the other hand, rotation periods of the secondaries observed by Pravec et al. (2016) range from ~14 h all the way up to ~37 h, which is within the range of what we find in our results for the stable cases. However, the mass ratio of the binary systems presented by Pravec et al. (2016) are lower than 0.125 (assuming equal bulk densities), while our stable cases are found when q > 0.2. Energy dissipation of the system is therefore required, for example from collision or secondary fission.

We compare rotation periods from our simulations with that of observed asteroids pairs by Pravec et al. (2019), and find that there is an overall agreement for systems with q < 0.3, as illustrated in Fig. 8. The primaries of low mass ratio asteroid pairs were observed to be rapidly rotating, which indicates that the secondary may have escaped very early after the initial fission process. However, some systems observed by Pravec et al. (2019) have mass ratios that are too high or have a primary with a rotation period that is too long. These systems are considered outliers, and the rotational fission theory is unable to explain their existence (Pravec et al. 2019). Kyrylenko et al. (2021) suggest that the mass ratio of the asteroid pair 1999 XF200 and 2008 EL40, which reside in the main belt, is q < 0.01. The rotation period of 1999 XF200 is estimated to be 4.903 h², which is within the range of the escape rotation periods of the primary for q = 0.01 in our models. Furthermore, Kyrylenko et al. (2021) estimate that the age of this asteroid pair is 265.8 kyr. Under this time period, the rotation period of the bodies have likely changed by a significant amount due to the YORP effect and possibly also via collisions with other bodies in the main belt.

Jacobson & Scheeres (2011) find that the separation between positive and negative energy regimes can be approximated to q ≈ 0.2, and that it should not change much if the bodies are more elongated. We find that this separation regime can go as high as q = 0.29 when the secondary has twice the density of the primary (model D3), and as low as q = 0.17 when the secondary is more elongated (model S3). This indicates that asteroid pairs formed through rotational fission may occur at higher mass ratios, up to q ~ 0.3, if the secondary has a higher density than the primary, or if it becomes less elongated.

If the secondary also fissions, then ternary or multiple systems may be formed. If any of the components escape or collide with the primary, this can stabilize the system (Jacobson & Scheeres 2011). We find that this process generally occurs at low mass ratios, as predicted by Jacobson & Scheeres (2011), and also fits the findings of Boldrin et al. (2016). On the other hand, unlike the work of Jacobson & Scheeres (2011), we find that secondary fission may still occur in systems where the total energy is negative. We also find that it is more likely for the secondary to disrupt when it has a lower density or when it is more elongated. The latter possibility is in agreement with the work of Sharma (2014), who shows that more elongated bodies are less stable to finite structural perturbations compared to the less elongated ones. Observations of Pravec et al. (2016) find that there is a scarce number of binaries with secondary elongations of a_s /b_s ≳ 1.5. This suggests that elongated secondaries may experience multiple fission events, and may thus reshape over time. The results of Davis & Scheeres (2020) also suggest a form of energy dissipation, such as secondary fission, is required to stabilize the 1994 KW4 and 2000 DP107 systems in their current state.

Boldrin et al. (2016) used second-order spherical harmonics to study the dynamical evolution of fissioned systems, while we use an exact expression. Higher order terms become more important when the bodies are more elongated (Hou et al. 2017) or when the bodies are close. A future study that compares an exact method with an approximation may give better insights into the importance of exact mathematical expressions used to study asteroid systems immediately after fission.

Acknowledgements

We want to thank the anonymous referee for their valuable feedback that improved the manuscript.

References

Appendix A Ellipsoid potential

For any general ellipsoid with semiaxes a > b > c and constant density ρ, the gravitational potential is given by (MacMillan 1930)

$$\begin{array}{lll}\text{Φ}\left(r\right)\hfill & =\hfill & \frac{2\pi \rho abc}{\sqrt{a^2-c^2}}(\left[1-\frac{x^2}{a^2-b^2}+\frac{y^2}{a^2-b^2}\right]F\left({\omega }_{\kappa },k\right)\hfill \\ \hfill & +\hfill & \left[\frac{x^2}{a^2-b^2}-\frac{\left(a^2-c^2\right)y^2}{\left(a^2-b^2\right)\left(b^2-c^2\right)}+\frac{z^2}{b^2-c^2}\right]E\left({\omega }_{\kappa },k\right)\hfill \\ \hfill & +\hfill & \left[\frac{c^2+\kappa }{b^2-c^2}{y}^2-\frac{b^2+\kappa }{b^2-c^2}{z}^2\right]\frac{\sqrt{a^2-c^2}}{\sqrt{\left(a^2+\kappa \right)\left(b^2+\kappa \right)\left(c^2+\kappa \right)}}),\hfill \end{array}$$(A.1)

where F(ω_κ, k) and E(ω_κ, k) are respectively the elliptic integrals of the first and second kind, κ is the largest root of the equation

$$\frac{x^2}{a^2+\kappa }+\frac{y^2}{b^2+\kappa }+\frac{z^2}{c^2+\kappa }=1,$$(A.2)

and

$${\omega }_{\kappa }={\sin }^{-1}\sqrt{\frac{a^2-c^2}{a^2+\kappa }}$$(A.3)

$$k=\sqrt{\frac{a^2-b^2}{a^2-c^2}}.$$(A.4)

The components of the gravitational field g = ∇Φ then become

$$g_x=\frac{4x\pi \rho abc}{\sqrt{a^2-c^2}}\frac{E\left({\omega }_{\kappa },k\right)-F\left({\omega }_{\kappa },k\right)}{a^2-b^2},$$(A.5)

$$\begin{array}{lll}{g}_y\hfill & =\hfill & \frac{4y\pi \rho abc}{\sqrt{a^2-c^2}}[\frac{F\left({\omega }_{\kappa },k\right)}{a^2-b^2}-\frac{\left(a^2-c^2\right)E\left({\omega }_{\kappa },k\right)}{\left(a^2-b^2\right)\left(b^2-c^2\right)}\hfill \\ \hfill & +\hfill & \frac{\left(c^2+\kappa \right)}{b^2-c^2}\frac{\sqrt{a^2-c^2}}{\sqrt{\left(a^2+\kappa \right)\left(b^2+\kappa \right)\left(c^2+\kappa \right)}}],\hfill \\ \hfill & \hfill & \hfill \end{array}$$(A.6)

$$g_z=\frac{4z\pi \rho abc}{\sqrt{a^2-c^2}}\left[\frac{E\left({\omega }_{\kappa },k\right)}{b^2-c^2}-\frac{\left(b^2+\kappa \right)}{b^2-c^2}\frac{\sqrt{a^2-c^2}}{\sqrt{\left(a^2+\kappa \right)\left(b^2+\kappa \right)\left(c^2+\kappa \right)}}\right].$$(A.7)

Despite being functions of x, y, and z, the variable k is treated as a constant when the partial derivatives are taken (see MacMillan 1930, for details).

Appendix B Verification of accuracy

The accuracy of the integration scheme can be demonstrated by inspecting the conservation of total energy E, total linear momentum p, and total angular momentum J. This is shown in Fig. B.1 for one of the models (D2 with q = 0.32, θ₀ = 6.1°). In the figure we plot, for each of these three quantities, the difference between the initial value at t = 0 and the value at each subsequent time step. For the energy and angular momentum, the difference is normalized by the initial values E₀ and J₀. We find that these quantities are conserved to the 11th decimal place; the error on the linear momentum fluctuates between the 4th and 7th decimal places.

Fig. B.1 Illustration of energy, angular, and linear momentum conservation for one of the simulations.

All Tables

All Figures

Fig. 1 Contact binary before separation (top) and after separation (bottom). Top: Configuration of the contact binary the moment of fission where the gravitational force FG and the centrifugal force FC are equal. The cross indicates the center of mass of the system. The long and short semiaxes (a and c) are aligned with the body-fixed x- and z-axes of the respective bodies. Bottom: Angular velocities of the bodies after the contact binary separates.

In the text

	Fig. 2 How the secondary moves closer to the primary when the angle θ0 of the secondary is increased.
In the text

Fig. 3 How angular velocity changes with the initial tilt angle θ0 of the secondary. Top panel: ω0 is affected when there is contact between the surfaces (blue) and when the surface-to-surface distance increases (red) as the secondary is tilted (see also Fig. 2). The mass ratio is q = 0.1 for the data in the top panel. Bottom panel: relative difference in ω0 between the two approaches shown in Fig. 2, but for three different mass ratios q. The mass ratio is defined as q = ms/mv.

In the text

	Fig. 4 Total energy of the system as a function of the mass ratio q, for each model, using θ0 = 0.001°. The gray dotted line shows the zero energy line.
In the text

	Fig. 5 Model outcomes as functions of mass ratio q and initial tilt angle θ0. In the top row, the left, middle, and right panels correspond to models D1, D2, and D3, respectively. In the bottom row, the left, middle, and right panels correspond to models S1, S2, and S3. The dashed lines indicate the value of q where the total energy is zero.
In the text

Fig. 6 Distribution of mass ratio (top) and initial tilt angle (bottom) as functions of time before collision. The left and right panels show the models with varied density and varied shape, respectively. Collisions that occur before 5 h have elapsed are excluded from this plot. For the remaining simulations shown in the figure, the average time it takes to collide is 143 h, 166 h, and 133 h for models D1, D2, and D3, respectively, while the average is 171 h, 193 h and 79 h for models S1, S2, and S3 models, respectively.

In the text

	Fig. 7 Escape times, averaged over the 45 initial values of θ0, as functions of the mass ratio. The left columns show the varied density models, while the right columns show the varied shape models. The error bars show the standard deviation of the escape times at the given mass ratio.
In the text

Fig. 8 Distribution of rotation periods, as functions of mass ratio, for non-collision systems with positive total energy. The top and bottom rows correspond to the rotation period of the primary and secondary, at t = 200 days, respectively. The first and second columns show the escape cases, while the third and fourth columns show the unstable cases. The D and S models correspond to models D1–D3 and S1–S3, respectively. The gray crosses are data from Pravec et al. (2019), and only asteroid pairs with q < 0.3 are included in the figure.

In the text

	Fig. 9 Distribution of the rotation state of the primary, at the end of the simulation, as functions of q and θ0. The white regions correspond to simulations that result in collisions. The dashed lines indicate the value of q where the total energy is zero.
In the text

	Fig. 10 Same as Fig. 9, but for the secondary.
In the text

	Fig. 11 Average rotation period of the primary (top) and the secondary (bottom), averaged over all the stable (left) and unstable (right) cases, as functions of time. The averages, for t > 0, are binned over 48 hour intervals.
In the text

Fig. 12 How the separation between bodies and rotation periods changes over time, for four selected simulations. The top row shows the separation, while the middle and bottom rows show the rotation period of the primary and secondary, respectively. The mass ratios of the collision, stable, escape, and unstable examples are q = 0.29,0.32,0.15, and 0.19 respectively, with corresponding initial rotation periods of T0 = 4.87,5.87,5.42, and 5.56 h.

In the text

	Fig. 13 Percentage of cases that experience secondary fission as a function of mass ratio.
In the text

	Fig. B.1 Illustration of energy, angular, and linear momentum conservation for one of the simulations.
In the text