Collisions between Gravity-Dominated Bodies: 1. Outcome Regimes and Scaling Laws
Collisions are the core agent of planet formation. In this work, we derive an analytic description of the dynamical outcome for any collision between gravity-dominated bodies. We conduct high-resolution simulations of collisions between planetesimals; the results are used to isolate the effects of different impact parameters on collision outcome. During growth from planetesimals to planets, collision outcomes span multiple regimes: cratering, merging, disruption, super-catastrophic disruption, and hit-and-run events. We derive equations (scaling laws) to demarcate the transition between collision regimes and to describe the size and velocity distributions of the post-collision bodies. The scaling laws are used to calculate maps of collision outcomes as a function of mass ratio, impact angle, and impact velocity, and we discuss the implications of the probability of each collision regime during planet formation.
Collision outcomes are described in terms of the impact conditions and the catastrophic disruption criteria, – the specific energy required to disperse half the total colliding mass. All planet formation and collisional evolution studies have assumed that catastrophic disruption follows pure energy scaling; however, we find that catastrophic disruption follows nearly pure momentum scaling. As a result, is strongly dependent on the impact velocity and projectile-to-target mass ratio in addition to the total mass and impact angle. To account for the impact angle, we derive the interacting mass fraction of the projectile; the outcome of a collision is dependent on the kinetic energy of the interacting mass rather than the kinetic energy of the total mass. We also introduce a new material parameter, , that defines the catastrophic disruption criteria between equal-mass bodies in units of the specific gravitational binding energy. For a diverse range of planetesimal compositions and internal structures, has a value of ; whereas for strengthless planets, we find . We refer to the catastrophic disruption criteria for equal-mass bodies as the principal disruption curve, which is used as the reference value in the calculation of for any collision scenario. The analytic collision model presented in this work will significantly improve the physics of collisions in numerical simulations of planet formation and collisional evolution.
Subject headings:planets and satellites: formation, methods: numerical
Planet formation is common and the number and diversity of planets found increases almost daily (e.g., Borucki et al., 2011; Howard et al., 2011). As a result, planet formation theory is a rapidly evolving area of research. At present, observations principally provide snapshots of either early protoplanetary disks or stable planetary systems. Little direct information is available to connect these two stages of planet formation, therefore, numerical simulations are used to infer the details of possible intermediate stages. However, the diversity of extrasolar planetary systems continues to surprise observers and theorists alike.
A complete model of planet formation has eluded the astrophysics community because of both incomplete physics in numerical simulations and computational constraints. In order to make the problem of planet formation more tractable, the process is often divided into separate stages, which are then tackled in isolation. This method has had some success. For example, -body simulations show that large ( km) planetesimals may grow into protoplanets of about a lunar mass on million year time scales (e.g., Kokubo & Ida, 2002). Other simulations, focusing on later stages of planet formation, created a variety of stable planetary systems from initial distributions of protoplanet size bodies (e.g., Chambers, 2001; Agnor et al., 1999). Recently, the distribution of stable planets has been investigated in population synthesis models (e.g., Mordasini et al., 2009; Ida & Lin, 2010; Schlaufman et al., 2010; Alibert et al., 2011). However, given the complexity of planet formation, it is unsurprising that the predictions from the first population synthesis models have been overturned by the rapidly growing catalog of exoplanets (Howard et al., 2010, 2011). Hence, the diversity of the extrasolar planets is still unexplained.
At the heart of the standard core-accretion model of planet formation is the growth of planetestimals. The evolution of planetesimals is dominated by a series of individual collisions with other planetesimals (e.g., Beauge & Aarseth, 1990; Lissauer, 1993). The outcome of each collision depends on the specific impact conditions: target size, projectile size, impact parameter, impact velocity, and some internal properties of the target and projectile, such as composition and strength. In the past, direct global simulations of planetesimal evolution have assumed very simplified collision models. In -body simulations, terrestrial planet embryos were shown to grow easily from an annulus of large planetesimals if the only outcome of collisions is merging (e.g., Kokubo & Ida, 2002). However, the computational demands of such numerical methods did not permit for the tracking of the very large numbers of bodies necessary to be able to include direct calculations of the erosion of planetesimals.
Statistical methods are required to describe the full population of bodies from dust size to planets. For example, Kenyon & Bromley (2009) conducted simulations that included fragmentation but still relied upon a simple collision model. Specifically, their simulations did not fully account for the effects of the mass ratio, impact velocity, or impact angle on the collision outcome. In order to overcome these simplifications some previous studies have employed a multi-scale approach that includes direct simulations of collision outcomes within a top-level simulation of planet growth (Leinhardt & Richardson, 2005; Leinhardt et al., 2009; Genda et al., 2011a). However, multi-scale calculations significantly increase the computational requirements. In addition, the numerical methods employed in previous studies were only valid for a specific impact velocity regime. In the case of Leinhardt & Richardson (2005) and Leinhardt et al. (2009), the collision model assumed subsonic collisions and could not be extended past oligarchic growth. In the case of Genda et al. (2011a), the technique assumed strengthless bodies and cannot be used in the early phases of planetesimal growth.
A general description of collision outcomes that spans the growth from dust to planets is required to build a self-consistent model for planet formation. In previous work, the description of collision outcomes drew upon a combination of laboratory experiments and limited numerical simulations of collisions between two planetary-scale bodies (see review by Holsapple et al., 2002). Collision outcomes themselves are quite diverse, and several distinct collision regimes are encountered during planet formation: cratering, merging/accretion, fragmentation/erosion, and hit-and-run encounters.
Individual collision regimes have been described in quite varying detail. In the laboratory, the erosive regimes (cratering and disruption) have been studied most comprehensively (Holsapple, 1993; Holsapple et al., 2002); however, even these regimes lack a complete description of the dependence on all impact parameters (particularly mass ratio and impact angle). Recently, numerical studies of collisions between self-gravitating bodies of similar size have identified new types of collision outcomes including hit-and-run and mantle-stripping events (Agnor & Asphaug, 2004a; Asphaug et al., 2006; Marcus et al., 2009, 2010b; Leinhardt et al., 2010; Asphaug, 2010; Kokubo & Genda, 2010; Benz et al., 2007; Genda et al., 2011b). Up to this point our understanding of these new regimes has not been sufficient to implement the diversity of collision outcomes in planet formation codes. In addition, the transitions between regimes are not clearly demarcated in the literature.
In the work reported here, we present a complete description of collision outcomes for gravity-dominated bodies. Using a combination of published hydrocode and new and published -body gravity code simulation results, we derive analytic equations to demarcate the transitions between collision regimes and the size and velocity distribution of the post-collision bodies. We describe how these scaling laws can be used to increase the accuracy of numerical simulations of collisional evolution without sacrificing efficiency. In a companion paper (Stewart & Leinhardt, 2011), we apply these scaling laws to the end stage of terrestrial planet formation by analyzing the range of collision outcomes from recent -body simulations.
This paper is organized as follows: §2 summarizes the numerical method for the new -body simulations. §3 derives a general catastrophic disruption scaling law. Then, we develop general scaling laws for the size and velocity distribution of fragments in the disruption regime. §4 defines the super-catastrophic and hit-and-run regimes. §5 presents the transition boundaries between collision outcome regimes from our numerical simulations and our analytic model. §6 discusses the range of applicability of our results, areas needing future work, and the implications of the scaling laws on aspects of planet formation. The Appendix summarizes the implementation of the scaling laws in numerical simulations of planet formation and collisional evolution. Table A.1 presents the definitions of variables and annotations used in this work.
2. Numerical Method
In this section, we describe the numerical method used in the new impact simulations presented in this work. Simulations of relatively slow subsonic impacts were conducted using a -body code with finite-sized spherical particles, PKDGRAV (Stadel, 2001), which has been extensively used to study the dynamics of collisions between small bodies (e.g., Leinhardt et al., 2000; Michel et al., 2001; Leinhardt & Richardson, 2002; Durda et al., 2004; Leinhardt & Richardson, 2005; Leinhardt & Stewart, 2009; Leinhardt et al., 2010).
Both the target and projectile are assumed to be rubble piles: gravitational aggregates with no bulk tensile strength (Richardson et al., 2002). The rubble pile particles are bound together purely by self-gravity. The particles themselves are indestructible and have a fixed mass and radius (for cases without merging). The equations of motion of the particles are governed by gravity and inelastic collisions. The amount of energy lost in each particle-particle collision is parameterized through the normal and tangential coefficients of restitution. The rubble piles are created by placing particles randomly in a spherical cloud and allowing the cloud to gravitationally collapse with highly inelastic particle collisions. Randomizing the internal structure of the rubble piles avoids spurious collision results due to crystalline structure of hexagonal close packing (see Leinhardt et al., 2000; Leinhardt & Richardson, 2002). The crystalline structure can cause large uncertainties in collision outcomes for super-catastrophic events.
All simulations had a target with radius of 10 km, mass of kg, bulk density of 1000 kg m, and escape velocity of 7.5 m s. The current study includes four projectile-to-target mass ratios (), four impact angles (), and a range of impact velocities spanning merging to super-catastrophic disruption. These results for a single size target body are used to derive scaling laws for any size body in the gravity regime. The target and projectile are initially separated by the sum of their respective radii to ensure that the impact angle of the impact is unchanged from the initial trajectory.
In order to resolve the size distribution after the collisions, each body needs a relatively high number of particles (, depending on the mass ratio). However, large numbers of particles are also time consuming to integrate, especially in a rubble pile configuration where there is a high frequency of particle-particle collisions. Each simulation here uses high resolution with inelastic particle collisions to resolve the initial impact. Once the velocity field is well established, the particles are allowed to merge with one another. Thus, our method has both accuracy and efficiency, resolving the size distribution to small fragments and completing the simulations as quickly as possible.
We considered the possible influence of the time of the transition from inelastic bouncing to perfect merging. Figure 1 presents a test case using a head-on catastrophic impact between equal-sized objects (mass ratio ). Each body has particles; in the inelastic bouncing phase, each particle has a normal coefficient of restitution, , and a tangential coefficient of restitution, , consistent with field observations and friction experiments on rocky materials (e.g., Chau et al., 2002). At a certain time, the colliding particles are allowed to merge, producing one particle with the same mass and bulk density as the two original particles. If merging is turned on too early, the mass of the largest remnant is overestimated (green and cyan lines) due to a geometric effect known as runaway merging (Leinhardt & Stewart, 2009). The results of this numerical test show that the mass distribution is stable if merging is turned on after 50 steps, where 1 step is one minute in simulation time in the frame of the particles. However, we choose to be conservative and merge after 250 steps of inelastic bouncing in all of the new simulations presented in this paper. All simulations were run for at least 0.2 years, at which point the size distribution had stabilized and clumps of rubble pile fragments were easily identifiable.
Previous studies using PKDGRAV did not have the numerical resolution to determine an accurate size or velocity distribution of the collisional remnants (e.g., , Leinhardt et al., 2000). In this work, we present more extensive simulations at an order of magnitude higher resolution (). Note that is high resolution for -body simulations of colliding rubble-pile bodies with bouncing particles. We conducted several resolution tests and find that the random error on the mass of the largest remnant is a few percent of the total system mass. Hence, the super-catastrophic impacts, where the largest remnant mass is a few percent of the total system mass, have the highest error. We achieve excellent reproducibility of the slope of the size and velocity distributions with the nominal resolution compared to the higher resolution tests.
Note that -body simulations are inherently higher resolution compared to smoothed particle hydrodynamics (SPH) simulations. Our simulations resolve over a decade in fragment size, comparable to SPH simulations using an order of magnitude more particles (Durda et al., 2004). Fragments of radius 0.5 km are considered the smallest usable fragments in these simulations.
In the following sections, we also include published results of subsonic and supersonic collisions from previous work (Leinhardt & Stewart, 2009; Agnor & Asphaug, 2004a, b; Marcus et al., 2009, 2010b; Durda et al., 2004, 2007; Jutzi et al., 2010; Benz & Asphaug, 1999; Benz, 2000; Stewart & Leinhardt, 2009; Korycansky & Asphaug, 2009; Nesvorný et al., 2006; Benz et al., 2007). Studies of supersonic collisions utilize shock physics codes, which include the effects of irreversible shock deformation. For computational efficiency, the shock code is generally used to calculate only the early stages of an impact event; after a few times the shock wave crossing time, the amplitude of the shock decays to the point where further deformation is negligible. After the hydrocode step, the gravitational reaccumulation stage of disruptive events has been calculated directly using PKDGRAV (e.g., Leinhardt & Stewart, 2009; Durda et al., 2004; Michel et al., 2002) or indirectly by iteratively solving for the mass bound to the largest fragment (e.g., Benz & Asphaug, 1999; Benz, 2000; Marcus et al., 2009).
3. Results: The Disruption Regime
In our model, the boundaries between collision outcome regimes are defined using the catastrophic disruption criteria, the specific energy required to gravitationally disperse half the total mass, because it provides a convenient means of calculating the mass of the largest remnant. Our definition of the disruption regime refers to collisions in which the energy of the event results in mass loss (fragmentation) between about 10% and 90% of the total mass. More quantitatively, the disruption regime is defined as collision that result in the largest remnant having a linear dependence on the specific impact energy. The rationale for this definition will become apparent in §3.1.1.
This section focuses on deriving the dynamical outcome (the mass and velocity distribution of post-collision fragments) in the disruption regime. Other collision regimes are discussed in §4. Before discussing the results of our new numerical simulations, we briefly review the catastrophic disruption criteria, as it is a fundamental part of our story.
In the literature on planetary collisions, traditionally denotes the specific energy of the impact (kinetic energy of the projectile/target mass) and indicates the catastrophic disruption criteria, where the largest remnant has half the target mass. Upon recognition that gravitational dispersal was important, and denoted the criteria for shattering in the strength regime and dispersal in the gravity regime, respectively. All of the previous definitions for assumed that the projectile mass, , was much smaller than the target mass, ; however, in several phases of planet formation it is expected that . Therefore, in previous work, we developed a disruption criteria in the center of mass reference frame in order to study collisions between comparably sized bodies (Stewart & Leinhardt, 2009). The subscript was added in the modification of the specific energy definition to denote reduced mass. The center of mass specific impact energy is given by
where , is the reduced mass , is the impact velocity, and and are the speed of the projectile and target with respect to the center of mass, respectively. At exactly the catastrophic disruption threshold,
where we explicitly define to be the critical impact velocity required to disperse half of the total mass for a specific impact scenario (total mass and mass ratio).
The catastrophic disruption criteria is a strong function of size with two components: a strength regime where the critical specific energy decreases with increasing size and a gravity regime where the critical specific energy increases with increasing size. The transition between regimes occurs between a few 100-m and few-km radius, depending on the strength of the bodies (see Figure 2, Stewart & Leinhardt, 2009). A general formula for as a function of size was derived by Housen & Holsapple (1990) using -scaling theory,
where the first term represents the strength regime and the second the gravity regime. is the spherical radius of the combined projectile and target masses at a density of kg m. The variable was introduced by Stewart & Leinhardt (2009) in order to fit and compare the disruption criteria for collisions with different projectile-to-target mass ratios and to account for bodies with different bulk densities (e.g., rock and ice). is the gravitational constant; and are dimensionless coefficients with values near 1. is a measure of the material strength in units of Pa s, and the remaining variables, and , are dimensionless material constants. is a measure of the strain-rate dependence of the material strength with values ranging from 6 to 9 (e.g., Housen & Holsapple, 1990, 1999). is a measure of how energy and momentum from the projectile are coupled to the target; is constrained to fall between 1/3 for pure momentum scaling and 2/3 for pure energy scaling (Holsapple & Schmidt, 1987). Note that the form of equation 3 assumes that the projectile and target have the same density.
In the strength regime, the largest post-collision remnant is a mechanically intact fragment. The catastrophic disruption criteria decreases with increasing target size because more flaws grow and coalesce during the longer loading duration in larger impact events (e.g., Housen & Holsapple, 1999). In the gravity regime, disruption requires both fracturing and gravitational dispersal (Melosh & Ryan, 1997; Benz & Asphaug, 1999); hence the disruption criteria increases with increasing target size. In this regime, the largest remnant is a gravitational aggregate composed of smaller intact fragments. In both regimes, the disruption criteria increases with impact velocity because more of the impact kinetic energy is dissipated by shock deformation at higher velocities (Housen & Holsapple, 1990). This work focuses on the gravity regime; the strength regime will be the subject of future studies.
Both equations 2 and 3 are satisfied by collisions at exactly the catastrophic disruption threshold. The general formula for catastrophic disruption given by equation 3 describes a family of curves that depend on size, impact velocity, and material parameters ( and ). Most previous work fit the material parameters in equation 3 to planetary bodies of a particular composition under various assumptions (e.g., a fixed impact velocity). Next (§3.1.2–3.1.4), we present a general method to calculate the values for the disruption energy and critical impact velocity for specific impact scenarios and materials.
3.1. Derivation of a general catastrophic disruption law in the gravity regime
3.1.1 The Universal Law
In previous work, using simulations of head-on impacts, we determined the value of for a particular pair of planetary bodies by fitting the mass of the largest post-collision remnant, , as a function of the specific impact energy, . The simulations held the projectile-to-target mass ratio fixed and varied the impact velocity. For a wide range of target masses, projectile-to-target mass ratios, and critical impact velocities, we found that the mass of the largest remnant is approximated by a single linear relation,
where was fitted to be the specific energy such that (Stewart & Leinhardt, 2009; Leinhardt & Stewart, 2009). We found that a single slope agreed well with results from both laboratory experiments and numerical simulations. Furthermore, the dimensional analysis by Housen & Holsapple (1990) supports the linearity of the largest remnant mass with impact energy near the catastrophic disruption threshold. Hence, we refer to equation 4 as the “universal law” for the mass of the largest remnant.
However, the most likely collision between two planetary bodies is not a head-on collision; a impact angle is most probable (Shoemaker, 1962). The impact parameter is given by , where is the angle between the centers of the bodies and the velocity vector at the time of contact (Figure 2). The impact parameter has a significant effect on the collision outcome because the energy of the projectile may not completely intersect the target when the impact is oblique. For example, in the collision geometry shown in Figure 2, the top of the projectile does not directly hit the target (above the dotted line). As a result, a portion of the projectile may shear off and only the kinetic energy of the interacting fraction of the projectile will be involved in disrupting the target. Thus, a higher specific impact energy is required to reach the catastrophic disruption threshold for an oblique impact.
The new PKDGRAV simulations conducted for this study were used to develop a generalized catastrophic disruption law, as previous work did not independently vary critical parameters. Table 1 presents the subset of the simulations discussed in detail below; for a complete listing see Table A.3 in the Appendix. The simulations are grouped by impact scenario: fixed mass ratio and impact angle. The value for the catastrophic disruption criteria, , is found by fitting a line to the mass of the largest remnant as a function of increasing impact energy in each group. The prime notation in the catastrophic disruption criteria indicates an impact condition that may be oblique ().
With the new data, we first consider how impact angle influences the universal law for the mass of the largest remnant. Figure 3 presents the normalized mass of the largest remnant versus normalized specific impact energy. Our previous simulations (all at ) are shown on the same universal law in Stewart & Leinhardt (2009). Note that comparable mass collisions with need to be considered carefully (offset for emphasis in Figure 3). Such collisions transition from merging to an inelastic bouncing regime (called hit-and-run, discussed in §3.1.2) before reaching the disruption regime. As a result, the mass of the largest remnant has a discontinuity between and with increasing impact energy. So, in the case of equal-mass collisions111A robust fit requires several points between and . A similar procedure was applied to fit the for the , and simulations from Marcus et al. (2010b) that are shown in Figure 4., only the fragments with are fit by a line of slope -0.5.
Our new results demonstrate that the same universal law for the mass of the largest remnant found for head-on collisions can be generalized to any impact angle:
In detail, the mass of the largest fragment for a specific subset of simulations may deviate slightly from the universal law (Figure 3). Note that the deviations vary between subsets, with some results systematically sloped more steeply and others sloped more shallowly. The deviations in from equation 5 are about 10% for near-normal impacts ( and ) and somewhat larger and more varied for highly oblique impacts.
Overall, the universal law provides an excellent representation for mass of the largest remnant for all disruptive collisions in the gravity regime. As a result, we have chosen to use the range of impact energies that satisfy the universal law for the mass of the largest remnant as the technical definition of the disruption regime. At higher specific impact energies, the linear universal law breaks down in a transition to the super-catastrophic regime (, see §4.1). At lower specific impact energies, the outcomes are merging or cratering (§5). Using our definition, the disruption regime encompasses less than a factor of two in specific impact energy. The outcomes in the disruption regime span partial accretion of the projectile onto the target to partial erosion of the target body.
Note that the derived values for are strong functions of both the mass ratio and the impact parameter (Table 1). The catastrophic disruption energy rises with smaller projectiles and larger impact parameters. Benz & Asphaug (1999) investigated the effect of impact parameter on the disruption criteria; however, their study fixed the impact velocity and varied the mass ratio of the bodies. Hence, the individual roles of the impact parameter and mass ratio cannot be discerned from their data. In the next two sections, the influence of each factor is isolated and quantified.
– mass of projectile normalized by mass of target; – impact parameter; – projectile impact velocity; – mass of largest remnant normalized by total mass; – mass of the second largest remnant; – slope of cummulative size distribution; – center of mass specific energy; – empirical critical center of mass specific energy for catastrophic disruption and gravitational dispersal derived from the simulations. In all cases, the target contained particles, kg, m; indicates models shown in blue in Figure 5; indicates and ; erosive hit-and-run regime, the disruption regime model does not apply; indicates an for which but thus ; – not enough material to fit a power law; indicates models shown in Figure 7.
3.1.2 Dependence of disruption on impact angle and derivation of the interacting mass
In order to describe the dependence of catastrophic disruption on impact angle, we introduce two geometrical collision groups (Figure 2): non-grazing – most of the projectile interacts with the target, and grazing – less than half the projectile interacts with the target. Following Asphaug (2010), the critical impact parameter,
is reached when the center of the projectile (radius ) is tangent to the surface of the target (radius ). Grazing impacts are defined to occur when .
When considering a non-grazing impact scenario with a particular and , the collision outcome transitions smoothly from merging to disruption as the impact velocity increases. For grazing impacts, however, the collision outcome transitions abruptly from merging () to hit-and-run () and then (less abruptly) to disruption (see §5). Thus, only collision energies that result in should be used in the derivation of in the grazing regime, as done for the results shown in Figure 3.
During oblique impact events, a significant fraction of the projectile may not actually interact with the target, particularly for comparable mass bodies. For gravity-dominated bodies, the projectile is decapitated and a portion of the mass misses the target entirely. As a result, only a fraction of the projectile’s total kinetic energy is deposited in the target, and the impact velocity must increase to reach the catastrophic disruption threshold.
Using a simple geometric model, we derive the fraction of the projectile mass that is estimated to be involved in the collision. First, we define as the projected length of the projectile overlapping the target. As shown in Figure 2,
where . Placing the origin at the bottom of the projectile on the center line and the positive z-axis pointing to the top of the page, the estimated projectile mass involved in the collision, , is determined by integrating cylinders of height and radius from to along the z-axis,
where is the bulk density of the projectile. The radius of each cylinder can be defined in terms of the radius of the projectile and the height from the origin,
Dividing by the total mass of the projectile, ,
Thus, is the mass fraction of the projectile estimated to be involved in the collision (see Table 1 for the values of in our simulations). The entire projectile interacts with the target when ; then and .
In order to account for the effect of impact angle on , we include the kinetic energy of only the interacting mass. The appropriate reduced mass is then
Now consider the difference between a head-on impact by a projectile of mass at and a head-on impact by a projectile of mass . At the same impact velocity, the impact energies between the two cases differ by the ratio of the reduced masses,
Next, in order to conserve the effective specific impact energy, the impact velocity must increase with increasing impact angle such that
However, the disruption criteria itself depends on the magnitude of the impact velocity (equation 3). In other words, when the effective projectile mass changes, the change in the impact velocity required for disruption varies by more than the factor presented in equation 14. Combining these two effects leads to the relationship between the oblique and head-on disruption energy for a fixed mass ratio collision,
By definition, the critical impact velocity for an oblique impact must satisfy equation 2:
The correction for changing the mass ratio is derived in the next section.
Our model for the effect of impact angle is used to derive equivalent head-on values from our new and previously published catastrophic disruption data. Using the values for and fitted to the oblique simulation results, the equivalent head-on impact disruption criteria are
We considered the catastrophic disruption of a wide variety of planetary bodies from the studies summarized in Table A.2. First, we fit the general expression for (equation 3222In the fitting procedure, the strength term is neglected in equation 3. For the lines plotted in Figure 4, the strength regime parameters are fixed at , Pa s, and based on the work in Stewart & Leinhardt (2009).) to the (equivalent) head-on disruption data to derive the values of and that best describe the entire data set, from planetesimals to planets. The same value for the material parameter is used in the angle correction and the fit to equation 3. A small number of data points were excluded from the global fit, which are discussed in §6.2.1. The best fit values for and were found by minimizing the absolute value of the log of the fractional error, .
In some cases, the impact angle correction is significant (e.g., the impact scenarios with small values of given in Table 1). With the exception of the constant-velocity results from Benz & Asphaug (1999) and Jutzi et al. (2010) and the mixed velocity data from (Benz, 2000), the disruption data were derived from simulations conducted with a constant mass ratio and the critical impact velocity for catastrophic disruption, , was found by fitting to the universal law. For the simulations described in Table 1, the model correction for impact angle usually yields an impact energy within a factor of 2 of the simulation results for head-on collisions (e.g., within the linear regime for the mass of the largest remnant). We restricted our fits to cases where to reduce any error contribution from a poor model correction for highly oblique impacts.
The compiled data and best fit model are presented in Figure 4A and B. The combined data are well fit by and (). Note the good match in the values for (colors) from the simulations with the lines of constant . Similarly good fits are found for and with . Amazingly, the compilation of catastrophic disruption data is well fit by equation 3 for single values of and for a wide variety of target compositions and over almost 5 orders of magnitude in size and 9 orders of magnitude in impact energy. The critical impact velocities span 1 m s to several 10’s km s. The best fit value for falls near pure momentum scaling ().
Upon closer examination, we found that the global fit with equation 3 systematically predicts a low disruption energy for small bodies ( km) and a high disruption energy for planet-sized bodies. Next, we consider separately the data for small and large bodies. A better fit is found for the small body data in Figure 4A with and with . The small body data includes hydrodynamic to strong bodies and different compositions. The planet data in Figure 4A are best fit with and with the very small error of . The planet size data includes three different target compositions. The data for collisions between small strong bodies have the largest dispersion; these data will be discussed in §6.2.1.
3.1.3 Dependence of disruption on mass ratio
By fitting such a large collection of data, it is clear that equation 3 describes a self-consistent family of possible values. For a specific impact scenario, the correct value for at each is ambiguous because depends on both a material property and the mass ratio. In studies that hold constant and vary the mass ratio, the derived value for the critical impact energy only applies for the corresponding critical mass ratio. As noted in previous work, the critical impact velocity falls dramatically as the mass ratio approaches 1:1 (Benz, 2000; Stewart & Leinhardt, 2009). As a result, collisions between equal-mass bodies require the smallest impact velocity to reach the catastrophic disruption threshold.
Because a mass ratio of 1:1 defines the lowest disruption energy for a fixed total mass, we derive the disruption criteria for different mass ratios with respect to the equal-mass disruption criteria, . We begin with the equality between the impact energy and gravity term in the disruption energy (equation 3),
Then, substituting for and ,
Then, for the same total mass, the relationship between the equal-mass disruption energy and any other mass ratio is determined by the difference in the critical impact velocities,
The equations for and are given in the next section.
In the compilation of catastrophic disruption data shown in Figure 4C and D, all the data have been converted to an equivalent equal-mass impact disruption energy and the colors denote . For example, the critical disruption energy from head-on PKDGRAV simulations with are a factor of three above the disruption energy for in Figure 4A (parallel sets of from Stewart & Leinhardt, 2009). The data lie on the same line after the correction in Figure 4C. The correction also brings together data from studies using different numerical methods and vastly different material properties. For example, the high-velocity for strong and weak basalt targets () fall on the same line as the PKDGRAV rubble piles after the conversion to an equivalent equal-mass impact. Similarly, studies of the disruption of Mercury (Benz et al., 2007) follow the same curve as disruption of earth-mass water/rock planets (Marcus et al., 2010b). The general form for the equal-mass disruption criteria is derived in the next section.
3.1.4 The principal disruption curve
In the previous two sections, we calculated the disruption criteria for head-on equal-mass collisions by adjusting the critical disruption energy to account for different impact angles and mass ratios. The head-on equal-mass data points, derived from the compilation of numerical simulations, fall along a single curve that we name the “principal disruption curve” (black lines in Figure 4C and D).
Then, substituting ,
Thus, along a curve with a fixed projectile-to-target mass ratio, the critical impact velocity has a linear dependence on . The linear dependence of on for a fixed mass ratio was confirmed by the numerical simulations in Stewart & Leinhardt (2009) ( in Figure 4).
Then, consider the dependence of the catastrophic disruption criteria on size (equation 3) and replace the velocity term with size,
Thus, the catastrophic disruption criteria scales as radius squared along any curve with a fixed projectile-to-target mass ratio.
Next, note the proximity of the gravity-regime equal-mass disruption energy to the specific gravitational binding energy,
shown as the grey line in Figure 4. We define a dimensionless material parameter, , that represents the offset between the gravitational binding energy and the equal-mass disruption criteria. Then, the principal disruption curve is given by
The parameter is a measure of the dissipation of energy within the target.
Finally, substituting into equation 25 gives
Hence, the critical velocity along the disruption curve for equal-mass impacts is solely a function of and .
The principal disruption curve (equation 28) is a simple, yet powerful way to compare the impact energies required to disrupt targets composed of different materials. Each material is defined by a single parameter . In Figure 4C and D, the best fit values are and for small bodies with a wide variety of material characteristics and and for the hydrodynamic planet-size bodies. These simulations span pure hydrodynamic targets (no strength), rubble piles, ice, and strong rock targets. Hence, for all the types of bodies encountered during planet formation, is limited to a small range of values. Note that the difference in between the small and large bodies is not simply because of the differentiated structure of the large bodies; two pure rock cases (, Marcus et al., 2009) fall on the same curve. Rather, the large bodies were all studied using a pure hydrodynamic model, whereas the small bodies were studied using techniques that incorporated material strength in various ways. A transition from a higher value for for small bodies to a lower value for planet-sized bodies is appropriate for planet formation studies, as discussed in §6.
Now it is clear that most of the differences in the catastrophic disruption threshold found in previous work are the result of differences in impact velocity and mass ratio (few studies varied impact parameter).
Here, we have derived a general formulation for the catastrophic disruption criteria that accounts for material properties, impact velocity, mass ratio, and impact angle. The forward calculation of for a specific impact scenario between bodies with material parameters and is described in the Appendix and in the companion paper (Stewart & Leinhardt, 2011).
3.2. Fragment size distribution
In the disruption regime, our new simulations resolve the size distribution of fragments over a decade in size (Figure 5). In general, the post-collision fragments smaller than the largest remnant form a smooth tail that can be fit well by a single power law. The second-largest remnant forms the base of this tail. For most collisions there is a significant separation in size between the largest and second largest remnant. However, if the collision is very energetic, the largest remnant joins the power-law distribution (e.g., in ). In addition, for the hit-and-run impacts with , the two largest remnants are comparable in size. Only the most energetic scenarios with fall in the disruption regime (e.g., black lines in and 0.7). In all disruption regime collisions, the slope of the cumulative power-law tail, (see Table 1), is effectively independent of the impact conditions (, , ).
Using the method from Wyatt & Dent (2002) and Paardekooper, Leinhardt, and Thebault (in preparation), the mass of the second largest remnant, , is fully constrained by knowledge of the mass/size of the largest remnant and the power-law slope for the size distribution of the smaller fragments. Let us consider a differential size distribution
where is the number of objects with a radius between and , is the slope of the differential size distribution, and is the proportionality constant. Integrating equation 31, the number of bodies between and is
Therefore, the number of objects larger than the second largest remnant () is . Assuming that ,
For spherical bodies with bulk density , the mass of material between and is
In order to enforce a negative slope of the remnants, must be less than 3. Mass is conserved in the impact; thus, the mass in the remnant tail must equal the total mass minus the mass in the largest remnant(s), , where is the number of objects with mass equal to the largest remnant (here, we allow for multiple largest remnants). Substituting for from equation 33, is given by
Substituting this expression for into equation 33 and assuming that all of the objects are spherical, the size and mass of the second largest remnant is expressed in terms of the total mass by
where . In the simulations presented here, the calculated diameter of the fragments is not informative because most PKDGRAV particles merge with other particles in gravitationally bound clumps; in these cases, the bulk density is assumed for the size of the merged particle. The mass of the remnants is accurate, however. Rewriting equation 36 in terms of mass,
where is given by the universal law and the catastrophic disruption criteria (equation 5).
In the last column of Table 1, the predicted mass of the second largest remnant (equation 37) is compared to the numerical simulations using the empirically fit , , , and . Since the analytic method presented here assumes an infinite size distribution in the fragment tail, we selected the value of to optimise the fit to the value of in the simulations. This simple method of predicting works well for all impact conditions.333Because the analytic model for the fragment size distribution assumes an infinite range of sizes in the tail, is constrained to be less than 3, which is slightly smaller than the slope of power laws that are best fit to the data (Table 1). The fragment size distribution may be modeled under different assumptions: e.g., choosing a minimum diameter in the integral of equation 31, which would represent the smallest constituent particles or grain size. We have chosen not the impose any assumptions about material properties in the model presented here, but there may be situations were more is known about the colliding bodies and the model for determining and may be modified. To illustrate the model, the predicted size distribution (blue line and triangles) is compared to selected numerical simulation results in Figure 5. In order to predict the size distribution of fragments in the hit-and-run regime impacts between comparable mass bodies ( and ), the model needs to be modified slightly. In this special case, we suggest adopting and because the target and projectile each have a nearly identical size distribution of fragments. In this example, the model is calculated for the same impact conditions as the cyan data set with , , and m s (see section 4.2 for more detailed discussion of the hit-and-run regime).
The fragment size distributions calculated using our subsonic -body simulations are consistent with shock code calculations investigating asteroid family formation via catastrophic impact events (Nesvorný et al., 2006; Jutzi et al., 2010). All of the asteroid family-forming simulations used a hybridized numerical technique, combining an SPH code with PKDGRAV in order to model the propagation of the initial shock wave and the subsequent gravitational reaccumulation of the collision remnants. The asteroid family-forming collisions have significantly different impact parameters compared to our simulations: was orders of magnitude larger, was an order of magnitude smaller than our smallest , and targets were larger (10’s km in diameter). These differences notwithstanding, we find that the range in the values of for the tail of the size distribution is very similar to our -body results (note that some published values for include the largest remnant in the fit, whereas we do not). Qualitatively, we also find a general trend in curvature of the size distribution consistent with Durda et al. (2007), with slightly convex size distributions for super-catastrophic impact events (section 4.1).
3.3. Fragment velocity distribution
Next, we consider the velocity of the collision remnants. The results are easier to interpret by separating the velocities of the largest remnant from the rest of the collision remnants.
We first consider the speed of the largest remnant with respect to the center of mass of the collision (Figure 6). For erosive events (), there is almost no change in the amplitude of the target velocity for impacts with . Even at , the velocity reduction is minimal for all fractions of mass lost. Because is the center of the probability distribution of impact angles, fully half of all erosive impacts have change in the target velocity amplitude. After head-on collisions (), the largest remnant moves with the center of mass velocity. Note there is significant scatter in the data from , which is due to the fact that there was a small number of particles delivering the impact energy to a localized region of the target; thus, the organization of the surface features on both objects become important. For disruptive impacts at , there is partial velocity reduction of the largest remnant. From these data, we cannot define a unique function for the dependence of on , and we suggest that a quasi-linear relationship for is a reasonable approximation. We stress that the specific dependence of on in the disruptive regime is likely to be sensitive to internal structure and composition, so extrapolation of these results beyond weak, constant density objects should be done with caution.
In complete merging events, of course, the post-impact velocity is zero with respect to the center of mass. The data with steadily approach the center of mass velocity with more mass accreted. The and 0.9 data points plotted near are primarily hit-and-run events, which will be discussed in §4.2.
The smaller remnants of disruptive collisions have a more complex behavior. Figure 7 presents mass histograms of fragments versus velocity with respect to the largest remnant from the simulations summarized in Table 1. The slowest simulations are not plotted for the and 1.0 grazing impacts because there are only a small number of fragments. A significant number of the fragments consist of 10 PKDGRAV particles or less; in Figure 7, mass associated 10 or less particles is shown as dotted histograms. The dotted histograms overlay the total mass histograms for all but the lowest velocity bins; thus, within most of the velocity bins, the simulations do not have the resolution to robustly constrain the size-frequncy distribution of the mass in the bin. The smallest (poorly resolved) fragments are found in all velocity bins, while the largest fragments tend to move slowly with respect to the largest remnant. For example, the second largest remnant falls in one of the lowest velocity bins, but that bin is also occupied by smaller fragments.
Hence, to describe the velocity field after a collision, we fit the velocity-binned mass of the collision remnants. The binned mass versus velocity is a fairly well defined exponential function for most of the simulations. In general, the lowest velocity bin in Figure 7 is of order . Using a least-squares fit of the subset of simulations in Table 1, we find the mass fraction in the lowest velocity bin is proportional to the largest remnant mass:
To determine the slope, , of the binned mass versus velocity exponential function, we integrate the differential mass function,
where , , is the bin width, and the total mass in the histogram is the total mass in the remaining remnants, .
The fragment velocity scaling law (equation 39) is shown in blue in Figure 7 for selected cases indicated by in Table 1. The velocity distributions of the remnants agree qualitatively with those found in hypervelocity simulations of asteroid family forming events, although previous workers have not fit any function to the velocity distribution of the fragments (e.g., Michel et al., 2002; Nesvorný et al., 2006).
4. Other Collision Regimes
4.1. Super-catastrophic regime
In both laboratory experiments in the strength regime (e.g., Kato et al., 1995; Matsui et al., 1982) and the few high resolution disruption simulations in the gravity regime (e.g., Korycansky & Asphaug, 2009), the relationship between the mass of the largest remnant and the specific impact energy shows a marked change in slope at around . We define the super-catastrophic regime when (e.g., when by the universal law, equation 5). In the super-catastrophic regime, the mass of the largest remnant follows a power law with rather than the linear universal law.
The slope of the power law for the largest remnant mass vs. impact energy shows some scatter in laboratory data, primarily in the range of -1.2 to -1.5. In Figure 8, our few simulations of super-catastrophic collisions (symbols) are compared to the range of outcomes from laboratory experiments (dotted and dashed lines). Based on the simulations in the gravity regime and laboratory experiments in the strength regime, we recommend a power law in the super-catastrophic regime,
where and the coefficient is chosen for continuity with the universal law (equation 5). The slope of the power law, about -1.5, is consistent with our gravity regime simulations and a wide range of laboratory experiments summarized in Figure 1 in Holsapple et al. (2002).
In Figure 8, the solid line is the combined universal law and the recommended super-catastrophic power law (equations 5 and 44). The dotted line is our fit to disruption data on solid ice from Kato et al. (1992), . The dashed line is our fit to disruption data on basalt from Fujiwara et al. (1977), . Note that lab data are available up to very high values of . The lab data spanning very weak to very strong geologic materials can be considered lower and upper bounds for the parameters in equation 44.
The general agreement between the gravity and strength regimes suggests that gravitational reaccumulation of fragments has a negligible effect in the super-catastrophic regime. In other words, the mass of the largest fragment is primarily controlled by the shattering process.
Based on the similarity of the size distribution of fragments in laboratory experiments to the gravity regime data presented here (Figure 5), we suggest that the dynamical properties of the smaller fragments in super-catastrophic collisions are similar to the disruption regime. Therefore, the size and velocity distributions described in §3.2 and §3.3 can be applied.
4.2. Hit-and-run regime
Non-grazing impacts in the gravity regime transition from perfect merging to the disruption regime with increasing impact velocity. However, for impact angles greater than a critical value, an intermediate outcome may occur: hit-and-run (Agnor & Asphaug, 2004a; Asphaug et al., 2006; Marcus et al., 2009, 2010b; Asphaug, 2010; Leinhardt et al., 2010). In a hit-and-run collision, the projectile hits the target at an oblique angle but separates again, leaving the target almost intact. Some material from the topmost layers of the two bodies may be transferred or dispersed. Depending on the exact impact conditions, the projectile may escape largely intact or may sustain significant damage and deformation (e.g., Figure 7 in Asphaug, 2010).
The hit-and-run regime is clearly identified by considering the accretion efficiency of a collision, defined by Asphaug (2009) as
In a perfect hit-and-run event (), . For a perfect accretion event (), . An erosive event in which leads to . Note that the negative value of that corresponds to catastrophic disruption () depends on the specific mass ratio of the two bodies ().
There is remarkably good agreement in the accretion efficiency and transitions from merging to hit-and-run and from hit-and-run to disruption between this work and previous simulations of higher velocity impacts between large planetary bodies (Agnor & Asphaug, 2004a, b; Marcus et al., 2009, 2010b). Figure 9 shows the accretion efficiency from our simulations in solid colored lines for four different projectile-to-target mass ratios and impact parameters. Data for collisions between protoplanets at supersonic velocities from Agnor & Asphaug (2004a, b) (and plotted in Asphaug, 2009) are shown in dashed lines for the common mass ratios (1:1 and 1:10). Hit-and-run collisions are indicated by a sudden drop from merging outcomes () to a nearly constant value of for a range of impact velocities. Note that the drop in is sharpest for equal-mass bodies. For smaller mass ratios, the transition is not as sharp, and partial accretion of the projectile occurs at energies just above perfect merging ( just above 0).
Outcomes that are defined by the disruption regime have steep negative sloped accretion efficiencies. The disruption regime equations apply for partial accretion () and for erosion of the target (). Note that for high impact parameters (e.g., ), there exists an intermediate regime where the accretion efficiency has a very shallow negative slope and values of just below 0. These impact events, termed erosive hit-and-run, lead to some erosion of the target and more severe deformation of the projectile. The erosive hit-and-run regime is eventually followed by a disruptive style erosive regime at sufficiently high impact velocities. The post-hit-and-run disruptive regime may be identified by finding the impact energy that leads to a linear relationship that satisfies the universal law. The required impact velocity increases substantially with increasing impact parameter; see §5 and Marcus et al. (2010b) for an example disruption regime after a hit-and run regime ( and ).
In an ideal hit-and-run event, the target is almost unaffected by the collision, and the velocity of the largest remnant (the target) is about equal to the initial speed of the target with respect to the center of mass. More commonly, there is a small velocity change in both bodies which increases the probability of merging in subsequent encounters (Kokubo & Genda, 2010). Agnor & Asphaug (2004a) referred to this collision outcome as inelastic bouncing. In our hit-and-run simulations with (green cluster of points in Figure 6 at ), the targets typically lose about 10% of their pre-impact velocity. For (red points), there is more significant slowing of the target. Our data does not provide a robust description of the dependence of the post-impact velocity on the impact parameter and impact velocity.
The projectile may be significantly deformed and disrupted during a hit-and-run event. The level of disruption of the projectile may be approximated by considering the reverse impact scenario: a fraction of the larger body impacts the smaller body. In this case, we estimate the interacting mass from the larger body with a simple geometric approximation. For the example geometry given in Figure 2, the cross-sectional area of the circular projectile interacting with the target is calculated. The apothem is given by , and the central angle is . Then, the projectile collision cross section is
The interacting length through the target is approximated by the chord at ,
And the interacting mass from the target is of order
Note that the interacting mass depends on the impact angle (through ).
To estimate the disruption of the projectile, we consider an idealized hit-and-run scenario between gravity-dominated bodies: the fraction of the target that does not intersect the projectile is sheared off with negligible change in momentum and gravitationally escapes the interacting mass. Hence, we ignore the escaping target mass and consider only the impact between and the projectile mass, .
The reverse impact is thus defined by