# Using a geometrical algorithm to provide N-body initial conditions for the gravitational phase of asteroid family formation

Ronald-Louis Ballouz, Kevin Walsh, Derek Richardson, Patrick Michel

► To cite this version:

Ronald-Louis Ballouz, Kevin Walsh, Derek Richardson, Patrick Michel. Using a geometrical algorithm to provide N-body initial conditions for the gravitational phase of asteroid family formation. Monthly Notices of the Royal Astronomical Society, 2019, 485 (1), pp.697-707. 10.1093/mnras/stz410. hal-02407853

HAL Id: hal-02407853

https://hal.science/hal-02407853

Submitted on 12 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Using a geometrical algorithm to provide *N*-body initial conditions for the gravitational phase of asteroid family formation

Ronald-Louis Ballouz<sup>1,2,3*</sup> Kevin J. Walsh,<sup>4</sup> Derek C. Richardson<sup>3</sup> and
Patrick Michel<sup>5</sup>

Accepted 2019 February 7. Received 2019 February 4; in original form 2018 October 3

## ABSTRACT

Numerical studies of asteroid family formation use the size-frequency and ejection-speed distributions of the reaccumulated fragments as the main constraints to validate asteroid collision models. However, when shape and spin are also considered, they add new constraints that must be matched simultaneously. While this poses new challenges, it also serves to increase the reliability of numerical models. Coupled with advances in shape and spin determination in ground and space-based observations, numerical simulations are poised to become a crucial tool for understanding asteroid-family-forming events and to infer internal properties of family members. Numerical simulations have typically relied on a two-stage process for modelling these events. First, a hydrocode models the progenitor's fragmentation. In this study, the considered hydrocode is based on smoothed-particle hydrodynamics (SPH) techniques. Secondly, the SPH output is fed to an *N*-body code to model the reaccumulation or escape of fragments. Here, we explore the upgraded capabilities for the second stage. We use a soft-sphere discrete element method to accurately model the contact forces as fragments reaccumulate. However, SPH simulations typically result in large particle overlaps, some of which have large speeds that drive the time-step requirement to prohibitively low values. We introduce a novel approach for handling the transition between SPH and *N*-body by using a computational geometry technique called *α*-shape modelling. This technique will enable future studies to explore the large parameter space of asteroid family formation in order to link observed asteroid family properties with specific collision scenarios and to probe asteroid material and internal properties.

**Key words:** methods: numerical – minor planets, asteroids: general.

## 1 INTRODUCTION

When large asteroids experience a catastrophic impact, a huge number of fragments can result that, with the help of self-gravity, can reaccumulate into a family of smaller asteroids. Previous research used numerical simulations to show that the reaccumulation of fragments following catastrophic impacts can match many properties of orbit- and spectra-linked asteroid families (Michel et al. 2001; Michel, Benz & Richardson 2003; Durda et al. 2007). These events have also been found to produce orbiting debris that can explain some of the satellites found around large asteroids (Durda et al. 2004).

One of the most important insights from these studies is the expectation that all but the largest asteroids are likely reaccumulated fragments from a previous catastrophic collision (Michel et al. 2001).

These numerical studies matched the resulting size-frequency distribution (SFD) and ejection-speed distribution of the reaccumulated fragments with observations of asteroid families in order to constrain models of asteroid collisions (Michel et al. 2001, 2003). Numerical simulations of asteroid collisions typically proceed in two stages. First, a hydrodynamics code is used to simulate the shock and fracture of the target asteroid. Secondly, a gravitational *N*-body code is used to simulate the gravitational reaccumulation of the fragments. These studies that initially linked the formation of asteroid families with catastrophic impacts used a simple contact model that

<footnote n="1"><sup>1</sup>Institute of Space and Astronautical Science, Japanese Aerospace Exploration Agency, Sagamihara, Kanagawa Prefecture, 252-5210, Japan</footnote>

<footnote n="2"><sup>2</sup>Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85705, USA</footnote>

<footnote n="3"><sup>3</sup>Department of Astronomy, University of Maryland, College Park, MD 20742, USA</footnote>

<footnote n="4"><sup>4</sup>Southwest Research Institute, 1050 Walnut St. Ste 300, Boulder, CO 80302</footnote>

<footnote n="5"><sup>5</sup>Lagrange Laboratory, Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, CS 34229, F-06304 Nice Cedex 4, France</footnote>

<footnote>* E-mail: rballouz@orex.lpl.arizona.edu</footnote>
merged colliding particles into spherical particles, conserving mass and momentum (Michel et al. 2001). This was done for simplicity and computational efficiency. However, fragment shapes and the effect of variations in the material properties of the individual particles could not be investigated with this approach. Retaining the shape and spin information after reaccumulation provides additional constraints that must be matched simultaneously. Increasing the number of constraints for numerical experiments makes these models more diagnostic of the collisional event and reaccumulation process.

More recently, Michel & Richardson (2013) used the rigid aggregate formulation of Richardson et al. (2009), which allows a cluster of particles to form non-spherical shapes, to reproduce the shape of the asteroid (25143) Itokawa (Fujiwara et al. 2006). This was a significant advancement from simple merging, but it had the issue of assuming material parameters and a primitive cohesion model. While this newer model was able to provide a qualitative idea of the possible shapes that can form after reaccumulation, it was not clear whether the contact model would be able to reproduce shapes and spin using realistic material properties. Hence, it would be difficult to tie observations of asteroid shape, spin, and geotechnical properties to a specific collision scenario.

In this study, we will demonstrate upgraded capabilities for the *N*-body calculation of a reaccumulation study. Rather than merging particles or allowing them to form rigid aggregates, we employ a soft-sphere discrete element method (SSDEM) that allows us to explicitly model the contact forces between realistic material types, enabling us to establish between a wider range of asteroid observables with their formation scenario. Here, we will explain and demonstrate the upgraded capabilities for the second stage, the *N*-body calculations.

## 1.1 Important physics

For our work, the collision phase of the impact between two asteroids is modelled with a 3D smoothed-particle hydrodynamics (SPH) code that models shock propagation in elastic solids. It is typically found that in the impact energy regimes experienced by asteroids, the target is fully damaged by the end of these simulations. However, new work (El Mir, Ramesh & Richardson 2018) using a material point method has shown that there are some scenarios in which the target is not fully damaged.

Following the impact phase of the simulation, essentially once fragmentation ceases, the simulations are handed off to the code *pkggrav* (Richardson et al. 2000; Stadel 2001; Leinhardt et al. 2002) for the *N*-body gravitational phase modelling. *pkggrav* is a parallelized hierarchical tree code, which has been modified to detect and resolve low-speed collisions between particles. This code has been used extensively for modelling gravitational aggregates, or ‘rubble-pile’ asteroids during disruption and reaccumulation (see Walsh, Michel & Richardson 2011 for a review). These codes, SPH and *N*-body, have been used extensively in tandem to simulate asteroid collision and disruption (Michel et al. 2001, 2002; Michel, Benz & Richardson 2003, 2004; Durda et al. 2004, 2007; Benavidez et al. 2011).

## 1.2 Issues

While these previous studies have successfully been able to reproduce the SFD and ejection-speed distributions of certain asteroid families, the simplicity of the collision models used has limited their ability to match other observables, such as the shape and spin of

small bodies. These works have relied on the perfect merging of particles and hard-sphere interactions to resolve collisions, or the rigid aggregate approach noted earlier. These types of interactions do not properly capture the real contact physics, making comparisons with realistic complex asteroid shapes uncertain. The implementation of SSDEM into *pkggrav* (Schwartz, Richardson & Michel 2012) led to more accurate modelling of asteroid geophysics in static (Zhang et al. 2017, 2018) and low-speed kinetic (Ballouz et al. 2014, 2015) regimes. In SSDEM, particles are allowed to slightly overlap in order for their collisions to be treated by a spring and dashpot model (Cundall & Strack 1979) that can account for frictional contact forces. Recent work by Schwartz et al. (2018) demonstrated that an *N*-body code with SSDEM for collision-handling can reproduce the shapes of bilobate comets. However, before an SSDEM implementation can be properly used for an asteroid reaccumulation simulation, a few issues in the handoff from the SPH simulation need to be resolved:

(i) Large overlaps: At the end of a typical hydrodynamics simulation of an asteroid collision there are significant overlaps between particles. In the SPH formulation, state variables such as pressure are smoothed over neighbouring particles so that the discrete approach mimics more continuum-like behaviour of the fluid. The particles therefore typically overlap, sometimes strongly. For a scheme that hands off the end-state of an SPH simulation to an *N*-body code, this is highly problematic. For a HSDEM code overlaps are unphysical and by default stop the integration. In SSDEM, large overlaps result in unrealistically strong restoring forces. Typical methods to resolve this issue for the handoff are particle mergers and/or shrinking particle radii. Decreasing particle radii for the handoff leads to unphysically high particle density (over 10 grams per cubic cm in some cases), and eliminating particles to make space leads to loss of linear and angular momentum. Somehow the state of the system–mass density and velocity field–needs to be maintained and represented with individual and roughly similarly sized particles with only minimal (1 per cent) overlaps.

(ii) Resolution (number of particles): The importance of resolution in a numerical simulation is dependent on the simulation method and application. For SPH simulations of hypervelocity impacts, higher resolutions translate to a better solution of the damage propagation and energy dissipation, improving realism. However, for the gravity-phase/reacummulation component of a simulation, an increase in the resolution of the simulation does not have an obvious or direct correlation to an improvement in realism. While resolving smaller particles to mimic a size distribution found in real asteroids can influence the equilibrium shapes determined in simulations (Walsh, Richardson & Michel 2012), the important contact physics in DEM simulations are typically parameterized through macroscopic dissipation, friction, and cohesion models that attempt to capture the grain-scale physics (Schwartz et al. 2012; Sánchez & Scheeres 2016; Zhang et al. 2017, 2018). For DEM simulations, this parameterization of grain-scale physics is necessary as simulating global asteroid geophysics (1–100 km scales) in relevant time-scales (days to weeks of simulated time) with realistic grain sizes (mm to cm) would be computationally infeasible with modern computational hardware and software. Therefore, transitioning from SPH to *N*-body requires a more thoughtful consideration of resolution that optimizes the divergent needs of each tool. A number of suites of hydrodynamic simulations have been done at very high resolution (Genda et al. 2015; Jutzi, Michel & Richardson 2019) that would be too large for a project that intends to survey a wide range of parameter space.
Furthermore, since these techniques and codes are relatively new, we aim to use techniques where we can test convergence at differing resolutions. See Asphaug et al. (2015) for further discussion of this topic.

These issues motivated the development of a technique that could transform the state of a hydrodynamics simulation into the initial conditions of an N-body model, eliminating overlaps and scaling down to the user's desired resolution.

In Section 2, we outline the methodology for the hydrocode and reaccumulation simulations, and then describe the new method of transitioning between them. In Section 3, we describe the results of the reaccumulation simulations that have been prepared using this new method. In Section 4, we discuss our results and provide some insights into the future of reaccumulation simulations as a way to investigate asteroid family formation and the structures of asteroid family parent bodies.

## 2 METHODOLOGY

### 2.1 SPH code

To compute the fragmentation phase of the collision, we use a smoothed particle hydrodynamics (SPH) code. The results of these simulations are described in Jutzi et al. (2019). The SPH simulations include a model adapted for porous materials (Jutzi, Benz & Michel 2008; Jutzi et al. 2009). This model is based on the P-alpha model initially proposed by Herrmann (1969). In the P-alpha model, the volume change in a porous material is treated in two parts: the pore collapse, and the compression of the material composing the matrix. This is done by introducing a distention parameter, which is the ratio of the density of the porous material and the corresponding solid material. The distention parameter is used in the computation of the pressure and the deviatoric stress tensor. By including this property, damage then increases as a result of both crack activation and change in the distention. Material parameters based on porous pumice were used as these were successfully validated by comparison with laboratory impact experiments (Jutzi et al. 2009). A detailed description of the SPH model and its implementation can be found in Jutzi et al. (2008). The SPH simulation resolution was fixed at 400 000 particles. The complete suite of results is presented in Jutzi et al. (2019).

### 2.2 N-body collision code: pkdgrav

Originally collisions in pkdgrav were treated as idealized single-point-of-contact impacts between rigid spheres. A SSDEM option was added by Schwartz et al. (2012). In SSDEM, particle contacts can last many time-steps, with reaction forces dependent on the degree of overlap and contact history. By allowing particles to overlap, multicontact and frictional forces can be modelled. The inclusion of these contact forces preserves shape and spin information of reaccumulated gravitational aggregates, which leads to better constraints on the geotechnical properties of the reaccumulated fragment.

This implementation of SSDEM uses a spring/dash-pot mechanism for the contact forces. In this model, two overlapping particles will experience a normal contact force, $F_N$, and a tangential stick-slip force, $F_T$. In the current implementation, contact forces are determined by restoring forces with optional damping and/or friction. The following formulation describes the normal and tangential

restoring forces with a static friction component included:

$$F_N = -k_N x \hat{\mathbf{n}} + C_N \mathbf{u}_N, \qquad (1)$$

$$F_T = \min(k_T \delta_T + C_T \mathbf{u}_T, \mu_S |F_N| \delta_T / |\delta_T|),$$

where $k_N$ and $k_T$ are the normal and tangential spring constants, respectively. The damping parameters ($C_N$, $C_T$) are related to the conventional normal and tangential coefficients of restitution used in hard-sphere implementations, $\epsilon_N$ and $\epsilon_T$ (see Schwartz et al. 2012 for further details). $x$ is the mutual overlap of the two particles, $\delta_T$ is the sliding displacement from the equilibrium contact point, and $\hat{\mathbf{n}}$ is the unit vector from the particle's centre to its neighbour. $\mathbf{u}_N$ and $\mathbf{u}_T$ are the normal and tangential components of the relative velocity between the two particles, respectively. The coefficient of static friction, $\mu_s$, determines the maximum amount of tangential force that can be supported by the contact point.

Previous studies (Iwashita & Oda 1998; Mohamed & Guitterez 2010) have shown the inclusion of rotational resistances into a contact model allows numerical simulations to better match laboratory experiments. In our implementation, two rotational resistances are modelled: twisting and rolling resistance. These are parameterized with the coefficients $\mu_T$ and $\mu_T$, respectively. Twisting resistance arises from the slip and friction at the contact region due to a difference in the rotation rate of the particles in a direction along the normal vector $\hat{\mathbf{n}}$. Rolling resistance can come about from a number of sources: slip and friction on the contact surface, viscous hysteresis, shape effects, and/or surface adhesion. Both resistances are modelled as functions of the relative twisting or rotational angular speed that can reach critical values based on a shape parameter, $\beta$, after which they maintain a maximum resisting torque (see section 2.2 of Zhang et al. 2017). $\beta$ represents a statistical measure of real particle shape, which has been found to be one of the primary physical mechanisms for rotational resistance (Jiang, Shen & Wang 2015) that also plays an important role in modelling interparticle cohesion (Zhang et al. 2018).

The SSDEM implementation of pkdgrav has been validated through comparison with laboratory experiments (Schwartz et al. 2014) and has been used for various Solar System applications: size-sorting on asteroids (Maurel et al. 2014), avalanche dynamics (Yu et al. 2014), the stability of asteroid rotations (Zhang et al. 2017, 2018), and collisions between rubble piles at low speeds (Ballouz et al. 2014, 2015).

### 2.3 Transitioning from SPH to N-body

Due to the nature of SPH simulations (see Section 1.2), and the need to transfer to a N-body gravity code in less than 1 dynamical time-scale (in those cases where no gravity is modelled in the SPH phase; see Section 2.4), there can be significant particle overlaps. Previous studies (Michel et al. 2001) used 'perfect merging' for a single time-step to handle these particle overlaps. However, this technique typically leads to large variations in particle sizes and densities. Often the resulting densities are unrealistically small for asteroidal material (<0.5 g cm⁻³). Here, we present a new method to resolve the issues currently present in the typical ad hoc handoff procedure to generate initial conditions for the N-body particle code.

#### 2.3.1 α-shape modelling of the disrupted core

The hydrocode phase of a hypervelocity impact models the mechanical response of the target to the impact-generated shockwave that propagates throughout the body. This typically lasts some tens
<figure><img src="image_3.png" /><figcaption>**Figure 1.** α-shape models of the asteroid Eros for different values of α, generated from a shape model of the asteroid derived from spacecraft data (Gaskell 2008). The vertices defined in the shape model are taken as input, and the α-shape algorithm attempts to find shapes that are bound by this set of points. The effect of changing the value of α is shown, with α decreasing from top to bottom and left to right. The top left shape is the convex hull (α = +∞). The bottom right shape is for a value of alpha smaller than the resolution of the original shape model. In this case, the shape is disassembled, and singular triangles, which do not bound the interior, are shown.</figcaption></figure>

of seconds in simulated time (corresponding to the time required for the shockwave to travel through the target and reflect back) and results in a disrupted body made up of two components: (1) a disrupted core with a cavity and (2) ejecta propagating away from the impact point. Our strategy is to substitute the core (which contains a majority of the particles and particle overlaps) with non-overlapping particles that fill the same shape and have the equivalent internal velocity distribution as the original core. In order to do this, we must first discriminate between SPH particles that make up the core and those that make up the ejecta.

For sub-catastrophic collisions, the disrupted core is easily identifiable since the spherical shape of the pre-impacted target remains mostly intact, except for a crater at the point of impact. For some catastrophic and super-catastrophic collisions, it is sometimes difficult to distinguish the boundary between the core-region and the ejecta; therefore, we adopt a core-selection criterion based on coordination number of particles (the number of contacts each particle has). The core is defined to be the region contained within a sphere which has a radius set equal to the distance from the centre of mass where the average particle begins to have a coordination number of 3 or less. This criterion essentially allows us to distinguish between the high-number-density core with its large overlaps and number of contacts, and the dilute ejecta region. We find that this selection process typically identifies a core region that has >95 per cent of all overlapping particles, regardless of the impact energy.

Once the core particles have been identified, we construct a geometric shape model of the core using a technique called ‘α-shape modelling’. α-shape modelling (Edelsbrunner & Mücke 1994) is a computational geometry algorithm that enables the construction of concave shapes that encompass a set of points, S, in Euclidean space. Similar to a convex hull, the shape that is generated is a triangulated tiling. The ‘α’ in the name denotes the single parameter the algorithm requires to generate the shape, and how tightly the boundary fits around the set of points. Intuitively, any facet of the triangulation only exists if a sphere with radius α can occupy its space and enclose a subset of the S. For α = +∞, the α-shape is identical to the convex hull of S. However, as α decreases, the shape shrinks, developing cavities. We demonstrate these properties of α-shapes for the shape model of Eros in Fig. 1. For modelling the shape of a gravitational aggregate, α is related to the ‘resolution’ of

the simulation, i.e. the typical particle size in the SPH simulation. In this work, we set α to 2.5 × the mean particle radius.

Once a shape model of the asteroid’s core is generated, it can be used to create an aggregate of particles, free of overlaps, with the same shape and with a user-specified resolution. In this work, we generate hexagonal-close-packed (HCP) and random close packing (RCP) aggregates of the same shape as that of disrupted-core’s alpha shape. The RCP aggregates are created by simulating the gravitational collapse of a cloud of randomly distributed monodisperse particles. In this manner, the heavily overlapped particles in the SPH-generated core are replaced with non-overlapping particles. In some instances, we find a small incidence of particles (fewer than 1 in 40 000 particles across all cases considered) that overlap in the ‘ejecta’ region. These rare remaining particle overlaps are resolved through perfect merging, conserving linear and angular momentum.

### 2.3.2 Recreating the velocity distribution

Finally, the velocity distribution of the original SPH output is reproduced in the initial conditions generated by this procedure. We do this by comparing the positions of the original SPH particles to the newly generated particles. We use a tree-code to perform nearest neighbour searches of particles in the new α-shape generated core (list-A) for each particle in the original list of SPH particles (list-B). We restrict the neighbor search to be within 1.1 × particle radius of each particle’s centre. Since these conditions force a monodisperse particle size distribution, this condition effectively finds all particle overlaps. However, for a polydisperse distribution, the average particle radius can be used, and a second conditional check can follow to ensure that the particles are indeed overlapping. The procedure continues with each particle in list-B transferring its linear and angular momentum to its nearest neighbour in list-A that is within 1.1 × particle radius of its centre. If a particle in list-B finds no neighbours within 1 particle radius, typically true for particles in the ‘ejecta’, it is added to list-A with its original velocity and spin. The tree is reconstructed whenever a new particle is added to list-A.

Fig. 2 shows a visual comparison of the speed of particles in the original initial conditions and the α-shape processed initial conditions for the case of a 5 km s⁻¹ impact on to a 100 km target at an angle of 60° by a 18 km-size impactor. Fig. 3 shows the distribution of particle speeds for the original SPH output compared with the α-shape processed initial conditions for two different resolutions. We demonstrate that the new technique can adequately reproduce the particle speeds of the original output. The differences between the original and processed data are most noticeable for low particle speeds, reflecting the resolution reduced core. Higher particle speeds are those of ejecta that are found at large distances from the core and are thus never pre-processed as part of this technique, hence a perfect match.

To summarize, the proposed new technique has the following protocol:

(i) Isolate ‘core’ from the rest of the ‘ejecta’.

(ii) Generate a shape model for core.

(iii) Generate an aggregate of particles through procedural-generation or gravitational collapse.

(iv) Use shape model to carve out the aggregate into the new resampled core.

(v) Use a nearest neighbour search to map the linear and angular momentum from original SPH file to the resampled core, while incrementally adding in particles from the ejecta and resolving the small number of overlaps in the ejecta through merging.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The α-shape processing technique can reproduce the speed distribution of particles. Here, we show the reduced resolution output (left) and the original (right) for the case of a 5 km s⁻¹ impact on to a 100 km target at an angle of 60° by a 18 km-size impactor. See the online version of this article for a color version of this figure.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Distribution of particle speeds for the case of a 5 km s⁻¹ impact on to a 100 km-diameter target at an angle of 60° by a 18 km-diameter impactor. The dashed blue curve shows the speed distribution of the particles at the end of the SPH simulation. The dotted green and the dot-dashed red curve represent the speed distributions of particles following the α-shape processing technique for two different replacement cores, resulting in final resolutions of 100 000 and 50 000 particles, respectively. We demonstrate that the technique can adequately reproduce the speed distribution of particles.</figcaption></figure>

We evaluated the change in the total energy after handoff and find that it varies by an order of 1 per cent for all resolution cases. This change reflects a variation in the gravitational potential energy as large particle overlaps in the SPH output lead to high values of the potential energy. We also evaluated the variability in the total energy across different resolutions and found this to be of the order of 0.1 per cent, which reflects the slight variations in the particle spatial distribution in the disrupted core that our method produces.

We caution special-handling for more energetic super-catastrophic impacts. Careful processing of the ‘core identification’ is required in order to isolate it from the rest of the ejecta. In these cases, the core is shifted from the centre of mass of the system, which makes it difficult to automate its identification. For such situations, we suggest identifying the core by visually inspecting the spatial distribution of particles. In the following section, we describe the results of N-body simulations that follow this newly developed handoff from SPH.

## 2.4 Simulation parameters

For the simulations presented here, the target had a radius of $R_{targ}$ ~ 50 km. The impactor had an average radius of $R_{imp}$ that ranged between 5 and 18 km (see Jutzi et al. 2019 for a more comprehensive account of the SPH simulation output that we use). The particle densities, $\rho$, were set to 2.7 g cm⁻³ and the bulk density of the target’s core varied depending on the type of packing that was used for the replacement core. We used HCP cores as a default for most of the simulations, but also looked at the effect of having a core made up of RCP particles of different packing fractions ($\phi$) that ranged from 0.53 to 0.59.

The collisional properties of the constituent particles are specified prior to each simulation. The coefficients of restitution were fixed at $\epsilon_n = 0.5$ and $\epsilon_t = 0.5$ (moderately dissipative collisions). Furthermore, since SSDEM models treat particle collisions as reactions of springs due to particle overlaps, the magnitude of the normal and tangential restoring forces are determined by the spring constants $k_N$ and $k_T = \frac{2}{7}k_N$. For rubble-pile collisions, the value of $k_N$ is constrained by

$$k_N = m \left( \frac{v_{max}}{x_{max}} \right)^2, \qquad (2)$$

where m corresponds to the typical mass of the most energetic particles, $v_{max}$ is the maximum expected speed in the simulation,
<table><caption>**Table 1.** SPH and *pkggrav* simulation parameters.</caption><thead><tr><th>SPH parameter</th><th>Value</th></tr></thead><tbody><tr><td>Equation of state</td><td>Tillotson, Basalt (Melosh 1989)</td></tr><tr><td>ρ<sub>material</sub></td><td>2.327 g cm<sup>-3</sup></td></tr><tr><td>P-alpha parameter</td><td>Pumice (Jutzi et al. 2009)</td></tr><tr><td>N</td><td>400 000</td></tr><tr><td>pkggrav parameter</td><td>Value or range</td></tr><tr><td>ρ</td><td>2.7 g cm<sup>-3</sup></td></tr><tr><td>ε<sub>n</sub></td><td>0.5</td></tr><tr><td>ε<sub>t</sub></td><td>0.5</td></tr><tr><td>k<sub>n</sub></td><td>8 × 10<sup>13</sup> kg s<sup>-2</sup></td></tr><tr><td>k<sub>t</sub></td><td>2.3 × 10<sup>13</sup> kg s<sup>-2</sup></td></tr><tr><td>μ<sub>s</sub></td><td>0–1.03</td></tr><tr><td>μ<sub>t</sub></td><td>0–1.31</td></tr><tr><td>μ<sub>r</sub></td><td>0–1.0</td></tr><tr><td>β</td><td>0.5</td></tr><tr><td>N</td><td>44 231–99 291</td></tr></tbody></table>

and $x_{max}$ is the maximum fractional particle overlap, which we set to be 1 per cent of the particle radius. The particle radius in the simulations are monodisperse, but vary by a factor of 1.25 from the highest resolution to the lower resolution simulation. Thus, we set the value of $k_n$ based on the particle radius of the highest resolution case, and keep it the same for consistency across simulations of different resolution. For our cases, $k_n \sim 8 \times 10^{13} \text{ kg s}^{-2}$. In order for the post-collision system to reach a steady state, the total run-time was set to at least 100 × the dynamical time for the system, $\text{sim } 1/\sqrt{G\rho} = 0.45 \text{ h}$, for a total simulation time of at least 45 h. Furthermore, a time-step $\Delta t = 0.03 \text{ s}$ was chosen on the basis of the time required to sample particle overlaps adequately, for the choice of $k_n$ and $x_{max}$ given above.

In order to evaluate the relative effects of interparticle friction through a DEM, we performed simulations for three sets of friction parameters. The nominal friction parameters are based on the values obtained by Jiang et al. (2015) for rough sand, and used by Zhang et al. (2017) for simulations of rubble-pile spin stability, and correspond to an internal angle of friction that ranges from 30° to 43° depending on the internal packing of the rubble pile. We take the coefficients of friction and shape factor used in those studies, but decrease the static friction by half. The set of parameter values for this material type are: $\mu_s = 0.5$, $\mu_r = 1.3$, $\mu_t = 1.0$, and $\beta = 0.5$. We also performed simulations with no interparticle friction (we term this the frictionless parameter set), and for intermediate values of friction, where all coefficients are set to zero except the static friction coefficient, which is set to 0.3 (we term this the $\mu_s = 0.3$ parameter set). We also performed single-case simulations for cases where only $\mu_s$ was varied from the [$\mu_s = 0.3$] case, and $\beta$ was changed from the nominal case. The simulation parameter values are summarized in Table 1.

## 3 RESULTS

Table 2 summarizes the results of a series of N-body simulations for a wide range of collision scenarios. These simulations used as input the final state of SPH impact simulations of two asteroids where the target had a diameter $d_{targ} = 100 \text{ km}$ and the impact speed was $v_{imp} = 5 \text{ km s}^{-1}$. The impact angle, $\theta_{imp}$, ranged between 15° and 60°, and the impactor diameter, $d_{imp}$, ranged between 5 and 18 km. The majority of the simulations presented here were for the specific impact scenario of an 18 km-diameter impactor striking the target

at an angle of 60°. For this case, we studied the effect of varying the simulation resolution (N) and the packing structure of the processed core in order to evaluate the sensitivity of the collisional outcomes to the α-shape processing of SPH simulations. Furthermore, we studied the effect of varying the interparticle friction properties on collision outcomes. The collision outcomes are characterized by the largest remnant's mass, $M_{lr}$, normalized by the total system mass, $M_{tot}$, and the spin period of the largest remnant, $P_{LR}$, reported in Table 2.

In Table 3, we describe the final shape of the largest remnant for each impact simulation. We report 5 metrics that describe the overall final shape and provide a text description of the shape based on visual inspection. The metrics we report are:

(i) intermediate-to-major axes ratio: $b/a$,

(ii) minor-to-major axes ratio: $c/a$,

(iii) the maximum projected sphericity, a measure of the roundness of the shape (Ehlmann, Viles & Bourke 2008), $\Psi = (\frac{c^2}{ab})^{1/3}$, where $\Psi \in [0, 1]$, and the shape tends towards a spherical shape as $\Psi \to 1$ and a flat shape as $\Psi \to 0$,

(iv) the form factor, a measure of the flatness of the shape (Ehlmann, Viles & Bourke 2008), $F = \frac{a-b}{a-c}$, where $F \in [0, 1]$, and the shape tends towards prolate as $F \to 1$ and oblate as $F \to 0$,

(v) the 3D solidity, a new metric we define that measures the largest remnant's relative concavity. Traditionally, 'solidity' is used for expressing grain morphology (e.g. Avery, Panter & Gorsevski 2017) and it represents the ratio of the area of the 2D projection of the grain to the area of its 2D convex hull. A high solidity value of 1 represents a smooth grain edge while a low solidity value of 0 represents a rough grain edge. Here, 3D solidity is defined as ratio of the 3D α-shape of the largest remnant to its 3D convex hull. A high 3D solidity value (tending towards unity) represents a smooth featureless object, and a low value represents larger and larger cavities.

### 3.1 Comparison of collision outcomes: effect of resolution

In Fig. 4, we show the collisional outcomes (largest remnant mass and final spin period) and shape properties of the largest remnants for cases 1, 6, and 8–18. These represent collisions with an impact angle of 60° and an impactor radius of 18 km for a range of N, $\phi$, and friction properties.

#### 3.1.1 Largest remnant mass and spin properties

In the top two panels of Fig. 4, we show that the mass of the largest remnant stays roughly constant across the resolution that we studied here. The largest deviations in the final remnant are for the cases with nominal friction, with the highest resolution cases having a slightly less massive largest remnant than lower resolution cases. For the two sets of HCP cases (N = 50 909 and 99 291), the increase in N leads to a decrease in the mass of the largest remnant by about 3 per cent $M_{tot}$. This variation is even less sensitive for the frictionless case, where the difference is less than 1 per cent $M_{tot}$. The slight decrease in largest remnant mass with resolution for higher friction cases is likely driven by the enhanced stochasticity introduced with interparticle-friction forces and larger N. However, this effect seems to diminish for $N > 60 000$.

We find a larger variation in the final spin period with resolution; however, the changes are still within reasonable limits. Moreover, the results show that friction more sensitively affects the final shape
<table><caption>**Table 2.** Simulation collisional outcomes: normalized largest remnant mass and spin period. The packing fraction, $\phi$, of the HCP cases were 0.70 ± 0.01.</caption><thead><tr><th>Case</th><th>θ<sub>imp</sub> (°)</th><th>d<sub>imp</sub> (km)</th><th>N</th><th>Friction properties</th><th>Packing</th><th>M<sub>lr</sub>/M<sub>tot</sub></th><th>P<sub>LR</sub> (h)</th></tr></thead><tbody><tr><td>1</td><td>60</td><td>18</td><td>50 909</td><td>Frictionless</td><td>HCP</td><td>0.792</td><td>9.02</td></tr><tr><td>2</td><td>60</td><td>18</td><td>50 909</td><td>μ<sub>s</sub> = 0.3</td><td>HCP</td><td>0.787</td><td>10.06</td></tr><tr><td>3</td><td>60</td><td>18</td><td>50 909</td><td>μ<sub>s</sub> = 0.6</td><td>HCP</td><td>0.782</td><td>10.38</td></tr><tr><td>4</td><td>60</td><td>18</td><td>50 909</td><td>μ<sub>s</sub> = 0.9</td><td>HCP</td><td>0.767</td><td>9.91</td></tr><tr><td>5</td><td>60</td><td>18</td><td>50 909</td><td>β = 0.2</td><td>HCP</td><td>0.777</td><td>10.70</td></tr><tr><td>6</td><td>60</td><td>18</td><td>50 909</td><td>Nominal</td><td>HCP</td><td>0.785</td><td>10.61</td></tr><tr><td>7</td><td>60</td><td>18</td><td>50 909</td><td>β = 0.8</td><td>HCP</td><td>0.778</td><td>10.93</td></tr><tr><td>8</td><td>60</td><td>18</td><td>99 291</td><td>Nominal</td><td>HCP</td><td>0.752</td><td>10.76</td></tr><tr><td>9</td><td>60</td><td>18</td><td>99 291</td><td>Frictionless</td><td>HCP</td><td>0.783</td><td>8.64</td></tr><tr><td>10</td><td>60</td><td>18</td><td>41 073</td><td>Nominal</td><td>RCP, φ = 0.59</td><td>0.782</td><td>11.47</td></tr><tr><td>11</td><td>60</td><td>18</td><td>41 073</td><td>Frictionless</td><td>RCP, φ = 0.59</td><td>0.786</td><td>8.97</td></tr><tr><td>12</td><td>60</td><td>18</td><td>41 073</td><td>μ<sub>s</sub> = 0.3</td><td>RCP, φ = 0.59</td><td>0.784</td><td>9.77</td></tr><tr><td>13</td><td>60</td><td>18</td><td>88 307</td><td>Nominal</td><td>RCP, φ = 0.64</td><td>0.756</td><td>12.57</td></tr><tr><td>14</td><td>60</td><td>18</td><td>88 307</td><td>Frictionless</td><td>RCP, φ = 0.64</td><td>0.778</td><td>9.73</td></tr><tr><td>15</td><td>60</td><td>18</td><td>88 307</td><td>μ<sub>s</sub> = 0.3</td><td>RCP, φ = 0.64</td><td>0.773</td><td>11.26</td></tr><tr><td>16</td><td>60</td><td>18</td><td>78 060</td><td>Nominal</td><td>RCP, φ = 0.53</td><td>0.757</td><td>11.19</td></tr><tr><td>17</td><td>60</td><td>18</td><td>78 060</td><td>Frictionless</td><td>RCP, φ = 0.53</td><td>0.778</td><td>8.81</td></tr><tr><td>18</td><td>60</td><td>18</td><td>78 060</td><td>μ<sub>s</sub> = 0.3</td><td>RCP, φ = 0.53</td><td>0.775</td><td>9.65</td></tr><tr><td>19</td><td>45</td><td>9</td><td>57 772</td><td>Nominal</td><td>HCP</td><td>0.759</td><td>15.28</td></tr><tr><td>20</td><td>45</td><td>13</td><td>76 639</td><td>Nominal</td><td>HCP</td><td>0.270</td><td>23.90</td></tr><tr><td>21</td><td>30</td><td>7</td><td>51 943</td><td>Nominal</td><td>HCP</td><td>0.682</td><td>27.74</td></tr><tr><td>22</td><td>30</td><td>7</td><td>51 943</td><td>μ<sub>s</sub> = 0.3</td><td>HCP</td><td>0.682</td><td>26.58</td></tr><tr><td>23</td><td>30</td><td>9</td><td>70 052</td><td>Nominal</td><td>HCP</td><td>0.275</td><td>35.18</td></tr><tr><td>24</td><td>15</td><td>5</td><td>44 231</td><td>Nominal</td><td>HCP</td><td>0.981</td><td>21.72</td></tr><tr><td>25</td><td>15</td><td>5</td><td>44 231</td><td>Frictionless</td><td>HCP</td><td>0.924</td><td>32.77</td></tr><tr><td>26</td><td>15</td><td>5</td><td>44 231</td><td>μ<sub>s</sub> = 0.3</td><td>HCP</td><td>0.925</td><td>33.67</td></tr><tr><td>27</td><td>15</td><td>7</td><td>48 603</td><td>Nominal</td><td>HCP</td><td>0.471</td><td>31.86</td></tr><tr><td>28</td><td>15</td><td>7</td><td>48 603</td><td>Frictionless</td><td>HCP</td><td>0.471</td><td>42.30</td></tr><tr><td>29</td><td>15</td><td>7</td><td>48 603</td><td>μ<sub>s</sub> = 0.3</td><td>HCP</td><td>0.467</td><td>35.82</td></tr></tbody></table>

(and hence spin) and not the remnant mass, which agrees with expectations. For the nominal friction cases, the variations in the period with resolution have a standard deviation of 0.78 h. This standard deviation is even lower for the frictionless case (0.4 h). Furthermore, this variation is likely more sensitive to the internal packing structure of the core, rather than the resolution. For the two sets of HCP cases, the difference in the final spin period for the nominal and frictionless cases are 1.3 per cent and 4.2 per cent, respectively. These results demonstrate the utility and efficiency of our α-shape method in decreasing the computational cost of reaccumulation simulations. For the same number of processors, the time taken for a `pkggrav` run to complete scales as $N \log N$. Thus, a reduction in the number of particles from 400 000 to 40 000 represents a speed-up of more than an order of magnitude.

Furthermore, our results reveal a relationship between interparticle friction and spin period. For this specific impact scenario, friction acts to diminish the final spin period. We find that across all resolutions, there is a strict correlation between interparticle friction and final spin period, with the frictionless case having the lowest $P_{LR}$, and the highest friction case having highest $P_{LR}$. Therefore, while the final spin period can likely be used to discriminate against the internal properties of an asteroid, these results suggest that there is an upper limit to the amount of angular momentum that can be imparted in a collision, which is set by a frictionless case. Note that the frictionless cases reported here are for particles with no static, sliding, and rolling friction; however, these particles still have some tangential resistance due to the tangential damping coefficient. The value of the internal angle of friction can be further lowered by

considering a polydisperse distribution of particles that may act as lubricant, increasing the likely limit of the largest remnant’s spin.

### 3.1.2 Shape properties

In the bottom four panels of Fig. 4, we show the influence of resolution and interparticle friction on the shape properties of the largest remnants. We find that for b/a, c/a, and Ψ, there is little to no variation across all resolution cases, demonstrating the capability of our method to reproduce consistent shapes for the largest remnant, even for much lower resolutions. The values for F are also consistent for different values of N for the nominal friction cases. However, there is a noticeable variation in F for the frictionless case for the larger values of N > 70 000, with the largest difference being 30 per cent. For these cases, the variation may be driven by the difference in the packing structure of the remnant rather than N, with higher packing structures leading to shapes that are relatively more prolate. For similar packing structure (HCP) in the frictionless cases, the variation in F is only 4.3 per cent in comparison. Overall, we find that our technique can robustly reproduce the shapes of largest remnant even for lower values of N.

## 3.2 Effect of friction and impact energy on shape

All 4 metrics suggest that friction strongly influences and determines the resulting shape of the largest remnant (across all resolutions studied here). For the nominal value of the friction
<table><caption>**Table 3.** Simulation collisional outcomes: largest remnant shapes.</caption><thead><tr><th>Case</th><th>b/a</th><th>c/a</th><th>Ψ</th><th>F</th><th>3D solidity</th><th>Description</th></tr></thead><tbody><tr><td>1</td><td>0.9104</td><td>0.8363</td><td>0.9159</td><td>0.5475</td><td>0.9655</td><td>Egg</td></tr><tr><td>2</td><td>0.7962</td><td>0.6633</td><td>0.8206</td><td>0.6053</td><td>0.9522</td><td>Veggie dog</td></tr><tr><td>3</td><td>0.7676</td><td>0.6054</td><td>0.7816</td><td>0.5890</td><td>0.9532</td><td>Veggie dog</td></tr><tr><td>4</td><td>0.8097</td><td>0.6253</td><td>0.7846</td><td>0.5078</td><td>0.9496</td><td>Veggie dog</td></tr><tr><td>5</td><td>0.8355</td><td>0.5184</td><td>0.6852</td><td>0.3415</td><td>0.9352</td><td>Pancake</td></tr><tr><td>6</td><td>0.8450</td><td>0.4844</td><td>0.6524</td><td>0.3006</td><td>0.9277</td><td>Pancake</td></tr><tr><td>7</td><td>0.8608</td><td>0.4843</td><td>0.6483</td><td>0.2699</td><td>0.9149</td><td>Pancake</td></tr><tr><td>8</td><td>0.8258</td><td>0.4928</td><td>0.6650</td><td>0.3436</td><td>0.9345</td><td>Pancake</td></tr><tr><td>9</td><td>0.9165</td><td>0.8541</td><td>0.9267</td><td>0.5719</td><td>0.9608</td><td>Egg</td></tr><tr><td>10</td><td>0.8187</td><td>0.4363</td><td>0.6149</td><td>0.3216</td><td>0.9181</td><td>Pancake</td></tr><tr><td>11</td><td>0.9031</td><td>0.8141</td><td>0.9020</td><td>0.5213</td><td>0.9671</td><td>Egg</td></tr><tr><td>12</td><td>0.7813</td><td>0.6750</td><td>0.8354</td><td>0.6729</td><td>0.9609</td><td>Veggie dog</td></tr><tr><td>13</td><td>0.8181</td><td>0.4554</td><td>0.6329</td><td>0.3339</td><td>0.9114</td><td>Pancake</td></tr><tr><td>14</td><td>0.9140</td><td>0.8186</td><td>0.9017</td><td>0.4741</td><td>0.9604</td><td>Egg</td></tr><tr><td>15</td><td>0.7579</td><td>0.7009</td><td>0.8654</td><td>0.8095</td><td>0.9493</td><td>Veggie dog</td></tr><tr><td>16</td><td>0.8268</td><td>0.4519</td><td>0.6274</td><td>0.3160</td><td>0.9103</td><td>Pancake</td></tr><tr><td>17</td><td>0.9246</td><td>0.8093</td><td>0.8914</td><td>0.3955</td><td>0.9634</td><td>Egg</td></tr><tr><td>18</td><td>0.7504</td><td>0.6898</td><td>0.8591</td><td>0.8048</td><td>0.9549</td><td>Veggie dog</td></tr><tr><td>19</td><td>0.7903</td><td>0.4244</td><td>0.6108</td><td>0.3642</td><td>0.9271</td><td>Pancake</td></tr><tr><td>20</td><td>0.7813</td><td>0.7355</td><td>0.8847</td><td>0.8270</td><td>0.8945</td><td>Rubber ducky</td></tr><tr><td>21</td><td>0.9286</td><td>0.5446</td><td>0.6836</td><td>0.1569</td><td>0.9405</td><td>Pancake</td></tr><tr><td>22</td><td>0.9050</td><td>0.6238</td><td>0.7548</td><td>0.2526</td><td>0.9479</td><td>Bean/pancake</td></tr><tr><td>23</td><td>0.7683</td><td>0.6409</td><td>0.8116</td><td>0.6451</td><td>0.8985</td><td>Lumpy bean</td></tr><tr><td>24</td><td>0.9013</td><td>0.5680</td><td>0.7100</td><td>0.2284</td><td>0.9236</td><td>Pancake</td></tr><tr><td>25</td><td>0.9007</td><td>0.8572</td><td>0.9344</td><td>0.6953</td><td>0.9745</td><td>Egg</td></tr><tr><td>26</td><td>0.9055</td><td>0.5915</td><td>0.7284</td><td>0.2313</td><td>0.933</td><td>Pancake</td></tr><tr><td>27</td><td>0.6577</td><td>0.6079</td><td>0.8252</td><td>0.8731</td><td>0.9161</td><td>Bean</td></tr><tr><td>28</td><td>0.9633</td><td>0.9256</td><td>0.9617</td><td>0.4934</td><td>0.9748</td><td>Spheroid</td></tr><tr><td>29</td><td>0.7999</td><td>0.7102</td><td>0.8575</td><td>0.6903</td><td>0.9665</td><td>Egg</td></tr></tbody></table>

parameter, the shapes tend to appear as slightly curved flat plates (or pancakes), with a relatively high *b/a* and low *c/a*. We show the cross-sections of these shapes in the upper half of Fig. 5. These cross sections also demonstrate that the overall shapes of the largest remnant are not sensitive to the change in resolution (with slight differences arising from changes in the orientation of the rubble pile due to the automated nature of our plotting scheme). The lower half of Fig. 5 illustrates that the independence of final shape on resolution is also true for different values of the friction coefficient.

We performed a series of simulations for a variety of impactor sizes and impact angles (Cases 19–29) in order to evaluate the range of shapes that can be produced for different specific impact energies. In Fig. 6, we show the shape properties as a function of the mass of the largest remnant. Overall, we find that for frictionless cases, the largest remnant shapes are spherical with values of *b/a*, *c/a*, and $\Psi$ close to 1.0, regardless of the impact energy.

When we included some non-zero friction, we found that the shape parameters become slightly dependent on the impact scenario. The dashed blue lines in Fig. 6 show the general trend for the nominal friction cases, and we report the goodness of fit through the coefficient of determination, $r^2$ value, for each shape parameter. We found a slightly positive, although weak, correlation between *b/a* and mass of the largest remnant, and a stronger negative correlation for the other three shape parameters. Of these, the value of *F* seems to be the most influenced by changes in the mass of the largest remnant, although the scatter in values does appear quite high. While the outcomes of the moderate friction runs ($\mu_s = 0.3$) were not included in the fitting, they do appear to follow the same general trends seen for the moderate friction cases. The trend in *F* would

suggest that objects become more oblate for less energetic impacts; however, looking at the cross-sections of these largest remnants gives a more complete story.

We show the cross-sections of some of our outcomes for nominal friction parameters in Fig. 7. As the mass of the largest remnant increases, its shape becomes rounder in the *a–b* plane ($F \to 0$), and flatter in the *a–c* plane ($\Psi \to 0$). This indicates flattening of the target along one direction. In other words, these largest remnants have ‘pancake’ shapes with some added complexity. Furthermore, the largest remnants of catastrophic collisions ($M_{lr}/M_{tot} < 0.5$) exhibit more complexity in their shapes, resulting in what are best described as bean-shaped remnants, similar to the asteroid (87) Sylvia (Marchis et al. 2005), and, in one case, the rubber-ducky type shape of comet 67P/Churyumov-Gerasimenko (Jorda et al. 2015). This is reflected in the values of the 3D solidity for these objects (see Table 3). We find that friction decreases the earlier used 3D solidity of objects as concavities are able to more easily form. Furthermore, the cases with high mass loss have the lowest values of 3D solidity, as the high-energy collisions allow the object to completely reshape.

## 4 DISCUSSION AND CONCLUSIONS

Through simulations of the reaccumulation process of asteroidal material following catastrophic disruption, studies have shown that hypervelocity impacts are an important process for the formation of asteroid families. In this study, we have demonstrated a new technique to execute the handoff of the output of the hypervelocity impact (modelled via SPH simulations) to the input of an *N*-body simulation that solves the longer term evolution of subsequent material that
<figure><img src="image_6.png" /></figure>

**Figure 4.** Mass, spin, and shape properties of the largest remnant following the reaccumulation simulation for the cases 1, 6, and 8–18 (see Table 2). Regardless of the friction parameters used, we see the resolution of a reaccumulation simulation has very little influence on the largest remnant’s mass (top left), spin period (top right), *b/a* (middle left), *Ψ* (middle right), *c/a* (bottom left), and *F* (bottom right). The red squares, green crosses, and blue circles represent simulations where the interparticle friction had zero, intermediate ($\mu_s$ = 0.3), and nominal values, respectively. See the text for discussion of these results.

either escapes (asteroid family members) or reaccumulates, helping to shape the largest remnant (asteroid family parent body) of a collision.

Previous studies have shown the ability of N-body simulations to predict the collisional scenario that formed particular asteroid families and match the SFD and ejection-speed distribution of asteroid families (Michel et al. 2001, 2003, 2004). With our new technique, we show that we can reconfigure the output of SPH simulations in order to: (1) use a SSDEM code to more accurately model particle collisions and (2) compute the longer reaccumulation process more efficiently.

The handoff is necessary for modelling both the gravitational dynamics and granular physics involved for the less-energetic reaccumulation process compared to the hypervelocity collision. By using SSDEM, we can begin to use the shapes and spins of asteroids as factors for discriminating the collision scenarios in asteroid family formation (Walsh et al. 2018) and also infer material and internal properties of asteroid families. For example, the Sylvia family has a parent body that is in a triple system, and was observed to have a distinct bean shape (Marchis et al. 2005). We can now

<figure><img src="image_7.png" /></figure>

**Figure 5.** Cross-sections of the largest remnants for cases with nominal friction values (top panels) and no friction (bottom panels). The simulation resolution increases from left to right. This shows that, for the range of resolution studied, the reaccumulation simulations lead to similar shape outcomes of the largest remnant regardless of the interparticle friction values. The cases in the top panels, from left to right, are: 11, 1, 17, 14, and 9. The cases in the bottom panels, from left to right, are: 10, 6, 16, 13, and 8.

<figure><img src="image_8.png" /></figure>

**Figure 6.** Shape properties of the largest remnant following the reaccumulation simulation for cases 19–29 (see Table 2) that vary in impactor size and impact angle, leading to variation in the final mass of the largest remnant. We find that the mass of the largest remnant, which may be regarded as a proxy for the specific impact energy, has some dependency on the shape parameters *b/a* (top left), *Ψ* (top right), *c/a* (bottom left), and *F* (bottom right). The red squares, green crosses, and blue circles represent simulations where the interparticle friction had zero, intermediate ($\mu_s$ = 0.3), and nominal values, respectively. The blue dashed lines represent the best linear fit to the nominal friction parameter outcomes, and we report the goodness of fit, $r^2$, on the top right of each individual panel. See the text for discussion of these results.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Cross-sections of the largest remnants for cases with different specific impact energies. The impact becomes less catastrophic from left to right. The cases shown here, in order from left to right, are: 20, 27, 21, 19, and 24.</figcaption></figure>

use the orbital, SFD, and shape of the parent body to constrain collision scenario and geotechnical properties of the asteroid. Future work will use the wealth of observational data on asteroid shapes through light curve inversion techniques (e.g. Durech, Sidorin & Kaasalainen 2010), radar shape modelling (e.g. Naidu et al. 2016), and spacecraft observations (e.g. Gaskell 2008; Gaskell et al. 2008), to investigate the origin of specific asteroid families.

We have demonstrated the sensitivity of the final shapes of the largest remnants to SSDEM parameters, which can be correlated to realistic material properties (e.g. Zhang et al. 2017). In general, frictionless cases (smooth values) tend to form oblate shapes, regardless of the impact energy. We showed that interparticle friction is necessary to reproduce the complexity in asteroid shapes. We find an indication that the shapes and masses of the largest remnants are correlated for simulations where friction is included. Sub-catastrophic collisions tend to produce flattened pancake shapes; in contrast, catastrophic and super-catastrophic collisions can produce more complex irregular shapes (beans and rubber ducky). The pancake shapes are produced in two scenarios: (1) a high-impact-angle collision that results in one hemisphere being sheared off the target and (2) a head-on collision that creates a large impact crater that effectively flattens the target in the direction of impact. In both cases, the relatively high interparticle friction (nominal values) allows mass loss to be localized to a limited region within the target. It is unclear whether these outcomes are representative of the true population, as there seem to be a dearth of actual asteroids with these shapes (though the amount is non-zero, see Walsh et al. 2018). The nominal friction parameters reflect those of sand grains used in laboratory experiments (Jiang et al. 2015) and were also used to simulate the rotational stability of asteroids (Zhang et al. 2017). These nominal friction simulations perhaps indicate that, while they can be used to study the dynamical evolution of a static system such as spin-state stability, they are likely too conservative for the interaction of large boulders feeling much higher shear stresses. With the eventual goal of tying observations of asteroid shapes and spins to collisional formation scenarios and asteroid material types, a concentrated study that can calibrate the friction properties of real materials with their simulated behaviour in high-energy regimes needs to first be completed.

Furthermore, parent body heterogeneity would influence the contact dynamics and gravitational evolution following a family forming collision. Currently, it is difficult to determine how differentiated a 100 km-diameter asteroid could or should be. Observational studies of asteroid families suggest that asteroid family members are homogeneous bodies, but this may be biased by the inherent selection criteria of asteroid family membership (Nesvorný 2015). Are

there dynamical members of a family with different spectral types that may indicate asteroid heterogeneity? If so, then simulations of asteroid family formation need to be performed with material heterogeneity. The friction properties of the interior, which can be parameterized or related to material viscosity, would be different than the exterior of the asteroid.

Finally, a key feature of our new technique is the reduction in *N* (resolution). Our study across different *N* showed that the overall shapes of the largest remnant can be reproduced well even at low resolutions (*N* = 40 000), which represents a reduction in the number of particles by an order of magnitude from the corresponding SPH simulation in this study. This demonstrates that one can reliably use this technique to study the reaccumulation process with a discrete-element method much more effectively. Furthermore, we argue that despite the reduction in resolution between the SPH and *N*-body segments, the outcomes should still robustly represent the asteroid family formation process. The resolution in DEM does not represent the same thing as it does in a hydrocode (Asphaug et al. 2015). For DEM, there is a ‘correct’ resolution to the problem, since every particle represents a discrete body. For these simulations, the correct resolution is set by the average boulder size of real rubble-pile asteroids. The influence of smaller grains (asteroid regolith) on the reaccumulation process is parameterized through the SSDEM contact parameters. In contrast, increasing resolution only serves to improve the realism of an SPH simulation, as damage propagation and energy dissipation can be better solved. As computing power inevitably improves, the resolution of SPH codes will reach levels that will be difficult for the slower DEM codes to match. Thus, we argue that special post-processing of hydrocode output, such as that demonstrated in this study, will be an essential tool for the adequate study of asteroid family formation and collisions in the future.

## ACKNOWLEDGEMENTS

RLB and KJW acknowledge support from National Aeronautics and Space Administration grant NNH12ZDA001N. RLB also acknowledges support from the Japanese Aerospace Exploration Agency's Aerospace Project Research Associate program.DCR acknowledges support from National Aeronautics and Space Administration grant NNX15AH90G. PM acknowledges support from the French space agency Centre National d'Études Spatiales. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. Simulations were run on the University of Maryland's Deepthought2 supercomputer, and the Bluecrab supercomputer cluster operated by the Maryland Advanced Research Computing Center under the auspices of the Johns Hopkins University and University of Maryland.

## REFERENCES

* Asphaug E., Collins G., Jutzi M., 2015, in Michel P., DeMeo F. E., Bottke W. F., eds, Asteroids IV. Univ. Arizona Press, Tucson, p. 661

* Avery M., Panter K. S., Gorsevski P. C., 2017, J. Volcanol. Geothermal Res., 332, 1

* Ballouz R.-L., Richardson D. C., Michel P., Schwartz S. R., 2014, Astrophys. J., 789, 158

* Ballouz R.-L., Richardson D. C., Michel P., Schwartz S. R., Yu Y., 2015, Planet. Space Sci., 107, 29

* Benavidez P. G., Durda D. D., Enke B. L., Bottke W. F., Nesvorný D., Richardson D. C., Asphaug E., Merline W. J., 2011, Icarus, 219, 57
Cundall P. A., Strack O. D. L., 1979, Geotechnique, 29, 47

Durda D. D., Bottke W. F., Enke B. L., Merline W. F., Asphaug E., Richardson D. C., Leinhardt Z. M., 2004, Icarus, 170, 243

Durda D. D., Bottke W. F., Nesvorný D., Enke B., Merline W. J., Asphaug E., Richardson D. C., 2007, Icarus, 186, 498

Durech J., Sidorin V., Kaasalainen M., 2010, A&A., 513, 13

Edelsbrunner H., Mücke E. P., 1994, ACM Trans. Graph., 13, 43

Ehlmann B. L., Viles H. A., Bourke M. C., 2008, J. Geophys. Res., 113, F02012

El Mir C., Ramesh K. T., Richardson D. C., 2019, Icarus, 321, 1013

Fujiwara A. et al., 2006, Science, 312, 1330

Gaskell R. W., 2008, Gaskell Eros Shape Model V1.0. NEAR-A-MSI-5-EROSSHAPE-V1.0. NASA Planetary Data System

Gaskell R. et al., 2008, Gaskell Itokawa Shape Model V1.0. HAY-A-AMICA-5-ITOKAWASHAPE-V1.0. NASA Planetary Data System

Genda H., Fujita T., Kobayashi H., Tanaka H., Abe Y., 2015, Icarus, 262, 58

Herrmann W., 1969, J. Appl. Phys., 40, 2490

Iwashita K., Oda M., 1998, J. Eng. Mech., 124, 285

Jiang M., Shen Z., Wang J., 2015, Comput. Geotechnol., 65, 147

Jorda L., Gaskell R. W., Capanna C., Hviid S., Gutierrez P. et al., 2015, SPC SHAP2 CARTESIAN PLATE MODEL DSK FOR COMET 67P/CG 6K PLATES, RO-C-MULTI-5-67P-SHAPE-V1.0 (all resolutions) NASA Planetary Data System and ESA Planetary Science Archive

Jutzi M., Benz W., Michel P., 2008, Icarus, 198, 242

Jutzi M., Michel P., Hiraoka K., Nakamura A. M., Benz W., 2009, Icarus, 201, 802

Jutzi M., Michel P., Richardson D. C., 2019, Icarus, 317, 215

Leinhardt Z. M., Richardson D. C., 2002, Icarus, 159, 306

Marchis F., Descamps P., Hestroffer D., Berthier J., 2005, Nature, 436, 822

Maurel C., Ballouz R.-L., Richardson D. C., Michel P., Schwartz S. R., 2017, MNRAS, 464, 2866

Melosh H. J., 1989, Impact Cratering: a Geologic Process. Oxford University Press, New York

Michel P., Richardson D. C., 2013, A&A, 554, L1

Michel P., Benz W., Tanga P., Richardson D. C., 2001, Science, 294, 1696

Michel P., Tanga P., Benz W., Richardson D. C., 2002, Icarus, 160, 10

Michel P., Benz W., Richardson D. C., 2003, Nature, 421, 608

Michel P., Benz W., Richardson D. C., 2004, Icarus, 168, 420

Mohamed A., Gutierrez. M., 2010, Granular Matter, 12, 527

Naidu S. P., Benner L. A. M., Margot J.-L., Busch M. W., Taylor P. A., 2016, AJ, 152, 99

Nesvorný D., Brož M., Carruba V., 2015, in Michel P., DeMeo F. E., Bottke W. F., eds, Asteroids IV 297-321. Univ. Arizona Press, Tucson

Richardson D. C., Quinn T., Stadel J., Lake G., 2000, Direct Large-Scale N-Body Simulations of Planetesimal Dynamics. Icarus, 143, 45

Richardson D. C., Michel P., Walsh K. J., Flynn K. W., 2009, Planet. Space Sci., 57, 183

Sánchez P., Scheeres D. J., 2016, Icarus, 271, 453

Schwartz S. R., Richardson D. C., Michel P., 2012, Granular Matter, 14, 363

Schwartz S. R., Michel P., Richardson D. C., Yano H., 2014, Planet. Space Sci., 103, 174

Schwartz S. R., Michel P., Jutzi M., Marchi S., Zhang Y., Richardson D. C., 2018, Nature Astron., 2, 379

Stadel J. G., 2001, Cosmological N-body simulations and their analysis in Thesis (PhD.), Univ. Washington, Washington, 141

Walsh K. J., Michel P., Richardson D. C., 2011, Adv. Sci. Lett., 4, 311

Walsh K. J., Richardson D. C., Michel P., 2012, Icarus, 220, 514

Walsh K. J. et al., 2019, Icarus, in review

Yu Y., Richardson D. C., Michel P., Schwartz S. R., Ballouz R.-L., 2014, Icarus, 242, 82

Zhang Y. et al., 2017, Icarus, 294, 98

Zhang Y., Richardson D. C., Barnouin O. S., Michel P., Schwartz S. R., Ballouz R.-L., 2018, AJ, 857, 15

This paper has been typeset from a *TEX/LATEX* file prepared by the author.