A New Database of Giant Impacts over a Wide Range of Masses and with Material Strength: A First Analysis of Outcomes

The simulation of the gravitational interactions of a planetary system (N -body simulations) depends on how collisions are treated. If the colliding bodies are assumed to be perfectly merging, the mass is conserved, and the new orbit is often placed at the center of mass of the colliding bodies (e.g., Duncan et al. 1998). Perfect merging may be a sufficient assumption for understanding the largest-scale architectures in the solar system (e.g., Kokubo & Genda 2010; Quintana et al. 2016; Walsh & Levison 2019). However, the approximation alters the dynamical evolution and the formation sequence.

Including the realistic collision outcomes increases the formation timescales (Chambers 2013; Quintana et al. 2016), because inefficient accretions—especially hit-and-runs—are common. If the runner returns to the target, or to another nearby accreting body (Emsenhuber et al. 2021), it is a "collision chain," a hit-and-run followed by another similar-sized collision (merger, hit-and-run, disruption, etc.). During planet formation at ∼1au around a solar-mass star, the recollision time for a collision chain can be 1 × 10³–1 × 10⁶ yr (Emsenhuber & Asphaug 2019a). Such multicollisional pathways could lead to mantle-stripped cores, the "stranded runner" hypothesis for the origin of Mercury (Asphaug & Reufer 2014; Chau et al. 2018), and for metallic asteroids and meteorites (Yang et al. 2007).

Also, a realistic treatment of collisions affects the resulting composition (Carter et al. 2015; Dwyer et al. 2015), and the mixing of material between planetary accretion zones (Burger et al. 2020), and leaves behind a diversity of smaller bodies (e.g., Mars-sized; Emsenhuber et al. 2020). This is because a growing fraction of the remainder end up surviving hit-and-run collisions as the growth of the largest bodies proceeds (Asphaug & Reufer 2014; Asphaug 2017). This could explain the observed increasing diversity of terrestrial planets with decreasing mass (Gabriel & Cambioni 2023).

Implementing realistic treatments of similar-sized collisions can be achieved in N-body codes in several ways. For example, scaling laws have been developed (e.g., Leinhardt & Stewart 2012; Reinhardt et al. 2022) from 3D giant impact simulations that characterize the outcomes of giant impacts. One application of our new database would be to validate or update the parameters to such scaling laws or develop new ones, as our parameter space includes friction in all simulations, providing more accurate results for impacts at smaller scales. Another approach is to apply machine learning to the outcomes of giant impact simulations (Timpe et al. 2020). Surrogate models (e.g., Cambioni et al. 2019, 2021) can be generated, which relate preimpact conditions to postimpact conditions such as the largest remnant mass m_L, the second remnant mass m_S (runner), dynamical properties, and compositional information. The present database, including an initial assessment of debris, is designed specifically for building machine-learning models and spans a larger range of collisions than an existing giant impact database used for machine-learning models (Cambioni et al. 2019; Emsenhuber et al. 2020). One important limitation, however, is that the impacting bodies are not rotating before impact in our database. Introducing rotation would expand the parameter space considerably, but is a factor to consider given its effect on postimpact outcomes (Timpe et al. 2020), and the profound effects it may have on the internal structure of the resulting bodies (Lock & Stewart 2017).

Certain aspects of giant impacts, such as the mass of remnants, can be modeled by analytical relationships (such as scaling laws). However, obtaining a unified model that can predict the outcomes of impacts across various regimes, such as in small asteroid-scale collisions and in shock-inducing Moon-formation events, is challenging. This is due to numerous physical complexities inherent to collisions at different scales and in different bodies (Gabriel et al. 2020). Transitions in the equations of state (EOSs; Stewart et al. 2020), vapor production (Benz et al. 2007; Carter et al. 2020; Davies et al. 2020), and the presence of shocks in large-scale collisions produce ample vapor and may alter the nature of erosion (Gabriel et al. 2020; Gabriel & Allen-Sutter 2021). At small scales, friction and strength make erosion less likely for a given scaled velocity (e.g., Jutzi 2015), and other forms of dissipation become substantial (Melosh & Ivanov 2018). Even in Mars-scale collisions where rheology has not classically been noted to alter the mass of the largest remnant (e.g., Benz & Asphaug 1999), other aspects of the collision such as heating and debris generation are influenced by the presence of strength (Emsenhuber et al. 2018). To make meaningful progress toward a unifying model across these complex regimes, a set of simulations that spans over a large range of preimpact conditions is required.

Still, giant impact outcomes do have systematic trends across vast ranges. For instance, Jutzi (2019) identified three basic regimes of giant impacts, or similar-sized collisions, for the velocity ranges that can potentially result in accretion. For bodies smaller than ∼100 km, there is the porosity regime, where the outcome is mainly affected by pore crushing (e.g., Housen & Holsapple 1999; Belton et al. 2007; Jutzi et al. 2015). Intermediate-mass bodies (up to ∼1000 km) are in the strength regime, where friction is important (for the link between strength and friction, see Section 2.4). The largest bodies are in the gravity regime, where strength and porosity can be ignored. For large enough bodies where mutual escape velocities (and thus impact velocities) exceed the sound speed, the shocks dominate and result in a fourth regime (Gabriel et al. 2020). These regimes are idealizations, and the transitions are not abrupt. Indeed, the transition from strength to gravity dominance may span the entire range of "oligarchic growth," the canonical late stage of Moon- to Mars-sized bodies (Kokubo & Ida 2002). For this reason, we include a strength model in all of our simulations, even in super-Earth collisions where the effect on the outcome is expected to be negligible. This allows the data to properly capture the smooth transition from strength to gravity dominance.

The end state of the last giant impact represents the beginning of a planet's long-term geophysical evolution (e.g., Zahnle et al. 2007). Accreted planets can end up with mantles and cores unstable to convection, especially for larger-scale collisions in which the gravitational potential released by the merging material exceeds the energy of melting (Lock et al. 2020). Smooth particle hydrodynamics (SPH), the technique most widely used to model giant impacts, is unable to properly represent long-term convective instability and other effects that happen on a much longer timescale than the collision. To study postimpact geodynamics requires a hybrid numerical approach, as has been applied by Golabek et al. (2018) to model consequences like geothermal overturn, geochemistry, and solidification after giant impacts. For solid final bodies, finite element analysis could be used to predict postcollisional relaxation from these SPH results.

Giant impacts are transformative events, and simulations of giant impacts have become the basis for making quantitative geochemical predictions about late-stage planet formation (e.g., Haghighipour et al. 2018), especially the origin of the Moon (Canup et al. 2021). But colliding planetary bodies in simulations are represented by idealized compositions, and these and other simplifying assumptions need to be recognized. We represent terrestrial planets and their progenitors as differentiated nonporous spheres of forsterite and iron composition, with varying core mass fractions. This simple interior structure allows us to focus on attaining accuracy and reliability of the model predictions, while making a sufficient sweep of the parameter space to obtain results that may be applied generally to models of terrestrial planet formation. Furthermore, there remains a limited availability of shock EOS models for a wider range of mantle materials, and the forsterite EOS (Stewart et al. 2020) is the most up-to-date and widely used for our purposes.

A two-component rock/iron interior is a justifiable approximation for terrestrial bodies in the modeled size range. For instance, 525 km diameter (4) Vesta has a rocky mantle and an iron core (e.g., Russell et al. 2012); so do Mercury, Venus, the Moon, and Earth, to good approximation. So we maintain this assumption from the smallest colliding bodies in our database, about 1000 km diameter, up to ∼5 M_⊕, beyond which point gas-free accretion is unlikely (e.g., Rogers 2015).

Numerical studies have revealed the complex thermodynamic evolution that occurs in giant impacts (Carter et al. 2020) and how outcomes can depend on the choice of initial thermal conditions, even for relatively simple two-component planets. For example, a massive magma ocean on the target can enhance the postimpact disk mass and its Earth-isotopic fraction in simulations of Moon formation (Hosono et al. 2019). That said, the material from which the Moon forms represents only ∼1% of the material of the colliding bodies in the Canonical model for Moon formation (Canup 2005). For our database, we thus implement only one thermal state across all colliding bodies: both the core and the mantle start slightly below the solidus, on the basis that convective cooling may be faster than the time between giant impacts.

We present a significant new database of terrestrial planet-forming giant impacts. To take into account the limitations and bottlenecks of previous works, our new database has the following characteristics:

• 1.  
  It is applicable to a wide range of impactor-target properties and impact parameters, from the sizes of the largest asteroids and rocky satellites, to terrestrial planets of several Earth-masses (M_⊕), which are not expected to be large enough to accrete significant atmospheres (e.g., Rogers 2015).
• 2.  
  It includes the transitions between collision regimes, especially between graze-and-merge and hit-and-run, i.e., accretion and nonaccretion.
• 3.  
  It provides the outcomes of collisions for bodies with various core mass fractions, thus enabling self-consistent treatments of core size evolution in collision chains.

Within this scope, we construct a database of 1250 simulations (split in two sets, one comprising 1000 simulations and one 250) suitable for the development of machine-learning models. To compute the final reduced properties of each simulation, we proceed in multiple steps. First, we define the setup of the study (Section 2.1) and determine the initial conditions of hydrodynamical simulations (Section 2.2). Then, we perform simulations using SPH (Section 2.3) and analyze their results (Section 2.3). We also present several results from our calculation: general results (Section 3.1); specific items relevant to low-velocity grazing collisions, satellite capture (Section 3.2), body shapes (Section 3.3), and debris (Section 3.4).

The first step to generating our database of giant impact simulations is to select the parameters that will be explored. Target mass m_tar and the impactor-to-target mass ratio γ = m_imp/m_tar are the primary parameters. For the composition, we use a single parameter, the core mass fraction Z_tar and Z_imp, which describes the fraction of iron in our two-layered bodies, the rest being forsterite. Core mass fraction is a free parameter in our database as it will allow us to understand the collision outcomes between planets that have been previously eroded, e.g., to know the outcomes of collisions between core-rich bodies. The impact velocity is given in terms of the mutual escape velocity v_coll/v_esc, which thus increases with total mass. The impact angle θ_coll is defined as the angle between the line connecting the centers of mass and the relative velocity vector at impact. We use the convention that θ_coll = 0° is a head-on collision, and θ_coll = 90° is a grazing collision.

For the present work, we have decided not to include the effect of preimpact rotation. Timpe et al. (2020) found that preimpact rotation had less influence on the giant impact outcomes (apart from final spin, which is strongly correlated) than the colliding mass, mass ratio, core mass fraction, impact angle, and impact velocity. An important advantage of leaving out preimpact rotation is the smaller dimensionality of the parameter space. Accounting for all possibilities of preimpact rotation when the target and impactor are similar in size would require six additional parameters, three for each body, to describe the spin angular momentum vectors.

We aim to study the effect of the target's mass, so we explore a large range from 1 ×10⁻⁴ M_⊕ to 5 M_⊕, which is sampled uniformly in logarithm space. The lower boundary of target mass corresponds to about the mass of (1) Ceres. This range was selected because it incorporates the size range where the effect of friction (included in every simulation in the generated database) is important (e.g., Emsenhuber et al. 2018), yet where porosity is likely unimportant for these massive, differentiated bodies. The upper boundary was selected so that the collisions can be applied to the formation of extrasolar planetary systems that contain super-Earths. We do not include simulations of planets with gas envelopes, so we limit our study to masses below 5 M_⊕ because more massive planets tend to have significant envelopes (e.g., Rogers 2015).

As for the mass ratio γ, our sampling is uniform so that unequal mass ratios are studied as much as more similar mass ratios, as suggested by Valencia et al. (2019). We selected the range of 0.05 < γ < 1, covering the ∼1:3 diameter range defining similar-sized collisions (Asphaug 2010) and as constrained by our capability to computationally resolve the smaller body. The least massive impactor in our database is 5×10⁻⁶ M_⊕, somewhat more massive than Main Belt asteroid (16) Psyche.

The core mass fractions of Z are sampled from a piecewise uniform distribution over the range of 10%–90%, where the range is selected so that the core and mantle are numerically resolved at least several particles across for the chosen resolution. Further, to account for the population of bodies around the average value of Fe/Mg for stars with planets from the Hypatia Catalogue (Hinkel et al. 2014; see also Figure 2 in Scora et al. 2020), and the average chondritic values of the terrestrial planets, we center the distribution so that half the bodies have a core mass fraction below 30%, and the rest are above. For this median value, which is around the value for Earth and Venus, the core radius is about half the body radius.

According to N-body studies, the majority of giant impacts during terrestrial planet formation occur in the range 1–2v_esc (e.g., Chambers 2013; Emsenhuber et al. 2020), although the median value depends on the damping effects of planetesimals (O'Brien et al. 2006). According to Asphaug (2010), the diversity of planet formation by giant impacts is due to the fact that the transitions in outcomes happen around this accretionary range of velocities. Genda et al. (2012) noted a significant variability at the low-velocity merging and graze-and-merge transition, for minor changes in simulation parameters. Similarly, Cambioni et al. (2019) found the highest rate of misclassification and prediction error to occur at the transition between hit-and-run and merging, in their machine-learning analysis of a prior database of giant impact outcomes. Therefore, we place emphasis on mapping out this transition by constructing a distribution of impact velocity that favors low-velocity collisions. We note that minor deviations in the transition as a function of different thermal conditions are not covered in this work. We define the cumulative distribution from the relative velocity at infinity, as scaled by the mutual escape velocity v_∞/v_esc. This is related to the impact velocity by

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{v}_{\mathrm{coll}}}{{v}_{\mathrm{esc}}}\right)}^{2}=1+{\left(\displaystyle \frac{{v}_{\infty }}{{v}_{\mathrm{esc}}}\right)}^{2}\end{eqnarray} \tag{ 3 }$$

according to energy conservation. The result is a piecewise uniform distribution as provided in Table 1. The resulting cumulative distribution, given for v_coll/v_esc, is shown in Figure 1.

Zoom In Zoom Out Reset image size

Table 1. Cumulative Distribution of Impact Velocities in the Grid of Hydrodynamical Simulations

v∞/vesc	vcoll/vesc	Cumulative Fraction
0	1	0
1	∼1.4	0.4
2	∼2.2	2/3
3.5	∼3.6	0.8
6	∼6.1	0.9
$\sqrt{99}$	10	1

Download table as:  ASCIITypeset image

As for the impact angle θ_coll, we follow the distribution expected for collisions between point impactors onto gravitational targets, ${\rm{d}}P\propto \sin (2{\theta }_{\mathrm{coll}}){\rm{d}}{\theta }_{\mathrm{coll}}$ (Shoemaker 1962). This ensures that the transitions between hit-and-run and graze-and-merge and accretion, which occur around the nominal range of impact angles at nominal velocities (e.g., Leinhardt & Stewart 2012), are well sampled.

To generate the specific combinations of the parameters, which will be used to perform the hydrodynamical simulations, we use the Latin hypercube sampling (LHS) method (e.g., McKay et al. 1979; Timpe et al. 2020). LHS divides each parameter into n intervals of equal probability based on the distribution, with n being the number of samples. Then, one sample is selected randomly from each interval. This ensures that the entire range of possible values for each parameter is sampled. In addition to that, LHS adopts criteria that ensure that the entire parameter space is well sampled. For instance, the algorithm avoids correlations between parameter values in the selected samples, so that the effect of each parameter can be disentangled. In this work, we used the pyDOE Python package with the minmax setting. This package returns all the samples in the [0, 1] range with uniform probability. We convert these to the actual collision parameters using the probability distributions discussed above.

For this work, we generated two lists of initial conditions using LHS, one with n = 1000 entries, and one with n = 250. The intention was to have separate sets for future machine-learning applications. For most of the analysis, we combine the two in a set of 1250 simulations.

To perform the hydrodynamical simulations, we use the SPH technique (e.g., Monaghan 1992; Rosswog 2009). SPH uses a Lagrangian description with continuum material divided into particles. Quantities are retrieved by performing a kernel interpolation and their spatial derivatives by interpolation of the underlying quantity with the kernel's gradient. Time evolution is given by the Euler's equations, except that the density is computed by performing kernel interpolation at each step and corrected for free-surface effects using the method of Reinhardt & Stadel (2017), which works by increasing the density of a particle if there is a spatial imbalance of neighbors around it. An artificial viscosity term inspired by Riemann solvers (Monaghan 1997) with a time-dependent factor (Morris & Monaghan 1997) is included to resolve shocks, the exact form of which is described in Emsenhuber et al. (2021). The interpolation is performed using a cubic spline kernel (Monaghan & Lattanzio 1985) with each particle having about 50 neighbors. We note that SPH tends to spuriously damp the subsonic turbulence (e.g., Cullen & Dehnen 2010; Bauer & Springel 2012; Deng et al. 2019); however, this should affect more the internal mixing within the final bodies rather than their global iron–silicate fraction. Finally, material strength is modeled according to the formulation discussed in the next section.

To close the Euler equations, we need to provide an EOS that provides the pressure as a function of the density and the specific internal energy p(ρ, u). The choice of EOS is limited by the requirement that they be thermodynamically reliable in all the energy regimes applicable to a giant impact (Melosh 2007). In this work, we use ANEOS for the iron core (Thompson & Lauson 1972) and a modified version of ANEOS for Mg₂SiO₄ (forsterite; Stewart et al. 2020) for the mantle. To improve performance, we precompute tables directly from ANEOS for the anticipated range of values as in previous works (e.g., Benz et al. 1989; Reufer et al. 2012). The code nevertheless retains the capability to call ANEOS for the few cases where the values lie outside of the tabulated range.

The computer code used in this study is SPHLATCH (e.g., Reufer et al. 2012; Asphaug & Reufer 2013, 2014) and includes additional updates and corrections presented in Ballantyne et al. (2023). An earlier version of the same code was used to generate a previous collision database(Reufer 2011) on which the first machine-learning derived surrogate models were based (Cambioni et al. 2019; Emsenhuber et al. 2020), and an updated semiempirical scaling law was developed (Gabriel et al. 2020). Those simulations did not include friction, and the reported database spans 1×10⁻²–1 M_⊕ planets (for more information, see Gabriel et al. 2020).

The constitutive model we adopt is similar to that from Emsenhuber et al. (2018): elastic perfectly plastic material (Benz & Asphaug 1994, 1995), with the deformation tracked by the deviatoric stress tensor, that is reduced at the Hugoniot elastic limit with $\sqrt{{J}_{2}}/Y$, where J₂ is the second invariant of the deviatoric stress tensor (not to be confused with the global gravity coefficient), and Y is the pressure-dependent yield strength (Collins et al. 2004; Jutzi 2015). We further include a correction tensor to compute the local velocity gradient to achieve angular momentum conservation (Bonet & Lok 1999; Speith 2006).

As a simplification, we assume that cohesion is zero, as its effect is only noticeable in similar-sized collisions in the diameter range ≲1 km (Jutzi 2015). The model thus assumes fully damaged material in the solid state, governed by a friction model, given by

$$\begin{eqnarray}&&{Y}_{{\rm{d}}}={\mu }_{{\rm{p}}}p,\end{eqnarray} \tag{ 4 }$$

where μ_p is the coefficient of friction, a material parameter, and p is pressure. The subscript "d" refers to damaged material. The strength does not increase arbitrarily with pressure; it is limited by the yield strength of intact (nondamaged) material,

$$\begin{eqnarray}&&{Y}_{{\rm{i}}}=\displaystyle \frac{{\mu }_{{\rm{i}}}p}{1+{\mu }_{{\rm{i}}}p/{Y}_{{\rm{m}}}},\end{eqnarray} \tag{ 5 }$$

where μ_i is the coefficient of internal friction, and Y_m is the von Mises plastic limit. The subscript "i" refers to intact material. The full form of the yield strength is ${Y}_{{\rm{p}}}=\min \left({Y}_{{\rm{d}}},{Y}_{{\rm{i}}}\right)$, which represents a material supported by friction subject to plastic yielding.

We adopt the constitutive model for "rock materials" in Table A.1 of Collins et al. (2004); here, μ_p = 0.8, and μ_i = 2.0 for all simulations. We assume the same coefficients of friction for the iron cores as well; although it is simplistic, this choice is unlikely to matter because the pressure at the core–mantle boundary in all our bodies exceeds Y_m for iron, 0.68 GPa. The forsterite mantle is subject to friction, when at lower pressure, and Y_m = 3.5 GPa. We study the effect of our choice for μ_p and μ_i in Appendix B.

Yield strength is also temperature-dependent. To capture this effect, we further modulate the yield strength as it approaches the melting point with

$$\begin{eqnarray}&&{Y}_{{\rm{T}}}={Y}_{{\rm{P}}}\tanh \left[{\zeta }_{{\rm{T}}}\left(\displaystyle \frac{{T}_{{\rm{M}}}}{T}-1\right)\right],\end{eqnarray} \tag{ 6 }$$

where ζ_T = 1.2 is the thermal softening parameter (Ohnaka 1995; Collins et al. 2004), and T_M is the melting temperature. The melting temperature is consistently recovered from the EOS by determining the melting temperature for both materials. This is possible because the ANEOS forsterite EOS (Stewart et al. 2020) provides the full phase information, while the older EOS for quartz (Melosh 2007) applied in previous work (e.g., Reufer 2011; Emsenhuber et al. 2018) does not.

In summary, all simulations are performed with the strength of a fully damaged material modulated by yielding beyond each material's plastic limit, and further by temperature. The strength thus applies to a thin or even unresolved layer in the larger bodies of the database, and to the entire mantle and perhaps some of the core in the smallest bodies in the database. While including strength adds computational cost for planet-scale collisions where strength is not expected to influence outcomes, it ensures constant treatment across the database. Moreover, the simulations involving the most massive bodies are the least computationally expensive, due to the Courant–Friedrichs–Lewy time step criterion being proportional to the spatial resolution (which is itself proportional to the body sizes with a constant number of particles, as in this work). Thus, including friction at the largest scales does not significantly slow down the computation of the database, since most computational effort is on the smaller bodies.

Colliding bodies need to have their initial thermal state specified. This is particularly important for smaller-scale collisions as our SPH model includes friction. The yield strength in the friction model is modulated by the ratio between the temperature of the material undergoing stresses and that of the melting point (see Equation (6)).

We use isentropic profiles for each of the core and mantle. We leave the initial specific entropy of the core and mantle, S_core and S_mant respectively, to be determined. We prescribe thermal profiles such that a large portion of the cores and mantles of the preimpact bodies are in a solid phase, but not far from their melting point. Due to the wide range of body masses that represent our starting conditions, this thermal state is not achievable with a single global value of entropy for the cores or mantles. Hence, we selected the relationship presented in Figure 2. These have two different regimes, with a constant entropy for all bodies whose masses are below 1×10⁻³ M_⊕ and a log-linear relationship above. The selected entropies result in core and mantle temperatures of 1651 and 1657 K respectively for the bodies at the low-mass end of the range. For larger bodies, a temperature drop occurs at the core–mantle boundary, which is expected. Thermal conduction is not included in our SPH model, so this temperature jump between materials persists through calculations. To highlight the jump in temperature, we show a profile on an Earth-like body (mass M = 1M_⊕, and core mass fraction Z = 30%) in Figure 3. Note that, for the chosen EOS, the outer radius of the body is 6577 km, roughly 3% larger than the Earth.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The body preparation scheme follows that of Emsenhuber et al. (2021): we begin by obtaining a 1D radial profile in hydrostatic equilibrium using the scheme presented in Benz (1991). The profiles are then mapped onto 3D profiles using the methodology described in Reinhardt & Stadel (2017) and further relaxed for 6 hr to damp particle motions. To enforce the prescribed initial value, we use a fixed-entropy SPH formulation (rather than evolving energy) during this relaxation step of problem setup.

We choose the resolution so that there are roughly 100,000 total particles in the simulation (in the target and the impactor). For the smallest mass ratio scenarios (γ = 0.05), this equates to roughly 5000 particles in the impactor, since the particles are roughly (but not exactly) equal mass in our simulations. Prior studies (e.g., Asphaug 2010) have shown that accretion efficiency is resolved to a few percent for SPH simulations with 20,000 and more particles, and Meier et al. (2021) found that ${Q}_{\mathrm{RD}}^{* }$ (the catastrophic disruption threshold) is converged for 1 ×10⁵ particles in SPH simulations using ANEOS, providing additional confidence that our results for lower-energy, hit-and-run scenarios are converged in terms of the largest remnant mass. We also performed simulations with a resolution increased by a factor of 10 (roughly 1 million total particles) that show only small differences compared to the nominal resolution (Appendix C).

We set up each collision by placing the bodies at a distance equal to 5 times the sum of the radii, well outside the region where tidal influence begins to deform the bodies (more than twice the Roche limit). The collisions are then evolved for 50 τ_coll past the initial contact, which is defined as

$$\begin{eqnarray}&&{\tau }_{\mathrm{coll}}=\displaystyle \frac{2\left({r}_{\mathrm{tar}}+{r}_{\mathrm{imp}}\right)}{{v}_{\mathrm{coll}}}.\end{eqnarray} \tag{ 7 }$$

50τ_coll corresponds to a bit more than 1 day for the low-velocity collisions (that is, at the mutual escape velocity). By inspection, we find that this duration is sufficient in most cases to determine the collision outcome. However, grazing collisions, where the impactor is captured on a bound orbit with a long duration loopback orbit, demand longer integrations. These collisions encompass both graze-and-merge, where the impactor is left on an orbit that intersects with the target, and graze-and-capture regimes. Graze-and-capture refers to collisions where the impactor remains as a bound satellite, akin to the scenario of Canup (2005) for the origin of Pluto–Charon.

For cases that end up near the transition between the hit-and-run and graze-and-merge regimes, the simulations and end-states must be analyzed separately (Emsenhuber & Asphaug 2019a). The runner's loopback orbit requires days of physical time, typically, so for returning runners with more than about ∼10% of the impactor's mass, we continue the SPH evolution until it has made at least one more passage of the pericenter. Afterwards, it may be tidally disrupted or partially accreted. We discriminate two cases in these scenarios. If the orbital period of the runner around the target is smaller than about 7 days, the simulation is further evolved until, usually, no such body remains. We checked that heating due to spurious activation of the artificial viscosity (for instance, due to residual motion inside the bodies) is minimal during the loopback orbit. However, if the orbital period is determined to be longer than 7 days, the simulation is kept in this intermediate state and marked as such. This is to avoid the buildup of numerical imprecision in the gravity solver during such an extended period, which could make the final state less realistic.

We proceed as in Emsenhuber & Asphaug (2019b), Emsenhuber et al. (2020) for the analysis of the SPH simulations, and perform additional analysis of body shapes, interior structures, and debris. We begin with a friends-of-friends search of SPH particles to find contiguous bodies. Each of these bodies is replaced by one superparticle with equivalent mass and momentum for the remaining steps. The second step is to find gravitationally bound material and identify them as unique clumps. In our work, we consider a minimum clump mass worthy of analysis to be 10 times the mass of the most massive SPH particle, although in practice our results do not depend on this minimum. The particles not part of any clump are considered unresolved debris. Identified clumps can include coorbiting pairs of superparticles and multiple-body systems. A last step is to find physical bodies. When a single contiguous body is found in a clump, all remaining particles in such a clump are attributed to that body. When there are more, the particles not part of a contiguous body are attributed to the body with which they have the lowest relative energy.

These bodies and clumps provide the basic metrics of collision outcomes, which are archived in machine-readable format. The most massive clump (from the second analysis step) is taken as the largest remnant, whose accretion efficiency ξ_L is given by Equation (1). In a perfect merger, ξ_L equals 1, and a value of 0 represents a target with a postimpact mass that is equal to its preimpact mass. The second most massive clump is defined as the second remnant, for which we define the accretion efficiency of the second remnant of mass m_S:

$$\begin{eqnarray}&&{\xi }_{{\rm{S}}}=\displaystyle \frac{{m}_{{\rm{S}}}-{m}_{\mathrm{imp}}}{{m}_{\mathrm{imp}}}=\displaystyle \frac{{m}_{{\rm{S}}}}{{m}_{\mathrm{imp}}}-1.\end{eqnarray} \tag{ 8 }$$

In a hit-and-run collision, this second remnant is usually made largely of an impactor continuing downrange (i.e., a runner). Under most circumstances, mass is eroded from the body, such that ξ_S < 0. Material not part of the two most massive remnants is considered debris; either resolved if it is part of a clump or unresolved if not. By analogy to the two largest remnants, debris production is characterized by its accretion efficiency ξ_D, which is its mass in units of impactor masses, or

$$\begin{eqnarray}&&{\xi }_{{\rm{D}}}=\displaystyle \frac{{m}_{{\rm{D}}}}{{m}_{\mathrm{imp}}}.\end{eqnarray} \tag{ 9 }$$

Mass conservation requires that ξ_L + ξ_S + ξ_D = 0.

We then derive the orbital–dynamical parameters that are sufficient to set up the postimpact bodies on their new orbits following an identified collision in an N-body simulation. In this step, we use the radius of an equivalent body with the same mass and core mass fraction obtained from 1D calculations with the initial thermodynamic profile, using the results shown in Appendix A. This is to use the same radii as in the N-body. Using the postcollision radii would either create additional complexity or lead to the use of inconsistent radii. Having a consistent radii is important to ensure that the relative velocity and impact parameter are correct.

For the collisions with no second remnant, the velocity change (postimpact orbit) of the target is computed, where, to conserve momentum in the center-of-mass frame, it has the equal but opposite momentum of the debris. For the collisions with two major remnants after the initial collision, whether they end up being hit-and-run or graze-and-merge, we determine the final (outgoing) orbits of both the bodies. Since we do not include the preimpact rotation, we can assume the outgoing orbits are in the same plane as the incoming colliding bodies. The true anomaly is not determined or considered herein, as the initial location of the remnants after the collision will be prescribed according to the capabilities and design of N-body codes that will use our results and simulate the evolution of the remnants.

As in Emsenhuber et al. (2020), there are three additional orbital–dynamical parameters from the postimpact remnants that need to be computed to describe their orbits in sufficient detail for N-body simulations. We compute the orbital energy as

$$\begin{eqnarray}&&\epsilon ^{\prime} =\displaystyle \frac{v{{\prime} }_{\mathrm{coll}}^{2}}{v{{\prime} }_{\mathrm{esc}}^{2}}-1=\displaystyle \frac{v{{\prime} }_{\infty }^{2}}{v{{\prime} }_{\mathrm{esc}}^{2}},\end{eqnarray} \tag{ 10 }$$

where the primed (') quantities are computed from the properties of the largest and second remnants instead of the target and impactor, respectively. The second equality is only valid in the case of an unbound orbit. Next, we compute the impact parameter

$$\begin{eqnarray}&&b^{\prime} =\displaystyle \frac{h^{\prime} }{{R}_{\mathrm{coll}}^{{\prime} }{v}_{\mathrm{coll}}^{{\prime} }}\end{eqnarray} \tag{ 11 }$$

where $h^{\prime}$ is the norm of the specific relative angular momentum vector. Finally, we have the shift of the orbit's pericenter

$$\begin{eqnarray}&&{\rm{\Delta }}\varpi =\varpi ^{\prime} -\varpi .\end{eqnarray} \tag{ 12 }$$

Together, these orbital–dynamical parameters are sufficient to set up the postimpact orbits in an N-body simulation.

The "regime" column in Table 2 is a flag that is automatically set from the analysis of a collision. It determines which properties of the simulation are returned. It is based on the number of significant remnants that are found, where a significant remnant is taken to be a body whose mass is at least 10% that of the impactor; considering the resolution of the simulations and the range of impactor-to-target mass ratios covered, this means that the significant renmants' classification is restricted to bodies of at least ∼500 SPH particles. The core mass fractions are only determined on significant remnants to avoid computing the statistics on bodies that have low numbers of SPH particles in the simulation and thus may be underresolved. To accommodate graze-and-merge collisions, we provide two sets of outcomes for the simulations results: the first for the end state of the simulation and the other after only a single encounter. Depending on the type of collision, these may or may not be at the same time.

The different regime labels are as follows:

• 1.  
  727 entries or 58.2% are hit-and-run collisions where two unbound significant remnants exist. Per our definition of a significant remnant, this includes erosive hit-and-runs. Here, all the values are computed at the end of the simulation.
• 2.  
  218 entries or 17.4% are finished graze-and-merge collisions, where two significant remnants were temporarily gravitationally bound after first contact, and collided again thereafter in an accretion. Here, all the properties are determined first at apocenter and then at the end of the simulation.
• 3.  
  10 entries or 0.8% are unfinished graze-and-merge collisions, where the loopback orbital period is too long to evolve back to apocenter using SPH. Here, only the properties after the initial grazing encounter are determined, from the end state of the simulation.
• 4.  
  281 entries or 22.5% are collisions that result in a single major remnant. Here, the accretion efficiency and core mass fraction are determined only for the largest remnant, at the end of the simulation.
• 5.  
  14 entries or 1.1% are fully erosive events where no significant remnant exists. Only the accretion efficiencies (negative) are determined in this case, at the end of the simulation, as there is insufficient resolution to determine the other properties.

We observe that hit-and-runs are the most represented collisions in our database. This is consistent with the fact that these collisions are found for a wide range of impact angles for our favored impact velocities (Cambioni et al. 2019). Only 14 collisions (or 1.1% of the total) are supercatastrophic. This is unsurprising as they only occur at large velocities and low impact angles (Leinhardt & Stewart 2012), which our selected velocity distribution does not favor.

For unfinished grazing collisions—those labeled 3—we refrain from providing the final values. The final outcome is less relevant in the context of planetary formation, where the loopback orbit would be perturbed by third bodies, such as the Sun (Emsenhuber & Asphaug 2019b).

Due to the sensitivity on the impact angle and velocity, grazing collisions can have a diversity of specific outcomes even for similar initial conditions. Example end states include the capture of a substantial part of the runner as a satellite, or producing bodies that are far from their original hydrostatic equilibrium. We provide a selection of these outcomes in Figure 4 that will be discussed in more detail in the following sections.

Zoom In Zoom Out Reset image size

Some of the peculiar outcomes of grazing collisions are the situations where most of the impactor (the runner) ends up forming a satellite orbiting the largest remnant. This is akin to a proposed scenario for the formation of the Pluto–Charon binary (Canup 2005) and some scenarios for the formation of our Moon (Benz et al. 1987). We remind the reader that a bound satellite—like any bound material—is considered part of the largest remnant in this study.

We also examine the long-term evolution of the satellite. If the satellite orbit is below the corotation radius, tidal friction will transfer angular momentum over time from the satellite's orbit to the spin of the central body, and it will spiral down on a timescale dictated by tidal dissipation and reimpact the central body. This results in the accretion of the majority of the satellite's mass and ejection of a fraction of the mass to release angular momentum. If the satellite is captured outside the corotation radius, then the opposite occurs, and angular momentum is transferred to the satellite's orbit, and it migrates farther out.

Examples of a satellite capture in our database are shown in Figure 4, left and center panels. One scenario is a collision with m_tar = 1.154×10⁻⁴ M_⊕, and γ = 0.657, and whereas the other is a higher-mass target m_tar = 1.618 M_⊕, with γ = 0.376. Both occur near the mutual escape velocity (v_coll/v_esc = 1.02 and 1.00, respectively) at a grazing incidence (751 and 795, respectively), which is characteristic of satellite capture scenarios.

Figure 5 shows how multiple loopback encounters can lead to satellite capture. The blue thin curve plots the relative distance between the largest remnants versus time for the low-mass case in Figure 4 (left panel). The orange curve shows the same information for the high-mass case (center panel). The gaps in the curves indicate the moment the bodies are in physical contact. An arm (continuous link of material) usually connects the two bodies until the bodies reach a distance of up to several mutual radii. During this time, the friend-of-friend algorithm cannot identify the central body and satellite as separate objects, and the orbit cannot be determined. The results show that the pathway toward capture differs between the two collisions.

Zoom In Zoom Out Reset image size

For the low-mass case, each close encounter leads to dissipation in the satellite, which decreases the apocenter. A little orbital angular momentum is transferred onto the target, likely due to the inclusion of material strength in our simulations, which limits its deformation. Consequently, the orbit becomes more circular. For the high-mass case, both the satellite and target are heavily deformed during each loopback encounter, and much of the impactor's material is transferred onto the target during successive encounters. This raises the pericenter when tidal deformation induces a torque on the remnant of the impactor, which increases its angular momentum.

In the high-mass scenario in Figure 5 (orange lines), it is the third encounter at ∼23 hr that is the most effective at establishing a relatively stable orbit for the impactor. The satellite's orbit at the end of our simulation is eccentric (see difference between thin curve and thick curve, representing radius and semimajor axis of the orbit respectively), resulting in the satellite spending most of the time many radii away from the target. On its subsequent close approaches, tidal destruction is avoided due to the sufficiently high pericenter. We evolve these simulations long enough for the satellite to make multiple pericenter passages to check that satellites are not destroyed during this phase.

For all capture scenarios, we report the mass ratio of the orbiting pairs in Figure 6. For this analysis, we ensure that there is at least one pericenter passage of the loopback orbit. We find that the satellites with more than ∼10% of the mass of the primary occur around the primaries with masses below 2 ×10⁻³ M_⊕. Friction likely aids in the capture of such massive moons, as exemplified by the low-mass capture scenario shown in Figures 4 and 5 where the satellite remains at about twice the mutual radii. Without friction, such a satellite would likely be disrupted by tidal forces, as it lies inside the Roche limit of a fluid body. Satellites are obtained around more massive bodies, but are of a small mass ratio with respect to the postimpact target. This outcome is consistent with prior work (Nakajima et al. 2022) on the limited ability of larger planets to form fractionally massive moons; however, the result here is for collisional capture and is a direct aspect of the giant impact phase itself. Also, as a general rule, fractionally massive surviving clumps are less common at larger scales of giant impacts because of the greater relative energies involved (Asphaug & Reufer 2014).

Zoom In Zoom Out Reset image size

An outcome of low-mass, graze-and-merge collisions is the emergence of nonspherical bodies. Friction plays a role in sustaining shapes that are far from hydrostatic equilibrium (e.g., Tanga et al. 2009; Sugiura et al. 2018; Jutzi et al. 2019). For example, in Mars-sized planets and smaller, friction forces can counteract gravity, allowing for the stranding of impact cores in mantle (Emsenhuber et al. 2018). On geologic timescales, these nonhydrostatic shapes and irregular core structures might relax, lowering the moment of inertia and changing the spin state.

An example is shown in the third panel of Figure 4, which shows a cross section of simulation 58 of the 1000 simulation set. These scenarios are analogous to hypotheses for the formation of contact binaries such as 486958 Arrokoth (Stern et al. 2019) and 67P Churyumov-Gerasimenko (Jutzi & Asphaug 2015; Jutzi & Benz 2017), although those scenarios invoke crushable solid bodies while our models are at much larger scales and consider bodies that start out with metallic cores. Here, it can be seen that the two bodies only slightly mix, and the cores remain separated.

The nonspheroidal shape of the final body, in these extreme cases, is maintained not only by material strength but also from the spin induced by the collision. The production of nonspherical shapes generally requires grazing mergers at velocities close to v_esc. In these cases, the two cores do not immediately intersect during the initial collision, yet there is enough dissipation in the initial encounter that the bodies remain bound, rotate quickly, and maintain a spin-supported shape. A caveat of our analysis, however, is that bodies are not spinning before the collision, so preimpact spin is not taken into account.

We detected bodies with nonspherical shapes by computing the moment of inertia of all massive remnants and comparing the values of an equivalent nonrotating spherical body of the same composition under a hydrostatic equilibrium. Figure 7 is a scatter plot of the ratio between these values as a function of mass. Ratios close to 1 indicate bodies whose internal structure is comparable to that of initial bodies.

Zoom In Zoom Out Reset image size

The analysis reveals several features. The remnants from hit-and-run collisions, where ξ_L ≈ 0 (in gray in the figure), in the end have their global structures barely affected by the collision. Accretionary or erosive collision usually leads to an increased moment of inertia, due to spin-induced deformation (equatorial bulge), for accretionary events, or pressure release, for erosive events. Finally, we see two regions that exhibit moment of inertia factors larger than 2: one for small bodies (below ∼2 × 10⁻³ M_⊕) and another for large bodies (around 1 M_⊕).

For small bodies below ∼2 × 10⁻³ M_⊕ (or 2000 km diameter) for the modeled terrestrial compositions, the high moment of inertia factor reflects the peculiar shapes. In Figure 4 (right panel), the contact binary would plot toward the top left of Figure 7. It is not a unique case in our database; however, this compound-core scenario only happens for small bodies. The specific compound-core outcome and other aspects of shape depend on the influence of friction, especially in the silicate material that supports the denser core material, where for consistency with past and ongoing research we have adopted the coefficients of friction in Table A1 of Collins et al. (2004). The upper limit of ∼2 ×10⁻³ M_⊕, which coincides with the limit for satellites with large mass ratios (Section 3.2), suggests that material strength has global-scale effects up to that mass. We note that this upper limit depends on the constitutive strength model and its parameters. As a consequence, the bodies with significantly higher initial temperature or lower von Mises plastic limit Y_m would exhibit the transition to more spherical shapes at smaller sizes.

In the case of the accretionary events in Figure 7, the increased moment of inertia factor is due to a combination of two effects. The predominant one is the spin, which causes a rotational bulge. Unlike the peculiar shapes we just discussed, however, the shape of these bodies is nearly symmetric about their spin axis. Other effects, limited to the masses above ∼ 1 × 10⁻² M_⊕, are thermal in nature. This is because mergers, especially those between similar-mass bodies (that is, for γ ≈ 1), release the most binding energy (Asphaug & Reufer 2013; Carter et al. 2020), increasing temperature. These bodies will also preferentially form distended vapor-rich outer structures due to shock heating in the largest remnant (Lock & Stewart 2017). This in turn expands the bodies relative to a starting configuration (Lock & Stewart 2017) and increases their moment of inertia.

Strongly erosive collisions (ξ_L ≈ − 1) also exhibit similar expanded states because of pressure release, as in these cases usually only the most central part of the body remains. Even-lower values of ξ_L are obtained in supercatastrophic disruptions, and show the same effect, but they are not plotted because the bodies are comprised of relatively few particles.

One of the main objectives of our new simulations is to better characterize the debris to represent it with known accuracy in planet formation studies. For this, we need to understand not only debris production efficiency ξ_D but also its velocity and the size distribution.

3.4.1. Fragment Sizes and Number

Resolvable fragments are clumps that are obtained in the debris field, identified in the same way as we identify the largest and second remnants, by performing friend-of-friend and gravity searches in that order. The smallest mass that can be resolved depends on the numerical resolution of the simulations; however, the size–frequency distribution as a whole is much less sensitive to numerical resolution (Appendix C.2). We assume that a fragment requires at least ten particles to be meaningful, for our database with roughly 1×10⁵ particles per simulation, and thus, only the clumps larger than ∼1×10⁻⁴ the total mass are regarded as significant. The particles in the smallest clumps have less than a full neighbor-list (approximately 50 particles), and while this means that SPH forces are not computed accurately, the self-gravity that causes clumping is accurate. We have therefore verified that the results discussed here remain if we change the fragment resolution limit from 10 to 50 to 100 SPH particles.

Fewer than half of our modeled collisions produce fragments 10 SPH particles in size or larger, besides the two largest remnants. No debris fragment is larger than 10% of the total mass, and only a few collisions produce any fragment larger than 1% of the total mass. The only exception is the second remnant, which we do not include in the debris size–frequency distribution. Most massive fragments are found in either graze-and-merge or disruptive hit-and-run collisions, where the impactor is broken in multiple fragments.

It is extremely difficult for N-body simulations to track debris directly, as the number of bodies (N ) would increase dramatically whenever there is a collision. Thus, debris need to be treated differently from the two largest remnants, with a more statistical approach, such as tracer particles (e.g., Levison et al. 2012; Walsh & Levison 2019), bodies of a fixed mass (e.g., Chambers 2013), or as a background surface density (e.g., Carter et al. 2015). Possible statistical approaches include looking for correlations between debris size and velocity (Section 3.4.2) or fitting a size–frequency distribution to the bodies in the debris field (as in Appendix C). Both require that at least a few clumps are reliably determined, and we thus investigated the conditions that are favorable for the production of many resolved debris clumps. We find that a well-populated size distribution of debris is generated in graze-and-merge and shallow angle, erosive collisions. Graze-and-merge ejects clumps due to the large angular momentum involved (Asphaug & Reufer 2013), and the shallow-angle erosive collisions (v_imp/v_esc ≥ 2) naturally produce ample amounts of debris.

We further find that low-mass bodies (m_tar ≲ 0.1M_⊕) are favored for the production of a large number of resolved debris clumps. We illustrate this in Figure 8, where the number of fragments for each collision is shown as functions of the target's mass m_tar and the relative collision velocity v_coll/v_esc. Here, we observe that there is a strong transition at v_coll ≈ 20 km s⁻¹, which corresponds to m_tar = 0.1 M_⊕ for v_coll/v_esc ≈ 4. While the maximum number of fragments identified in any collision is larger than 400, none of the collisions above 20 km s⁻¹ has more than 40 resolvable fragments. At even larger velocities (v_coll ≈ 50 km s⁻¹), only few collisions report any resolvable fragment; nearly all the debris are found in unresolved particles. This effect could be due to stronger shocks or vaporization being more common at a larger absolute velocity, both of which can prevent the formation of debris clumps.

Zoom In Zoom Out Reset image size

Our debris results should be taken with caution. First, we use a constant number of particles for our simulations, which means that the minimum mass of resolved fragments scales with the target mass. Therefore, fragments of a given absolute size may be resolved in collisions involving low-mass bodies and not with high-mass bodies. Second, high-velocity collisions involve extreme amounts of heating, where numerical effects in SPH may be exacerbated.

3.4.2. Debris Dynamics

Collision after collision, the fraction of the total mass that becomes debris may become significant with time (Emsenhuber et al. 2020). These debris are reaccreted by the large bodies on timescales that are millions of years (Jackson & Wyatt 2012), although some works suggest that the mass fraction is never significant at any given time (Carter et al. 2015). As such, debris orbital evolution must be modeled consistently during N-body calculations. For this purpose, we report postimpact velocities and their correlation with remnant sizes.

First, we compute the rms of the debris velocity distribution relative to the center of mass, reported as v_rms in Table 2. This value is computed from their specific kinetic energy. Because the debris are in the gravitational potential of the largest remnants, this value is time-dependent. To account for this, we also report the velocity of the debris calculated using the total orbital energy (kinetic and potential) of the debris field (v_RMS,∞).

We then check for any correlation between debris size and their velocity using the Spearman's rank correlation coefficient. We find that, over the 1250 simulations presented here, 591 have a sufficient number of resolved debris such that a correlation can be computed. Of those, 62 have p-values consistent with there being no correlation lower than 1%, indicating that the correlation is statistically significant. Of those, 10 have a positive correlation with debris size and velocity, and 52 have a negative correlation. Thus, only a minority (about 10%) of our simulated collisions with resolved debris fields show correlations between debris size and velocity. In sum, we find that debris velocity and mass are not correlated in a systematic way.

To broadly describe debris velocity fields in our simulations, we compute the mean velocity of the debris at infinity in terms of the relative velocity of the precursor bodies at initial contact (see Figure 9). The left panel shows that the collisions producing the most debris (shown toward the right of the panel) have velocities that trend near a value of about 0.45, with a significant dispersion around the mean due to the other parameters (mainly the mass ratio γ). The collisions that produce a smaller fraction of debris exhibit the greatest dispersion and a lower mean velocity; however, capturing the debris velocity accurately is less important for those events, from the point of view of planet formation, as their debris constitutes a much smaller fraction of the total mass. The right panel of Figure 9 shows that the debris velocity in terms of v_coll is largely independent of v_coll/v_esc of the incoming bodies. This indicates that v_coll is a better reference for the debris velocity rather than, for instance, v_esc.

Zoom In Zoom Out Reset image size

As a rule of thumb, we find that the mean debris velocity at infinity v_RMS,∞ is about half that of the collision velocity v_coll. This can be significantly greater than v_esc of the largest remnant. Thus, the assumption that the debris always leave at near-escape velocities, which is made by Chambers (2013) and other works based on their prescription (e.g., Clement et al. 2019), severely underestimates debris velocity. The consequence is an artificial damping of the velocities of a late-stage accreting planetary system and, as a potential consequence, an artificial limitation to the exchanges of materials between formation regions. The debris velocities we obtain differ significantly from the ejecta velocity distributions of collisions reported in Leinhardt & Stewart (2012); this is understandable for two reasons. First, their results rely on a discrete element model instead of solving the Euler equations, and thus are idealized representations of the collisional physics at the high velocities of planet formation. Second, and more fundamentally, their results are computed for the catastrophic disruption regime (their Figure 5), in which the escaping velocities are governed by impact ejection mechanics rather than gravitational accretion dynamics. Our focus is on the regimes relevant to planetary accretion, so our database barely extends into their regime of catastrophically disruptive collisions (Table 1). Debris are not blasted from the target so much as they are flung out by the complex gravitational interactions of two similar-mass bodies, which explains why the debris velocities in our database are governed by v_esc.

If the relative velocity at infinity v_RMS,∞ is half of v_coll, then for collisions where ${v}_{\mathrm{coll}}/{v}_{\mathrm{esc}}\lt 2\sqrt{3}/3\simeq 1.15$, the debris velocity at infinity v_RMS,∞ can be larger than the relative velocity of the incoming bodies at infinity v_∞. This is not paradoxical, as it happens in the case of strong gravitational torques by the more massive bodies in a close approach, and the subsequent ejection of arms of material. Collisions with such low impact velocities usually produce relatively little debris in any case, as can be seen by the color scale on the right panel of Figure 9.

Several simulations exhibit debris with negative total energy and are not represented in Figure 9. This is confined to weakly interacting hit-and-runs resulting in two large remnants with debris that resides between them in space. These debris experience the gravitational potential of both remnants, which are in opposite directions. Thus, it is not a contradiction that the debris are not bound to either remnant while they have negative energy. These cases result in debris with low relative velocity (lower than that of the runner). This peculiar behavior should not be affecting the overall analysis, as they only occur in simulations with a low-mass fraction of debris. While the major characteristics of the largest remnants are well determined in this study, a more informed assessment of debris will require higher-resolution simulations running to a much later time.

We provide a comparison of the global outcomes of five simulations in Tables 3 and 4. The former represents the outcomes after the first contact (i.e., the first phase of a graze-and-merge collision), and the latter represents the final outcome. The columns provided are the same as those in Table 2, except they include a new parameter, N, for the SPH resolution (particle count), and we removed the outcome n_P since satellites are not produced in these collisions. Furthermore, we did not include the number of resolved clumps N_res in these tables, as we analyze this later in this section.

Table 3. Intermediate Outcomes of Five Representative Collision Scenarios Computed at Two Different SPH Resolutions, N = 1×10⁵ (from the Database) and N = 1×10⁶, Shown after First Contact

Collision Parameters							After First Enc.
mtar	γ	Ztar	Zimp	$\tfrac{{v}_{\mathrm{coll}}}{{v}_{\mathrm{esc}}}$	θcoll	N	ξL	ξS	ξD	ΩL	ZL	IZ,L	IA,L	$\epsilon ^{\prime}$	$b^{\prime}$	Δϖ	ΩS	ZS	IZ,S	IA,S
(M⊕)					(deg)					(rad hr−1)						(rad)	(rad hr−1)
1.72 × 10−4	0.961	0.371	0.284	1.169	51.2	1 × 105	0.004	−0.004	0.000	0.840	0.369	1.127	1.043	−0.196	0.892	−0.687	0.834	0.281	1.135	1.047
1.72 × 10−4	0.961	0.371	0.284	1.169	51.2	1 × 106	0.004	−0.004	0.000	0.734	0.368	1.114	1.036	−0.147	0.881	−0.640	0.717	0.281	1.111	1.037
7.14 × 10−4	0.319	0.192	0.294	1.837	30.3	1 × 105	0.209	−0.321	0.112	0.844	0.200	1.107	1.039	0.334	0.719	−0.979	1.579	0.312	1.174	1.063
7.14 × 10−4	0.319	0.192	0.294	1.837	30.3	1 × 106	0.208	−0.354	0.146	0.764	0.201	1.108	1.032	0.350	0.699	−0.994	1.239	0.324	1.120	1.045
1.01 × 10−3	0.052	0.398	0.130	1.026	51.7	1 × 105	0.464	−0.469	0.005	0.318	0.393	1.024	1.016	−0.261	0.940	−1.030	6.720	0.160	11.162	7.585
1.01 × 10−3	0.052	0.398	0.130	1.026	51.7	1 × 106	0.412	−0.414	0.002	0.273	0.394	1.013	1.009	−0.230	0.882	−1.035	10.048	0.158	20.187	13.636
4.80 × 10−2	0.812	0.731	0.270	2.225	12.5	1 × 105	−0.081	−0.764	0.844	0.713	0.725	1.049	1.214	−0.392	0.245	−2.556	1.023	0.623	7.755	5.539
4.80 × 10−2	0.812	0.731	0.270	2.225	12.5	1 × 106	−0.145	−0.725	0.870	0.591	0.724	0.954	1.241	−0.423	0.262	−2.538	0.632	0.626	1.488	1.273
3.42 × 10−1	0.458	0.542	0.520	5.210	23.5	1 × 105	−1.249	−1.000	2.249	0.018	...	1.853	1.733	...	...	...	...	...	...	...
3.42 × 10−1	0.458	0.542	0.520	5.210	23.5	1 × 106	−1.263	−1.000	2.263	0.040	...	1.686	1.636	...	...	...	...	...	...	...

Note. Intermediate outcomes are especially useful for comparing the graze-and-merge cases where the bodies recollide later in the simulation. Each group of two rows represents one scenario simulated at both resolutions: graze-and-merge, hit-and-run, another graze-and-merge, near-head-on erosive, and slightly off-axis highly erosive/disruptive.

Download table as:  ASCIITypeset image

Table 4. The Final Outcomes of the Same Five Collision Scenarios in Table 3, Again as a Function of SPH Resolution

Collision Parameters							End State
mtar	γ	Ztar	Zimp	$\tfrac{{v}_{\mathrm{coll}}}{{v}_{\mathrm{esc}}}$	θcoll	N	ξL	ξS	ξD	ΩL	ZL	IZ,L	IA,L	ξP	γP	$\tfrac{{v}_{\mathrm{RMS}}}{{v}_{\mathrm{esc}}}$	$\tfrac{{v}_{\mathrm{RMS},\infty }}{{v}_{\mathrm{esc}}}$
(M⊕)					(deg)					(rad hr−1)
1.72 × 10−4	0.961	0.371	0.284	1.169	51.2	1 × 105	1.000	−1.000	0.000	4.366	0.326	2.474	1.801	1.000	0.000	...	...
1.72 × 10−4	0.961	0.371	0.284	1.169	51.2	1 × 106	1.000	−1.000	0.000	4.367	0.326	2.394	1.746	1.000	0.000	...	...
7.14 × 10−4	0.319	0.192	0.294	1.837	30.3	1 × 105	0.209	−0.321	0.112	0.844	0.200	1.107	1.039	0.209	0.203	0.838	0.828
7.14 × 10−4	0.319	0.192	0.294	1.837	30.3	1 × 106	0.208	−0.354	0.146	0.764	0.201	1.108	1.032	0.165	0.179	0.956	0.947
1.01 × 10−3	0.052	0.398	0.130	1.026	51.7	1 × 105	0.917	−0.928	0.011	0.534	0.389	1.040	1.018	0.875	0.004	0.561	0.543
1.01 × 10−3	0.052	0.398	0.130	1.026	51.7	1 × 106	0.894	−0.942	0.048	0.506	0.388	1.032	1.011	0.836	0.003	0.357	0.320
4.80 × 10−2	0.812	0.731	0.270	2.225	12.5	1 × 105	0.278	−0.953	0.675	1.559	0.696	1.410	1.268	0.206	0.032	1.050	1.078
4.80 × 10−2	0.812	0.731	0.270	2.225	12.5	1 × 106	0.248	−0.923	0.675	1.570	0.697	1.236	1.176	0.170	0.052	1.134	1.127
3.42 × 10−1	0.458	0.542	0.520	5.210	23.5	1 × 105	−1.249	−1.000	2.249	0.018	0.868	1.853	1.733	−1.249	0.000	2.398	2.395
3.42 × 10−1	0.458	0.542	0.520	5.210	23.5	1 × 106	−1.263	−1.000	2.263	0.040	0.868	1.686	1.636	−1.263	0.000	2.403	2.400

Download table as:  ASCIITypeset image

The simulations for this resolution study were selected to represent a variety of outcomes from the overall database. Most of the metrics that we are studying are similar whether the number of particles N is 1×10⁵ or 1×10⁶. The perfect merger collision, shown in the first two rows in each table, shows nearly identical results in all outcomes as a function of resolution. For the others, the accretion efficiencies show differences in the range of 0.01–0.03, which can equate roughly to less than a few percent of total mass depending on the impactor mass. The final rotation rate and core mass fraction of the largest remnant are also very similar across resolution in the final state. For the intermediate state, the rotation rate of the largest remnant seems to be somewhat variable, showing ∼20% relative deviation in some cases. The debris velocities are also more uncertain, in two cases with differences of up to 0.2v_esc. These differences, however, do not generally exceed the statistical variations of these quantities in the database under small variations in the impact angle or velocity. Additionally, as noted, in cases where the core mass fractions of the initial bodies are very large, the mantle materials are less well resolved, and thus may be subject to SPH challenges resolving material discontinuity across the core–mantle boundary. While these cases might be expected to show greater variations as a function of N, in our resolution comparisons with the largest core mass fractions, we do not observe large deviations in the composition or other properties of the remnants. Thus, while the higher-resolution studies produce some differences, which are important to explore further as resources permit, overall, we conclude that our findings are robust at the N = 1×10⁵ resolution of the database.

To check for the robustness of the size–frequency distribution with resolution, we selected one simulation where several hundreds of clumps are produced as debris: simulation 566, with m_tar = 4.80×10⁻² M_⊕, γ = 0.81, v_coll/v_esc = 2.23, and θ_coll = 125. The size–frequency distributions at the end of the evolution for two simulations, one from the database (labeled "100k") and a similar one but with 1 million SPH particles instead (labeled "1M"), are shown in Figure 12. Overall, the two distributions appear to match quite well; the mass of the largest remnant is similar in both simulations, and the total number of clumps at the resolution limit of the nominal simulation is also similar. There are nevertheless a few differences at intermediate remnant masses. The high-resolution simulation produces a more massive second remnant, as expected (e.g., Genda et al. 2015). The high-resolution run also only has two remnants comparable in mass to the second-largest, whereas the nominal-resolution simulation has three.

Zoom In Zoom Out Reset image size

To understand whether debris masses and velocities are influenced by common choices in SPH resolution used in the giant impact studies, we perform a simulation at 10⁶ particles and compare it to the nominal resolution in the data set (10⁵ particles). We select a run from the database where debris represents roughly 70% of the total mass, and where no clumps are found (simulation 561, m_tar = 3.42×10⁻¹ M_⊕, γ = 0.46, v_coll/v_esc = 5.21, θ_coll = 235). As shown in Figure 13, the amount of debris produced and their velocities are very similar in both simulations; the debris are traveling at roughly 2.5 times the mutual escape velocity of the colliding bodies, or nearly half of their relative velocity at initial contact. In addition, no clumps (minimum of 10 particles) are found within the debris even when using 1 million particles, so the debris size remains unresolved.

Zoom In Zoom Out Reset image size