© ESO 2019

1 Asteroid collisions in the main belt

The main belt of asteroids is a collisional system. The breakups of asteroids have been recorded in the form of asteroid families (Nesvorný et al. 2015; Vinogradova 2019). We can also see impact features, such as craters or boulders, on surfaces of asteroids. These features can be observed directly during a satellite flyby or even with ground-based instruments (Vernazza et al. 2018; Fétick et al. 2019).

Physical processes during asteroid collisions are rather complicated for purely analytical estimates to yield precise quantitative results, as it is necessary to model a propagation of shock wave in the target, crack growth and consequent fragmentation, gravitational reaccumulation of ejected fragments, etc. Fragmentation of targets due to hyper-velocity impacts has been studied using laboratory experiments (Nakamura & Fujiwara 1991; Morris & Burchell 2017; Wickham-Eade et al. 2018). While the experiments can produce valuable constraints, the results cannot be directly compared with the breakups of asteroids, as the sizes of targets and kinetic energies of the impact would have to be extrapolated over several orders of magnitude. Experiments also do not take into account a gravitational reaccumulation. The collisions of asteroids are therefore commonly studied using numerical methods; the experiments then provide the calibration for the respective numerical codes.

Common methods for studying the collisions can be divided into particle-based (for example the N-body code pkdgrav, see Richardson et al. 2000) and shock-physics ones, such as mesh-based methods (used by code iSALE, see Suetsugu et al. 2018), or the smoothed particle hydrodynamics (Jutzi et al. 2015), used in this work. This is a Lagrangian, gridless method, which makes it suitable for impact simulations, as the computational domain is not a priori known. Especially in a hit-and-run impact, fragments of the projectile can travel to considerable distances from the target. In the smoothed particle hydrodynamics (SPH), the distant fragments do not require any special handling (they only might affect the performance of the code). SPH is also versatile, allowing the relatively easy implementation of new physics. The model of fragmentation can be straightforwardly incorporated into SPH, but it would be a difficult task for grid-based methods.

An outcome of a collision depends on a number of parameters, namely the diameter D_pb of the target, the diameter d_imp of the projectile, the specific impact energy Q, the impact angle ϕ_imp, but also the rotational periods of the colliding bodies, their shapes, material properties, etc. For completeness, we should also include parameters introduced by the numerical scheme, such as the spatial resolution, the time step, the artificial viscosity Π_AV, etc. The extent of this parameter space prohibits a thorough analysis of every collision as a function of all the mentioned parameters, we thus have to restrict ourselves to a particular set of simulations, varying some parameters and keeping the others constant.

Asteroidal targets have been considered non-rotating in most previous studies of asteroid families (Durda et al. 2007; Benavidez et al. 2012; Ševeček et al. 2017; Benavidez et al. 2018; Jutzi et al. 2019). Jutzi et al. (2013) considered a rotating Vesta, though rotating bodies have also been studied by Jutzi & Benz (2017). Ćuk & Stewart (2012) and Canup (2008) take the rotation into account for simulations of the Moon-forming impact and Canup (2005) for the impact event forming Pluto and Charon. Ballouz et al. (2015) used an N-body code to simulate collisions of rotating rubble-pile asteroids, while Takeda & Ohtsuki (2007, 2009) studied the angular momentum transfer for both stationary and rotating rubble-piles. In this work, we study formation of asteroid families from monolithic targets, extending the parameters of the simulation by the initial rotational period P_pb of the target, including cases close to the spin barrier.

The paper is organized as follows. In Sect. 2, we describe our numerical model. Section 3 analyzes differences between synthetic families created from parent bodies with various rotational periods. Section 4 is focusedon largest remnants, specifically on their accelerations or decelerations and the angular momentum transfer. Finally, we summarize our results in Sect. 5.

2 Numerical model

We developed a new SPH and N-body code. The code is publicly available, see Appendix C. In this section, we do not attemptto present a thorough review of the SPH method (as e.g., Cossins 2010), but instead we summarize the exact equations used in the code, emphasizing the modifications introduced in order to properly deal with rotating bodies.

2.1 Equations

The set of hydrodynamical equations is solved with a smoothed particle hydrodynamics (Monaghan 1985). The continuum is discretized into particles comoving with the continuum, with the density profile of the particles given by a kernel W, which is a cubic spline in our case: $$W(r,h)=\frac{\sigma }{h^3}\{\begin{array}{ll}\frac{1}{4}{(2-q)}^3-{(1-q)}^3,\hfill & 0\le q<1,\hfill \\ \frac{1}{4}{(2-q)}^3,\hfill & 1\le q<2,\hfill \\ 0\hfill & q\ge 2,\hfill \end{array}$$(1)

where q ≡ r∕h.

Below, we denote indices of particles with Latin subscripts (usually i, j, ...), the components of vector and tensor quantities with Greek superscripts (usually α, β, ...). We also use Einstein notation to sum over components (but not for particles).

The equation of motion for ith particle reads: $$\begin{array}{lll}\frac{\text{d}{v}_i^{\alpha }}{\text{d}t}\hfill & =\hfill & {\sum }_j{m}_j\left(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2}+{\Pi }_{ij}{\delta }^{\alpha \beta }+{\zeta }_{ij}^{\alpha \beta }{f}^n\right)\frac{\partial W_{ij}}{\partial x^{\beta }}\hfill \\ \hfill & \hfill & +\nabla \Phi -{\left[\omega \times (\omega \times r_i)\right]}^{\alpha }-{\left[2\omega \times v_i\right]}^{\alpha },\hfill \end{array}$$(2)

where σ^αβ = −Pδ^αβ + S^αβ is the total stress tensor, Π is the artificial viscosity (Monaghan & Gingold 1983), ζ^αβfⁿ is the artificial stress (Monaghan 2000), Φ is the gravitational potential (including both external fields and self-gravity of particles), and ω is the angular velocity of the reference frame.

The respective terms in the Eq. (2) correspond to the stress divergence, gravitational acceleration, centrifugal force, and Coriolis force. Inertial forces are only applied if the simulation is carried out in a non-inertial reference frame, corotating with the target asteroid (see Appendix B).

We used the standard artificial viscosity with linear and quadratic velocity divergence terms and coefficients α_AV and β_AV, respectively. This term is essential for a proper shock propagation and thus was always enabled in our simulations. Optionally, the code allows enabling the Balsara switch (Balsara 1995), which reduces the artificial viscosity in shear motions in order to reduce an unphysical angular momentum transfer. Additionally, the code includes the artificial stress term ζ^αβ fⁿ, which reduces the tensile instability, meaning an unphysical clustering of particles due to negative pressure. We tested the effects of this term using the “colliding rubber cylinders” test (cf. Schäfer et al. 2016).

We used a different discretization of the equation than in Ševeček et al. (2017), as we found the above equation to be more robust, avoiding a high-frequency oscillation in the pressure field. This is a recurring problem in high-velocity impact simulations and while it can be suppressed by a larger kernel support, additional modifications of the method have been suggested to address the issue, for example the δ-SPH modification (Marrone et al. 2011). The density is evolved using the continuity equation: $$\frac{\text{d}{\rho }_i}{\text{d}t}={\sum }_j{m}_j\frac{\partial v_i^{\alpha }}{\partial x^{\alpha }}.$$(3)

We solved the evolution equation instead of direct summation of particle masses to avoid the artificial low-density layer at the surface of the asteroid (Reinhardt & Stadel 2017). The velocity gradient at the right-hand side of Eq. (3) is computed as: $${\rho }_i\frac{\partial v_i^{\alpha }}{\partial x^{\beta }}={\sum }_j{m}_j(v_j^{\alpha }-v_i^{\alpha })C_i^{\beta \gamma }\frac{\partial W_{ij}}{\partial x^{\gamma }},$$(4)

where the correction tensor C^αβ is defined as(Schäfer et al. 2016): $$C_i^{\alpha \beta }\equiv {\left[{\sum }_j\frac{m_j}{{\rho }_j}(r_j^{\alpha }-r_i^{\alpha })\frac{\partial W_{ij}}{\partial x^{\beta }}\right]}^{-1}.$$(5)

In the case the bracketed matrix was not invertible, we used the Moore–Penrose pseudo-inverse instead. The correction tensor is further set to 1 for fully damaged material.

The correction tensor has been introduced to tackle the linear inconsistency of the standard SPH formulation. It is a fundamental term in the velocity gradient that allows for a stable bulk rotation of the simulated body and significantly improves the conservation of the total angular momentum. We evolved the internal energy u using the energy equation: $$\begin{array}{lll}\frac{\text{d}{u}_i}{\text{d}t}\hfill & =\hfill & -\frac{{\sigma }^{\alpha \beta }}{{\rho }_i}\frac{\partial v_i^{\alpha }}{\partial x^{\beta }}+{\sum }_j\frac{1}{2}{m}_j{\Pi }_{ij}\left(v_i^{\alpha }-v_j^{\alpha }\right)\frac{\partial W_{ij}}{\partial x^{\alpha }}\hfill \\ \hfill & \hfill & +{\sum }_j\frac{1}{2}{m}_j{\zeta }_{ij}^{\alpha \beta }{f}^n\left(v_i^{\beta }-v_j^{\beta }\right)\frac{\partial W_{ij}}{\partial x^{\alpha }}.\hfill \end{array}$$(6)

In this equation the velocity gradient is also corrected by the tensor C^αβ. While this is required for a consistent handling of rotation, the inequality of kernel gradients used in the energy equation (Eq. (6)) and in the equation of motion (Eq. (2)) implies the total energy is generally not conserved in the simulations. This is usually not an issue, as the total energy does not increase by more than 5%.

However, in some cases (e.g., weak cratering impacts or exceedingly long integration time), the energy growth can be prohibitive. For such cases, the code also offers an alternative way to evolve the internal energy, using a compatibly-differenced scheme (Owen 2014). Instead of computing the energy derivative, the energy change is computed directly from particle pair-wise accelerations $a_{ij}^{\alpha }$ and half-step velocities $w_i^{\alpha }=v_i^{\alpha }+\frac{1}{2}{a}_i^{\alpha }\Delta t$, using the equation: $$\Delta u_i={\sum }_j{f}_{ij}\left(w_j^{\alpha }-w_j^{\beta }\right)a_{ij}^{\alpha }\Delta t,$$(7)

where the energy partitioning factors f_ij can be chosen arbitrarily, provided they fulfill constraint f_ij + f_ji = 1. With this form of SPH, the total energy can be conserved to machine precision, at a cost of performance overhead. However, this does not solve the inconsistency mentioned above.

The listed set of equations is supplemented by the Tillotson equation of state (Tillotson 1962). To close the set of equations, we have to specify the constitutive equation. We used the Hooke’s law, evolving the deviatoric stress tensor in time using: $$\frac{\text{d}{S}_i^{\alpha \beta }}{\text{d}t}=2\mu \left(\frac{\partial v_i^{\alpha }}{\partial x^{\beta }}-\frac{1}{3}{\delta }^{\alpha \beta }\frac{\partial v_i^{\gamma }}{\partial x^{\gamma }}\right),$$(8)

where μ denotes the shear modulus. To account for plasticity of the material, we further used the von Mises criterion, which reduces the deviatoric stress tensor by the factor: $$f=\min \left[\frac{Y_0^2{(1-u/u_{\text{melt}})}^2}{\frac{3}{2}{S}^{\alpha \beta }{S}^{\alpha \beta }},1\right],$$(9)

where Y₀ is the yield stress and u_melt is the specific melting energy. While more complex, pressure-dependent rheology models exist (Jutzi et al. 2015), von Mises rheology is reasonable for monolithic asteroids and still consistent with laboratory experiments (Remington et al. 2018). The effects of friction have been studied by Jutzi et al. (2015) or Kurosawa & Genda (2018) so we do not discuss such effect in this work.

Additionally, we integrated the fragmentation model to model an activation of flaws and a propagation of fractures in the material. Following Benz & Asphaug (1994), we define a scalar quantity damage D, modifying the total stress tensor as $${\sigma }_D^{\alpha \beta }=-\left(1-DH(-P)\right)P{\delta }^{\alpha \beta }+(1-D)S^{\alpha \beta },$$(10)

where H(x) is the Heaviside step function. A fully damaged material (D = 1) has no shear nor tensile strength, it only resists compressions.

Smoothing lengths of particles are evolved to balance the changes of particle concentration. We thus derived the equation directly from the continuity equation: $$\frac{\text{d}{h}_i}{\text{d}t}=\frac{h_i}{3}\frac{\partial v_i^{\alpha }}{\partial x^{\alpha }}.$$(11)

Since it is also an evolution equation, we need to specify the initial conditions for smoothing lengths: $$h=\eta {\left(\frac{V}{N}\right)}^{\frac{1}{3}},$$(12)

where V is the volume of the body, N is the number of particles in the body and η is a free non-dimensional parameter, which we set to η = 1.3. This corresponds to an average number of neighboring particles N_neigh ≃ 65.

2.2 Gravity

Beyond hydrodynamics, the code also computes accelerations of SPH particles due to self-gravity. To compute it efficiently (albeit approximately), we employed the Barnes-Hut algorithm (Barnes & Hut 1986). Instead of computing pair-wise interactions of particles, we first clustered the particles hierarchically and evaluated gravitational moments of particles in each node of the constructed tree. The accelerations were then computed by a tree-walk; if the evaluated node was distant enough, the acceleration could be approximated by evaluating the multipole moments up to the octupole order, otherwise we descended into child nodes. The precision of the method is controlled by an opening radius r_open. For an extensive description of the method, see Stadel (2001).

As our SPH particles are spherically symmetric, they can be replaced by point masses, provided they do not intersect each other (the corresponding kernel W_ij is zero). However, it is absolutely necessary to account for softening of the gravitational potential for neighboring particles. We follow Cossins (2010) by introducing a gravitational softening kernel ϕ (associated with the SPH smoothing kernel W) using the equation: $$\frac{\partial \varphi }{\partial r}=\frac{4\pi }{r}{\displaystyle {\int }_0^r{r^{\prime }}^2}W(r^{\prime })\text{d}{r}^{\prime }+\frac{h}{r^2}.$$(13)

The gravitational kernel ϕ corresponding to our M4 spline kernel W is then: $$\varphi (r,h)=\{\begin{array}{ll}\frac{2}{3}{q}^2-\frac{3}{10}{q}^4+\frac{1}{10}{q}^5-\frac{7}{5},\hfill & 0\le q<1,\hfill \\ \frac{4}{3}{q}^2-q^3+\frac{3}{10}{q}^4-\frac{1}{10}{q}^5\hfill & \hfill \\ -\frac{8}{5}+\frac{1}{15q},\hfill & 1\le q<2,\hfill \\ -\frac{1}{q},\hfill & q\ge 2,\hfill \end{array}$$(14)

where q ≡ r∕h. However, this kernel does not have compact support.

2.3 Temporal discretization

Using derivatives computed at each time step, the equations were integrated using explicit timestepping. The scheme used in this work is the predictor-corrector, however other schemes are implemented in the code, such as the leapfrog, fourth order Runge–Kutta, or Bulirsch–Stoer.

The value of the time step is determined by the Courant–Friedrichs–Lewy (CFL) criterion: $$\Delta t\le C_{\text{CFL}}{\min }_i\frac{h_i}{c_{\text{s}}},$$(15)

where h_i is the smoothing length of ith particle, c_s the local sound speed, and C_CFL is the Courant number. In our simulations, we usually used C_CFL = 0.25, as higher values can lead to instabilities in some cases. Moreover, we restricted the time step by the value-to-derivative ratio of all time-dependent quantities in the simulation to control the discretization errors. The upper limit of the time step is therefore: $$\Delta t\le C_{\text{d}}\frac{\left|x\right|+s_x}{\left|\dot{x}\right|},$$(16)

where s_x is a parameter with the same dimensions as x, assigned to each quantity in order to avoid zero time step if the quantity x is zero. Constant C_d is 0.2 for all quantities.

2.4 Equilibrium initial conditions

Setting up the initial conditions for the impact simulation is not a trivial task. It is necessary to assign a particular value to the density ρ, internal energy u, and deviatoric stress tensor S to each particle, so that the configuration is stable when the impact simulation starts and the particles do not oscillate.

This problem is not restricted to simulations with rotating targets. A proper handling of initial conditions is essential in simulations of the Moon formation, collisions of planetary embryos, etc. If neglected, the initial gravitational compression would introduce macroscopic radial oscillations in the target.

For small and stationary asteroids with D ≃ 10 km, the self-gravity is much less important, in fact it was often completely neglected during the fragmentation phase (as in Ševeček et al. 2017). These asteroids are assumed to be undifferentiated, hence it was reasonable to set up a homogeneous bulk density of ρ₀ = const. For these stationary targets, the stress tensor of particles is zero in equilibrium.

For rotating targets, however, such initial conditions are unstable due to the emerged centrifugal force (in the corotating frame). To prevent any unphysical fractures in the target, the configuration of particles has to be set up carefully, especially for asteroids rotating close to the critical spin rate. For this reason, we run a stabilization phase before the actual impact simulation, with an artificial damping of particle velocities: $$v_{\text{damp}}=\frac{v-\omega \times r}{1+\delta \Delta t}+\omega \times r,$$(17)

where v is the undamped velocity, ω the angular frequency of the target, r the position vector of the particle, δ an arbitrary damping coefficient (gradually being decreased during the stabilization phase), and Δt the actual time step. In this equation, we need to subtract and re-add the bulk rotation velocity, otherwise the damping would cause the target to slow down. We also did not integrate the fragmentation model during this phase, as the oscillations of the particles might damage the target prematurely.

While more general approaches for setting up the initial conditions exist (Reinhardt & Stadel 2017), the presented method is simple, robust, and sufficient for our purposes. A disadvantage of our method is a significant computational overhead, as for some simulations the time needed to converge into a stable solution is comparable to the duration of the actual impact simulation. However, here we performed many simulations with fixed target diameter D_pb and period P, so we had to precompute the initial conditions only once and then used the cached particle configuration for other runs.

2.5 Reaccumulation phase

As our numerical model contains both the hydrodynamics and the gravitational interactions, it could be used for the entire simulation – from the stabilization, pre-impact flight, fragmentation phase until the gravitational reaccumulation of all fragments. However, the time step is often severely limited by the CFL criterion (Eq. (15)).

We can increase the time step by several orders of magnitude and hence speed up the simulation considerably by changing the numerical scheme from SPH to N-body integration. This is a common approach in studies of asteroid families (Durda et al. 2007; Michel et al. 2015; Ševeček et al. 2017; Jutzi et al. 2019). Once the target is fully fractured and the fragments start to recede, we converted all SPH particles into hard spheres and replaced the complexity of hydrodynamic equations with a simple collision detection. This allows us to overcome the time step limitation.

We further optimized the simulation by merging the collided particles into larger spheres. By doing so, we lost information about the shape of the fragments; to preserve the shapes, it is necessary to either form rigid aggregates of particles instead of mergers (Michel & Richardson 2013), or simulate the entire reaccumulation using the SPH (as in Sugiura et al. 2018). Here, as we are mainly interested in distribution of sizes and rotational periods, merging the particles is thus a viable option.

Merging not only affects the shape, but also the dynamics of fragments. As it modifies the moment of inertia, the merger has a generally different rotational period than a real non-spherical fragment would have. Merging also removes higher gravitational moments, thus altering motion of near fragments. This is a slight limitation of the presented model.

Hard spheres are created directly from SPH particles. Their mass is unchanged, and the radius of the formed spheres is computed as: $$r_i=\sqrt[3]{\frac{3m_i}{4\pi {\rho }_i}},$$(18)

so that the volume of spheres is equal to the volume of SPH particles. As the total volume is conserved, created spheres inevitably overlap. Appendix A describes how the code handles such overlaps.

When two spheres collide, they are merged only if their relative speed is lower than the mutual escape speed $$v<v_{\text{esc}}\equiv \sqrt{\frac{2G(m_i+m_j)}{r_i+r_j}}$$(19)

and if the rotational angular frequency of the merger does not exceed the critical frequency $$\omega <{\omega }_{\text{crit}}\equiv \sqrt{\frac{G(m_i+m_j)}{r_{\text{merger}}^3}}.$$(20)

In this way we prevent formation of unphysically fast rotators. If any of these conditions is not fulfilled, particles undergo an inelastic bounce. The damping of velocities is determined by the normal η_n and tangential η_t coefficient of restitution, which we set to 0.5 and 1, respectively.

When merging the particles, we determined the mass, radius, velocity, and angular frequency of the merger, so that the total mass, volume, momentum, and angular momentum are conserved. As the tangential components of velocities are not damped by a bounce, merging is the only way to spin up fragments in our simulations.

3 Synthetic families created from rotating targets

To better understand how the rotation influences impact events, we decided to compute a matrix of simulations for various impact parameters. We ran simulations for two different target sizes, D_pb = 10 km and D_pb = 100 km, in order to ascertain the scaling of the rotational effect. We tested head-on impacts, having the impact angle ϕ_imp = 15°, the intermediate cases with ϕ_imp = 45°, and oblique impacts with ϕ_imp = 75°. We have to further distinguish prograde events, meaning impacts where the orbital velocity has the same direction as the impact velocity, and retrograde events, where the orbital velocity has the opposite direction. In the following, the prograde impacts have positive values of impact angles, while the retrograde impacts have negative values.

In all of our simulations, we set the impact velocity to v_imp= 5 km s⁻¹, which is close to the mean velocity for Main-belt collisions (Dahlgren 1998). The simulation matrix covers both the cratering and the catastrophic events. We ran simulations with relative impact energies $Q/Q_{\text{D}}^{\star }=0.1,0.3,1,\text{ and }3$, where the critical energy $Q_{\text{D}}^{\star }$ is given by the scaling law of Benz & Asphaug (1999). As $Q_{\text{D}}^{\star }$ is defined as the specific impact energy (relative to mass of the target) required to eject 50% of the target’s mass as fragments, it necessarily depends on the rotational period of the target. However, we consider $Q_{\text{D}}^{\star }$ to be independentof rotation and used the same value for all performed simulations, as it provides a convenient dimensionless measure of the impact energy.

We assume that both the target and the impactor are monolithic bodies, the material parameters are summarized in Table 1. The spatial resolution of the target was approximately N = 500 000 SPH particles and the number of projectile particles was selected to match the particle density of the target. Three simulations with different periods P_pb are shown in Fig. 1.

Fig. 1 Impacts into Dpb = 100 km targets without rotation (first row), rotational period Ppb = 3 h (second row), and Ppb = 2 h (third row). The period of two hours is approximately the critical period of the target. The color brightness corresponds to specific internal energy of the particles.

3.1 Coordinate system

Due to the rotation, the impact geometry is more complex compared to the stationary case, where it was determined by a single parameter – the impact angle ϕ_imp between the normal at the impact point and the velocity vector of the impactor. To describe the impact into a rotating target, we first defined a coordinate system of the simulations. We placed the target at origin with zero velocity. The impactor has velocity [−v_imp;0;0] and its position in x-y plane is given by ϕ_imp, specifically: $$r_0=\left[\begin{array}{l}{x}_0+0.5(D_{\text{pb}}+d_{\text{imp}})\cos {\varphi }_{\text{imp}},\hfill \\ 0.5(D_{\text{pb}}+d_{\text{imp}})\sin {\varphi }_{\text{imp}},\hfill \\ 0\hfill \end{array}\right],$$

where x₀ is the distance of the impactor from the impact point. These initial conditions have a mirror symmetry in z.

The rotation vector ω_pb of the target adds an additional three free parameters into the simulation setup. Generally, the vector does not have to be aligned with any coordinate axis. We reduced the number of free parameters and thus simplified the analysis by aligning the vector with z-axis, meaning we only consider impacts in the equatorial plane of the target. We expect these impacts will be affected by the rotation the most, as the centrifugal force is the largest on the equator. Furthermore, the angular momentum of the target is aligned with the angular momentum of the impactor. We can thus expect the largest changes in the angular momentum. These expectations have been confirmed for rubble-pile bodies by N-body simulations of Takeda & Ohtsuki (2009).

3.2 Size-frequency distributions for D_pb = 10 km targets

The first set of simulations was carried out with the target size D_pb = 10 km. The diameters of impactors were d_imp = 394, 570, 850, and 1226 m, respectively. We ran a number of simulations for different P_pb, d_imp, and ϕ_imp and compared the size-frequency distribution (SFD) of a family created by an impact into a rotating target with corresponding impact into a stationary target. The resulting distributions are plotted in Fig. 2.

At firstglance, the differences between the targets rotating with a period of P_pb = 3 h and the non-rotating targets are relatively small. The slope of the SFD is almost unchanged in most simulations, it is only shifted as more mass is ejected from the rotating target. In several simulations, like for ϕ_imp = −45° and d_imp = 0.85 km, we can see a larger number of fragments in the middle part of the SFD for rotating targets; fragments that would reaccumulate to the largest remnant in the stationary case now escape due to the extra speed from rotation and contribute to the family.

Much larger differences in SFDs can be seen for the target with period P_pb = 2 h, which is rather expected; for ρ₀ = 2700 kg m⁻³ the critical period is P_crit ≃ 2.009 h, so a P_pb = 2 h target actually rotates very slightly supercritically, although it is held stable by the material strength. The difference is the most prominent for oblique ϕ_imp = ±75° impacts. In several cases, the rotation seems to make the SFD less steep, although this might be partially attributed to a numerical artifact, as the synthetic families of non-rotating targets are already very close to the resolution limit.

On the other hand, energetic impacts produce practically the same SFDs regardless of P_pb. In this regime, angular momentum of projectiles is larger than the rotational angular momentum of the target; for $Q/Q_{\text{D}}^{\star }=3$ and P_pb = 2 h, the angular momentum of a projectile is five times larger. Ejection velocities are also considerably larger than orbital velocities, hence it is not surprising that the rotation does not make a substantial difference.

It is not probable that these differences come from different fragmentation patterns, as targets are fully damaged by the impact. In our model, such damaged material is strengthless and it essentially behaves like a fluid. Since there is no internal friction nor a mechanism to regain the material strength, this model is insufficient to determine shapes of the fragments; however, here we are only interested in size distributions and using a simplified model is therefore justified.

Fig. 2 Cumulative size-frequency distributions N(>D) of synthetic families for Dpb = 10 km targets. Stationary targets are plotted in black; red and yellow plots correspond to targets with rotational period P = 3 h and P = 2 h, respectively.Columns 3–5 of the plot show prograde impacts (i.e., positive impact angles); Cols. 1 and 2 are retrograde impacts (i.e., negative impact angles).

3.3 Size-frequency distributions for D_pb = 100 km targets

It is a priori not clear how rotation affects targets of different sizes. To preliminarily estimate the importance of initial rotation, we computed the ratio of the angular frequency ω_pb of the targetand ω_imp of the impactor with respect to the target: $$\frac{{\omega }_{\text{pb}}}{{\omega }_{\text{imp}}}~\frac{D_{\text{pb}}}{v_{\text{imp}}{P}_{\text{pb}}\sin {\varphi }_{\text{imp}}}.$$(21)

Because the ratio scales linearly with the target size D_pb, we expect that the rotation will play a bigger role for impacts into larger targets; however, this back-of-the-envelope computation is by no means definite proof and it needs to be tested.

To thispoint, we ran a set of simulations with target size D_pb = 100 km. The set is analogous to the one in Sect. 3.2: we used the same impact angles and rotational periods, the impactor diameters were d_imp = 11.170, 16.110, 24.064, and 34.706 km in order to obtain the required relative energies $Q/Q_{\text{D}}^{\star }$. The size-frequency distributions of the synthetic families are plotted in Fig. 3.

As expected, the differences between rotating and non-rotating targets are indeed substantially larger than for D_pb = 10 km. The rotation can completely change the impact regime from cratering to catastrophic; see for example the impact with ϕ_imp = 15° and d_imp = 16.110 km, where a cratering gradually changes to a catastrophic disruption as we decrease P_pb. Focusing on P_pb = 3 h targets, they produce very shallow SFD in case of oblique prograde impacts. For ϕ_imp = ±45° cratering impacts, we see numerous intermediate-sized fragments (somewhat separated in the SFD) if the target is non-rotating, but the SFD becomes continuous when rotation is introduced.

The effects are even stronger for critically rotating bodies with P_pb = 2 h, of course. Generally, SFDs of formally cratering events are more similar to catastrophic ones. It also seems that oblique retrograde craterings produce more fragments than prograde ones. For the impact ϕ_imp = 15° and d_imp = 24.064 km, the SFD is well below the non-rotating case and most of the mass is contained in the smallest fragments.

Although large (D ≫ 10 km) asteroids typically rotate much slower than smaller bodies, there are a few that rotate close to the critical spin rate for elongated bodies, such as (216) Kleopatra (Hirabayashi & Scheeres 2014). Rotation in collisional simulations of such bodies should therefore not be neglected.

3.4 Total ejected mass

While Figs. 2 and 3 clearly show the differences between the SFDs, it is quite difficult to read the total mass ejected from the target during the impact. Even when the SFDs of a rotating and a stationarytarget seem to differ only negligibly, the total integrated mass of fragments may be significantly different.

To show the effect of the initial rotation on the ejected mass clearly, we performed over 400 simulations with the target size D_pb = 10 km, various impact angles ϕ_imp, projectile diameters d_imp, and initial periods P_pb of the target. These simulations have a lower spatial resolution compared to the simulations of family formation in previous sections, as here we do not need to resolve individual fragments in detail. The target is resolved by approximately N = 100 000 particles.

We ran simulations for ϕ_imp ranging from 15°–75° (both prograde and retrograde). To capture the dependence on P_pb, we selected nine different values from P_pb = P_crit to 50 P_crit. The impact energies of the simulations were $Q/Q_{\text{D}}^{\star }=0.03,0.1,0.3,\text{ and }1$, meaning thesimulations range from cratering events to mid-energy events.

Our goal is to compute the total mass of the fragments as a function of the impact angle ϕ_imp, the initial rotational period P_pb and the diameter d_imp. We are actually not interested in the absolute value of the ejected mass, but rather in the ejected mass relative to the mass that would be ejected if the targets were stationary. Therefore, we computed the ratio: $${\mu }_{\text{ej}}=\frac{M_{\text{ej}}({\varphi }_{\text{imp}},d_{\text{imp}},P_{\text{pb}})}{M_{\text{ej}}({\varphi }_{\text{imp}},d_{\text{imp}},\infty )}$$(22)

and plot the result in Fig. 4. Values μ_ej < 1 would mean that the impact into the rotating target ejected fewer fragments, compared to the stationary target; no such result was found in the performed simulations.

Generally, the rotation amplifies the ejection by several tens of percent. However, the increase is significantly higher if the following conditions are satisfied:

• Target rotates near the critical period. As expected, the effect of rotation decreases rapidly with the increasing period of the target.

• The impact results in a cratering rather than a catastrophic event. While high-energy impacts eject more fragments in an absolute measure, the initial rotation does not affect the value notably in this regime.

• The impact is oblique and has a prograde direction. Head-on impacts and the impacts in retrograde directions are not affected by the rotation to the same degree.

In extreme cases, the rotation can amplify the ejected mass by a factor of five. On the other hand, the ejection ratio μ_ej does not exceed 1.6 for rotational periods P_pb > 2P_crit in any of the performed simulations.

Although it is a different rheology, rubble-pile bodies also exhibit a minor effect of rotation (on $Q_{\text{D}}^{\star }$ as well as μ_ej) in this range of P_pb (Takeda & Ohtsuki 2009, see Fig. 2 therein). However, their strength is an order of magnitude lower than for monoliths of Benz & Asphaug (1999), so the comparison is not straightforward.

Fig. 3 Cumulative size-frequency distribution for Dpb = 100 km targets. Thenotation is the same as in Fig. 2.

Fig. 4 Total mass of fragments Mej(P)∕Mej(∞) ejected by collisions, normalized by mass of fragments from corresponding collision into stationary target. The four figures correspond to different relative impact energies, from cratering (left) to mid-energy (right) events. The ejected mass is plotted as a function of the impact angle ϕimp and initial period Ppb of the target.

Fig. 5 Damage D of the target at time t = 10 s after impact for various impact angles, from left to right, ϕimp = 75°, 45°, 15°, −45°, −75°. Simulations were carried out with the impactor of size dimp = 314 m and speed vimp = 5 km s−1; target was not rotating. The red outline shows the original position of the target and the impactor. There is an undamaged cavity only for oblique impacts, otherwise the target is fully damaged by the impact.

Fig. 6 Change of spin rate Δω of target (or largest remnant in case $Q/Q_{\text{D}}^{\star }\simeq 1$) as function of initial period Ppb of target (here in units of critical period Pcrit) and impact angle ϕimp. The four images correspond to different sizes of the impactor. The energy of the impact is from left to right $Q/Q_{\text{D}}^{\star }=0.03,0.1,0.3\text{ and }1$. The impact velocity was vimp = 5 km s−1 in all cases.

4 Embedding and draining the angular momentum

Impact into a rotating target can cause either an acceleration or a deceleration of the target’s rotation. This can be immediatelyseen from two limit cases: a stationary target is always spun by the impact, on the other hand a target rotating at the breakup limit cannot be accelerated any further and the collision thus always causes a deceleration.

It has been proposed that rotating asteroids are decelerated over time by numerous subsequent cratering collisions, as a fraction of the angular momentum is carried away by fragments. Coined the angular momentum drain (Dobrovolskis & Burns 1984), this effect could explain the excess of slow rotators in the main belt. In this section, we examine whether this effect emerges in our simulations and we determine the functional dependence of the deceleration on the impact parameters.

We analyze the angular momentum transfer as a function of impact parameters, using the set of simulations described in Sect. 3.4. The impacts range from cratering ($Q/Q_{\text{D}}^{\star }~0.03$) to mid-energy ($Q/Q_{\text{D}}^{\star }~1$) events. For $Q/Q_{\text{D}}^{\star }~1$, the whole target asteroid is disintegrated by the collision and fragments with mass of about 0.5 M_pb are reaccumulated later, forming the largest remnant. This can no longer be viewed as a cratering event that merely modifies the rotational state of the target, nevertheless we can still formally compute the relation between the period of the target and the largest remnant.

In a majority of performed simulations, the target is completely damaged by the impact. There remained an undamaged cavity only for the weakest oblique impacts with $Q/Q_{\text{D}}^{\star }=0.03$, as shown inFig. 5. Otherwise all particles of the target have damage D = 1 after the fragmentation phase.

The change of spin rate Δω of the target is plotted in Fig. 6. We plot the change of frequency rather than period, as the period is formally infinite for anon-rotating body; a change of period is thus not a meaningful quantity.

For cratering events, the prograde events (denoted with positive ϕ_imp) mostly accelerate the target, while retrograde events cause deceleration. The two exceptions from this rule are:

• 1.

  a prograde impact into a critically rotating body, in which case it cannot be accelerated any more and some deceleration is expected; and

• 2.

  a retrograde impact to an almost stationary body.

Impacts with higher energies show a different pattern. It seems that the two regions described above expand. Prograde impacts into fast rotators actually decelerate the target, while retrograde impacts start to accelerate it. For the most energetic event we studied, $Q/Q_{\text{D}}^{\star }\simeq 1$, it seems that the two regions completely swapped: most prograde impacts cause a deceleration while retrograde impacts cause an acceleration.

The despinning (or angular momentum draining) on rubble-piles in catastrophic disruptions was confirmed by Takeda & Ohtsuki (2009). In our simulations, the pattern is more complex, likely because the parameter space ($Q/Q_{\text{D}}^{\star }$, P_pb) is significantly more extended; additionally rubble-piles cannot initially rotate critically.

Fig. 7 Dimensionless effectivity γ of angular momentum transfer. The other quantities are the same as in Fig. 6.

4.1 Effectivity of the angular momentum transfer

Let us define the effectivity γ of the angular momentum transfer as: $$\gamma \equiv \frac{L_{\text{lr}}-L_{\text{pb}}}{L_{\text{imp}}},$$(23)

where L_pb is the rotational angular momentum of the target before the impact, L_pb is the rotational angular momentum of the largest remnant, and L_imp is the angular momentum of the impactor with respect to the target. As these values are scalars, we assign a negative sign to the value L_imp for retrograde impacts.

We emphasize that the effectivity γ is not necessarily in the unit interval ⟨0;1⟩. Specifically, it may be significantly larger than one for head-on impacts, as the delivered angular momentum is very low; in fact, L_imp approaches zero for ϕ_imp = 0. The effectivity can also be a negative number for impacts to critically rotating targets, as in these cases, the target cannot accelerate over the breakup limit, so a zero or even negative values of γ are expected.

The effectivity γ as a function of the initial period P_pb, the impact angle ϕ_imp, and the impactor diameter d_imp is plotted in Fig. 7. We can see that cratering impacts have generally higher effectivity than high-energy impacts. This result might have been expected, as the cratering impacts eject less mass and thus transfer less angular momentum to fragments, compared to the catastrophic impacts. A less expected outcome is the negative effectivity for the high-energy impacts. We predicted the negative values only for prograde impacts into critically rotating targets, but for d_imp = 850 m the effectivity is negative for the majority of performed simulations.

Finally, the highest effectivity is achieved for retrograde impacts into critically rotating targets. However, this results is slightly false, because in this regime the target is always decelerated; γ can therefore exceed one as the angular momentum lost in the collision is higher than L_imp (in absolute value). Impacts into slower rotators have values of γ around 0.5.

4.2 Angle-averaged ejection and momentum transfer

To express an overall effect of collisions on a rotational state of a target, it is useful to consider a large number of collisions at random impact angles and compute the average change of spin rate: $$\overline{\Delta \omega }\equiv {\displaystyle {\int }_0^{\pi /2}\Delta }\omega \sin 2\varphi \text{d}\varphi \approx \frac{{\sum }_i\Delta {\omega }_i\sin 2{\varphi }_i}{{\sum }_i\sin 2{\varphi }_i},$$(24)

where sin2ϕ is the probability of the impact with the impact angle ϕ. The averaged change $\overline{\Delta }\omega$ is plotted in Fig. 8, together with a plot of the relative mass ejection ${\overline{\mu }}_{\text{ej}}$, defined in Sect. 3.4.

The figures show that targets are on average decelerated for both low-energy and mid-energy impacts. A target is only accelerated if it was originally a very slow rotator, since it cannot be decelerated any more. The transition between these two regimes (plotted in white) is where the target does not change its angular frequency upon the impact. It seems to depend on the energy of the impact; for $Q/Q_{\text{D}}^{\star }~0.1$, the transition occurs at around P ~ 20P_crit, while for $Q/Q_{\text{D}}^{\star }~1$, it is shifted to around P ~ 6P_crit.

5 Conclusions and future work

In this paper, we showed that a fast initial rotation of targets may significantly affect resulting synthetic families. The effect is more prominent for larger target bodies and for oblique impact angles. Generally, more fragments are ejected from prograde (ϕ_imp > 0) compared to the retrograde (ϕ_imp < 0) targets.

In extreme cases, the mass ejection can be amplified by a factor five. Neglecting the rotation would therefore introduce a considerablebias. Other parameters of the simulation do introduce similar (or sometimes larger) uncertainties, for example the Weibull parameters of the fragmentation model (Ševeček et al. 2017), the rheological model of the target, etc. As shown in Jutzi &Benz (2017), the initial shape can also have a significant effect.

Throughout this paper, we assumed that both the targets and impactors are monolithic bodies. It is a priori not clear whether the rotation would be more important for rubble-pile bodies (with macro-porosity), or when a rheological model with crushing (microporosity) is used in the simulations (Jutzi et al. 2019). It should also be explored how initial shapes relate to spin rates of fragments. We postpone such studies to future works.

In the future, we also plan to determine the scaling law as a function of both D_pb and P_pb. It is clear that the critical energy $Q_{\text{D}}^{\star }$ is a steep function of P_pb close to the critical spin rate. Finding a functional dependence $Q_{\text{D}}^{\star }=Q_{\text{D}}^{\star }(P_{\text{pb}})$ might be a valuable result for studies of main belt evolution, as it could be used to construct a combined model that includes both collisions and rotations.

Fig. 8 Quantities averaged over impact angles using Eq. (24) as function of period Ppb of target andimpactor diameter dimp. The left figure shows the change of spin rate Δω, the right one the total mass of fragments, normalized by the mass of fragments from corresponding collision into a stationary target.

Acknowledgements

The work of P.Š. and M.B. has been supported by the Czech Science Foundation (GACR 18-04514J) and Charles University (GAUK 1584517). M.J. acknowledges support from the Swiss National Centre of Competence in Research PlanetS.

Appendix A Handling particle overlaps

Since we treat particles as solid spheres during the reaccumulation phase, particle overlaps are unavoidable and need to be handled by our N-body integrator. There are two main reasons why particle overlaps occur.

First, the spheres overlap initially after the hand-off (Eq. (18)). In SPH, the particles naturally overlap as they describe a continuum rather than point masses. After converting them to solid spheres, particles belonging to the same body will necessarily overlap, unless their radius is decreased significantly.

Second, overlaps occur when particles are being merged. When two spherical particles collide, they merge into a larger particle with volume equal to the sum of volumes of the colliders. This merging is an atomic operation, particles are converted into the merger in an instant rather than over several timesteps, so any other particles located close to the colliders potentially overlap the particle merger.

Our code allows for several options to resolve overlaps. One straightforward solution is to always merge the overlapping particles. While this is a simple and robust solution, it can potentially create unphysical, supercritically rotating bodies. Alternatively, we can repel the overlapping particles, so that they are in contact rather than overlap. However, this causes an “inflation” of the largest remnant after the hand-off. Even worse, the angular momentum is no longer conserved.

Another option is to abandon the 1:1 conversion of spheres and instead construct a new set of spheres inside the alpha-shape of the largest remnant (Ballouz et al. 2018). Such an approach allows the placement of spheres onto a regular grid and to thus avoid overlaps by construction. However, it is more suitable when collided particles form rigid aggregates instead of mergers. As spheres never fill the entire volume (a filling factor of hexagonal close packing is about 0.74) and the merging conserves volume, fragments would shrink considerably.

We decided to merge particles only if the spin rate of the would-be merger is lower than the critical spin rate, otherwise we allow particles to pass through each other. Of course, such handling is only applied to resolve overlaps; particles that collide are always treated as solid.

Appendix B Comparison of inertial and co-rotating reference frames

We can choose two different approaches to implementing the rotation of the target:

• 1.

  rotate the particles around the center of the target, and

• 2.

  perform the simulation in the coordinate system co-rotating with the target.

From a numerical point of view, the second approach is easier to handle, as the particles of the target initially have zero velocities and we thusavoid numerical problems with the bulk rotation outlined in Sect. 2.1. The rotation is taken into account by introducing inertial accelerations, meaning the last two terms of Eq. (2).

However, it only solves the issue partially. Even though the target is stationary (in the co-rotating frame) before the impact, the projectile can spin up the target and the impact can also eject rotating fragments. To properly handle rotating bodies in SPH, it is necessary to introduce the correction tensor in Eq. (5). This allows us to perform simulations in the inertial frame, which is a natural choice.

Ideally, these two approaches should produce identical results. To test it, we ran two simulations with D_pb = 200 km parent bodies rotating critically. The results are plotted in Fig. B.1. We observe some stochastic differences, but the spatial distribution of the ejected fragments is similar in both simulations. This test confirms the consistency of both approaches.

Fig. B.1 Impact into Dpb = 200 km target. This was computed with a rotating target (upper row) and with a rotating coordinate frame in which the target is stationary (lower row), and plotted at times t = 40, 80, 120, and 160 min after the impact. The color scale represents the specific energy of the particles (in SI units).

Appendix C Implementation notes

The code OpenSPH used throughout this work is open-source, available under MIT license. As of August 7, 2019, it can be downloaded¹. It includes ahelpful graphic interface, allowing the user to interactively visualize the simulation (see Fig. C.1). It is also a stand-alone viewer of OpenSPH output files and potentially files generated by other particle-based codes, provided their file formats are implemented.

Fig. C.1 Screenshot of graphic interface of our code.

Our code can be used as both an SPH solver and an N-body integrator, as we separated computations of SPH derivatives and gravitational accelerations. In each time step, accelerations due to hydrodynamics and gravity are computed independently and summed up.Even though we miss some optimization possibilities with this approach, it allows us to use the same code for both the fragmentation and the reaccumulation phase; the hand-off is thus only an internal change of a solver, replacing SPH hydrodynamics with a collision detection.

We used the Barnes–Hut algorithm to evaluate gravitational accelerations (Barnes & Hut 1986). The code uses the same functionality in SPH and N-body solver. It only differs in the softening kernel ϕ; in SPH, it is defined by Eq. (14), while for N-body it corresponds to a homogeneous sphere. Our implementation uses a k-d tree, which is also used to find particle neighbors.

In the fragmentation phase, we found that the time step criterion that uses stress tensor derivatives (see Eq. (16)) is often unnecessarily restricting, as the stress tensor changes rapidly inside the projectile. However, as we are not very interested in remnants of projectiles, simulations can be thus be sped up by applying the criterion only for particles of the target. No such optimization is applied for CFL criterion, as it determines stability of the method and it is thus essential for all particles.

In the reaccumulation phase, collided particles are merged, provided their relative velocity does not exceed the escape velocity v_esc (as in Eq. (19)) and the spin rate does not exceed the critical spin rate ω_crit (as in Eq. (20)). In practice, both v_esc and ω_crit are multiplied by user-defined factors (i.e., a merging limit). It may be useful to tune it in such a way that a simplified N-body model of reaccumulation matches a full SPH simulation (in terms of resulting SFDs).

The total computation time needed depends on a type of simulation. Generally, catastrophic impacts take longer to compute than cratering impacts (quantities change more rapidly, hence smaller time steps are needed).For smaller targets, SPH particles are also smaller, which in turn implies smaller time steps due to the CFL criterion. The computation time is thus longer. The code is parallelized using a custom thread pool, utilizing the native C++11 threads, or optionally using the Intel Thread Building Blocks library. For D_km = 10 km target and N = 500 000 particles, a single simulation takes about ten hours on a AMD Ryzen Threadripper 1950X 16-Core processor.

References

All Tables

All Figures

	Fig. 1 Impacts into Dpb = 100 km targets without rotation (first row), rotational period Ppb = 3 h (second row), and Ppb = 2 h (third row). The period of two hours is approximately the critical period of the target. The color brightness corresponds to specific internal energy of the particles.
In the text

	Fig. 2 Cumulative size-frequency distributions N(>D) of synthetic families for Dpb = 10 km targets. Stationary targets are plotted in black; red and yellow plots correspond to targets with rotational period P = 3 h and P = 2 h, respectively.Columns 3–5 of the plot show prograde impacts (i.e., positive impact angles); Cols. 1 and 2 are retrograde impacts (i.e., negative impact angles).
In the text

	Fig. 3 Cumulative size-frequency distribution for Dpb = 100 km targets. Thenotation is the same as in Fig. 2.
In the text

	Fig. 4 Total mass of fragments Mej(P)∕Mej(∞) ejected by collisions, normalized by mass of fragments from corresponding collision into stationary target. The four figures correspond to different relative impact energies, from cratering (left) to mid-energy (right) events. The ejected mass is plotted as a function of the impact angle ϕimp and initial period Ppb of the target.
In the text

Fig. 5 Damage D of the target at time t = 10 s after impact for various impact angles, from left to right, ϕimp = 75°, 45°, 15°, −45°, −75°. Simulations were carried out with the impactor of size dimp = 314 m and speed vimp = 5 km s−1; target was not rotating. The red outline shows the original position of the target and the impactor. There is an undamaged cavity only for oblique impacts, otherwise the target is fully damaged by the impact.

In the text

Fig. 6 Change of spin rate Δω of target (or largest remnant in case $Q/Q_{\text{D}}^{\star }\simeq 1$) as function of initial period Ppb of target (here in units of critical period Pcrit) and impact angle ϕimp. The four images correspond to different sizes of the impactor. The energy of the impact is from left to right $Q/Q_{\text{D}}^{\star }=0.03,0.1,0.3\text{ and }1$. The impact velocity was vimp = 5 km s−1 in all cases.

In the text

	Fig. 7 Dimensionless effectivity γ of angular momentum transfer. The other quantities are the same as in Fig. 6.
In the text

	Fig. 8 Quantities averaged over impact angles using Eq. (24) as function of period Ppb of target andimpactor diameter dimp. The left figure shows the change of spin rate Δω, the right one the total mass of fragments, normalized by the mass of fragments from corresponding collision into a stationary target.
In the text

	Fig. B.1 Impact into Dpb = 200 km target. This was computed with a rotating target (upper row) and with a rotating coordinate frame in which the target is stationary (lower row), and plotted at times t = 40, 80, 120, and 160 min after the impact. The color scale represents the specific energy of the particles (in SI units).
In the text

	Fig. C.1 Screenshot of graphic interface of our code.
In the text