THE FREQUENCY OF GIANT IMPACTS ON EARTH-LIKE WORLDS

For each collision event throughout an N-body simulation, the specific impact energy in the center of mass reference frame, Q, can be computed by

$$\begin{eqnarray}&&Q=\displaystyle \frac{\mu \ {v}_{\mathrm{imp}}^{2}}{2\ ({m}_{\mathrm{imp}}+{m}_{\mathrm{tar}})},\end{eqnarray} \tag{ 1 }$$

where ${m}_{\mathrm{imp}}$ is the mass of the impactor, ${m}_{\mathrm{tar}}$ is the mass of the target, μ is the reduced mass, where

$$\begin{eqnarray}&&\mu =\displaystyle \frac{{m}_{\mathrm{imp}}\ {m}_{\mathrm{tar}}}{{m}_{\mathrm{imp}}+{m}_{\mathrm{tar}}},\end{eqnarray} \tag{ 2 }$$

and ${v}_{\mathrm{imp}}$ is the impact velocity.

Leinhardt & Stewart (2012) performed a large number of high-resolution two-body collision simulations to explore the outcomes using a wide range of initial conditions, including impactor mass (${m}_{\mathrm{imp}}$), target mass (${m}_{\mathrm{tar}}$), impact angle θ, and impact velocity ${v}_{\mathrm{imp}}$. By comparing Q with the criteria for catastrophic disruption between two colliding bodies, ${Q}^{*}$ (which by definition is the specific impact energy needed to disrupt half of the total mass), and evaluating the outcomes for all of their collision simulations, they derived analytic equations to map out the transitions between collision regimes. For a given collision, this model provides the outcome (for example, a grazing impact with partial erosion) as well as the size and velocity distributions of the post-impact bodies.

Chambers (2013) implemented these collision outcomes into Mercury and introduced a new input parameter m_frag that defines a minimum fragment mass for all collisions. In this modified version of Mercury, if a collision leads to fragmentation, the mass of the largest remnant remaining from the target, m_lr, is recorded. The remaining mass from the target (or the impactor in some scenarios), m_fragtot, is broken up into x equal-sized fragments, where x = m_fragtot/m_frag (with any remainder mass less than m_frag left as is). The fragments are gravitationally dispersed at slightly above the escape velocity in uniform directions within a plane that includes the center-of-mass of the system. The evolution of these new fragments is then followed from that time forward while they gravitationally interact with other bodies in the disk. Note that the collision model of Leinhardt & Stewart (2012) does not describe the fragments as a set of equal size bodies. The simplification of equal size fragments implemented in Mercury is intended to keep simulations computationally tractable, as the true fragments would have a size distribution with many smaller fragments.

The outcome regimes of two-body collisions in our simulations with fragmentation include:

• 1.  
  A collision with the central star (or outer giant planet) in which case the body is removed from the simulation.
• 2.  
  Ejection from the system if a body orbits beyond 100 au from the central star.
• 3.  
  Perfect accretion (inelastic collisions) where the mass of the new body m_lr is equal to ${m}_{\mathrm{imp}}+{m}_{\mathrm{tar}}.$
• 4.  
  A grazing collision that can lead to accretion or erosion. Grazing collisions occur when less than half of the projectile interacts with the target. In this case the impact parameter b is greater than the critical impact parameter b_crit = r_tar/(r_tar + r_imp), where r_tar is the radius of the target and r_imp is the radius of the impactor.
• 5.  
  A non-grazing collision, where b < r_tar, that can lead to accretion or erosion.
• 6.  
  A "hit-and-run" collision, discussed in detail in Asphaug et al. (2006) and Genda et al. (2012), in which the target does not gain or lose a significant amount of mass, but the impactor may suffer fragmentation, so m_lr = ${m}_{\mathrm{tar}}$ and m_fragtot$\quad \leqslant$ ${m}_{\mathrm{imp}}$.

At each collision event the initial and final masses, positions, and velocities are recorded for the target, impactor, and any fragments. In addition, parameters to characterize the collision are output such as the fate (type of collision), impact angle, and impact velocity. All of the information needed to evaluate the collision history for planets that form in an N-body simulation is made available in the output of Mercury.

Our simulations begin at the epoch of planet formation in which runaway and oligarchic growth have already taken place (Kokubo & Ida 1998), resulting in a bimodal disk mass distribution around the Sun. In total, 26 embryos and 260 planetesimals reside in our initial disk. The embryos are approximately Mars-sized (r_em = 0.56 ${R}_{\oplus }$ $\approx$ 3500 km; m_em = 0.093M${}_{\oplus }$) and are spread among approximately Moon-sized planetesimals (r_em = 0.26 ${R}_{\oplus }$ $\approx$ 1600 km; m_em = 0.0093 M${}_{\oplus }$). All bodies are assumed to have a density of 3 ${\rm{g}}\;{\mathrm{cm}}^{-3}$ to determine their radii and the bodies are placed between 0.35 and 4 au from the Sun. The surface density of solids varies as ${a}^{-3/2}$, consistent with Solar Nebula models (Weidenschilling 1977), yielding a total disk mass of 4.85 M${}_{\oplus }$. The bodies begin on nearly circular and coplanar orbits with all other orbital elements (argument of periastron, longitude of ascending node, and mean anomaly) selected at random. At the start of our simulations we include a Jupiter-mass planet at 5.2 au and a Saturn-mass planet at 9.6 au from the Sun at their present orbits and assume that these giant planets have already formed and that any gas in the disk has been dispersed. The problem is then purely gravitational and collisional, and any damping of the eccentricities and inclinations of embryos arises from dynamical friction from the swarm of smaller bodies and not by gas drag.

We implemented the same initial disk model as was used by Quintana & Lissauer (2014), in which several dozen simulations were performed with the perfect accretion model to examine the role of Jupiter and Saturn on water delivery to the terrestrial planets. This disk model was adapted from Chambers (2001) (which was also used by Chambers (2013) to test the fragmentation model), but we extrapolated our disk out from 0.35 to 4 au, whereas the disk in Chambers (2001, 2013) was truncated at 2 au. The sizes of embryos and planetesimals remained the same in all studies, but there were only 14 embryos and 140 planetesimals in Chambers (2001) compared to our 26 embryos and 260 planetesimals.

Because N-body simulations are highly stochastic, we examine a large number of simulations in order to interpret our results statistically. We performed 280 N-body simulations—140 with fragmentation and 140 without—of accretion from the disks of bodies around a Sun-like star with Jupiter and Saturn analogs. The same sets of elements for all bodies in the initial disk were used in each simulation apart from minuscule variations. Starting with the second simulation, we iteratively displaced a planetesimal near 1 au by 1 m along its orbit, thereby creating a set of 140 initial disks that were all (slightly) different. In each simulation the evolution of all bodies in the disk was followed for 1 Gyr (note that in Section 5 we extended the fragmentation simulations to 2 Gyr).

An additional parameter that is needed for the integration is the minimum fragmentation mass, which ultimately constrains the total number of bodies in a given simulation. We set the minimum fragmentation mass to be 0.0047 M${}_{\oplus }$, which is r_frag = 0.2 ${R}_{\oplus }$ $\approx$ 1300 km, which is a little under half a lunar-mass. This is the same value used in Chambers (2013) and was adopted here for comparison. While smaller minimum fragmentation mass values would allow a better representation of the protoplanetary disk that formed our solar system planets, this choice allowed us to constrain the total number of bodies in order to perform a large number of simulations. With this choice of m_frag, the average number of fragments created in our simulations was 270, comparable to the number of initial embryos and planetesimals. The sensitivity of the results to this parameter are discussed in Section 5.3.

To examine the collisions experienced by the Earth-analogs formed in our models, we first determined how to quantify the impacts. Stewart et al. (2014) performed a suite of 3D hydrodynamic simulations of two-body gravitational collisions to study the types of collisions that could blow off a portion or all of an atmosphere or ocean, assuming the bodies involved had accreted an atmosphere and/or ocean. In their study, when the two bodies collide, the shock waves propagate and affect their atmospheres. Using a wide range of target and impactor masses, impact speeds, and impact angles, Stewart et al. (2014) found that variations in the shock pressure fields, and therefore the fraction of atmosphere loss, are highly sensitive to the impact parameter and mass ratio of the impactor to the target. Using the results from these hydrodynamic simulations, S. Lock & S. Stewart (2016, in preparation) empirically derived a modified version of the specific impact energy Q (as defined in Equation (1)), designated Q_s, that better encapsulates variations in the shock pressure field and includes only the interacting mass of the two bodies:

$$\begin{eqnarray}&&{Q}_{s}=Q^{\prime} \left(1+\displaystyle \frac{{m}_{\mathrm{imp}}}{{m}_{\mathrm{tar}}}\right)(1-b).\end{eqnarray} \tag{ 3 }$$

At the point of contact, assuming the target is stationary, the impactor (the bluebody in Figure 7) is traveling with impact velocity ${v}_{\mathrm{imp}}$ toward the target body (shown in brown). The impact parameter is defined as b = sin θ and the projected distance between the center of the two bodies at contact is B = (${r}_{\mathrm{imp}}$ + ${r}_{\mathrm{tar}}$) b, where ${r}_{\mathrm{imp}}$ and ${r}_{\mathrm{tar}}$ are the radii of the impactor and the target, respectively. The projected length of the impactor that intersects the target is l, and depending on whether the impactor fully interacts with the target, can take a value of

$$\begin{eqnarray}l=\left\{\begin{array}{ll}2\ {r}_{\mathrm{imp}} & {\rm{if}}\;{\rm{}}B+{r}_{\mathrm{imp}}\leqslant {r}_{\mathrm{tar}}\\ {r}_{\mathrm{imp}}+{r}_{\mathrm{tar}}-B & {\rm{if}}\;{\rm{}}B+{r}_{\mathrm{imp}}\gt {r}_{\mathrm{tar}}.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

These two cases are shown in the upper and lower panels of Figure 7.

Zoom In Zoom Out Reset image size

The specific impact energy, Q_s, is a function of the ratio of the impactor mass (${m}_{\mathrm{imp}}$) to the mass of the target (${m}_{\mathrm{tar}}$), the impact parameter b, and ${Q}^{\prime }$, which is a modified version of Q that takes into account only the interacting mass of the two bodies. As shown in S. Lock & S. Stewart (2016, in preparation), ${Q}^{\prime }$ can be calculated by scaling Q by the ratio of the reduced mass for the interacting mass, ${\mu }_{\alpha }$, to the reduced mass μ as defined in Equation (2)

$$\begin{eqnarray}&&Q^{\prime} =\displaystyle \frac{{\mu }_{\alpha }}{\mu }Q.\end{eqnarray} \tag{ 5 }$$

The interacting reduced mass is

$$\begin{eqnarray}&&{\mu }_{\alpha }=\displaystyle \frac{\alpha \ {m}_{\mathrm{imp}}\ {m}_{\mathrm{tar}}}{\alpha \ {m}_{\mathrm{imp}}+{m}_{\mathrm{tar}}},\end{eqnarray} \tag{ 6 }$$

where α is the fraction of the impactor mass that interacts and is defined as

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{3\ {r}_{\mathrm{imp}}\ {l}^{2}-{l}^{3}}{4\ {r}_{\mathrm{imp}}^{3}}.\end{eqnarray} \tag{ 7 }$$

The original and full derivations of these impact geometries are given in Leinhardt & Stewart (2012) and S. Lock & S. Stewart (2016, in preparation).

To associate Q_s to physical events on planets, Stewart et al. (2014) showed that, for an atmosphere-to-ocean ratio of 1:300, an impact with Q_s > 10⁷ J kg⁻¹ was needed to strip 100% of an atmosphere, and an order of magnitude greater energy (Q_s > 10⁸ J kg⁻¹) was needed to remove an entire ocean. Evaluating Q_s for collisions in previously published N-body simulations (ones that assumed perfect accretion), Stewart et al. (2014) found that impacts that could cause partial atmospheric loss were common but impacts with energies sufficient to blow off an entire atmosphere or ocean were rare. Note that the presence of an ocean can significantly enhance the loss of atmosphere during a giant collision due to both evaporation, which can push out an atmosphere, and a lower shock impedance of the ocean compared to a rocky surface (Genda & Abe 2005).

Herein, we adopt Q_s as our parameter of interest to evaluate the collisions in our fragmentation simulations, particularly those that lead to the formation of an Earth-analog. Using the Stewart et al. (2014) empirical relations (that considered the 1:300 atmosphere:ocean case), we chose the value of Q_s that is sufficient to remove about half of a planet's atmosphere to define a "giant impact" threshold. Collisions with Q_s that reach or exceed this value, $2\times {10}^{6}$ J kg⁻¹ (equivalent to ${\mathrm{log}}_{10}\;{Q}_{s}=6.3$), are considered a giant impact for our purposes of evaluating the frequency of giant impacts onto the growing terrestrial planets.

As shown in Stewart et al. (2014), the values of Q_s for collisions that occur during the late stages of terrestrial planet formation span a wide range that overlap with the estimated Q_s values from various Moon-formation models. Leading theories for the formation of the Moon suggest that a giant collision occurred between the proto-Earth and another object, resulting in a remnant disk from which the Earth's Moon formed (Canup & Asphaug 2001; Canup 2004, 2008, 2012; Ćuk & Stewart 2012). Based on estimates of the age of the lunar melt, this collision likely occurred some time between 30 and 110 Myr after the solar system began to form (Taylor 1975; Halliday 2008); however, it could have occurred later (see Section 2). Nonetheless, the Moon-forming giant impact was likely the final significant giant impact during the Earth's formation.

The full range of Q_s for the canonical Moon-forming impact models is about 0.3 × 10⁶–2 × 10⁶ J kg⁻¹ (Canup & Asphaug 2001; Canup 2004). For models that require higher angular momenta (Canup 2012; Ćuk & Stewart 2012), the range of Q_s is about 2 × 10⁶–2.5 × 10⁷ J kg⁻¹ (Stewart et al. 2014). Our choice of Q_s to define a giant impact therefore corresponds to significant (∼50%) atmosphere loss and also to the minimum impact energy from the more recent high angular momentum Moon-forming impact models (Canup 2012; Ćuk & Stewart 2012). Our results on the frequency and timing of giant impacts presented herein can be compared with models of the Moon's formation and can provide insight into the frequency of these types of impacts during late stage planet formation. However, our results are not meant to predict the abundance of actual moon-forming systems, as many impacts with energies comparable to the Moon-forming impact may not have enough angular momentum to produce a Moon-forming disk.

We computed Q_s for all impacts that were involved in the formation of an Earth-analog for all of our fragmentation simulations. As illustrated in Figure 8, impacts with enough energy to strip entire atmospheres (${Q}_{s}\geqslant {10}^{7}$ J kg⁻¹) and oceans (${Q}_{s}\geqslant {10}^{8}$ J kg⁻¹) are rare, a result consistent with those of Stewart et al. (2014). The percentage of impacts with energies comparable to the high angular momentum Moon-forming impact theories is just 5% of all impacts, while less than 0.2% of impacts are energetic enough to fully strip an atmosphere of an Earth-analog planet.

Zoom In Zoom Out Reset image size

Figure 9 shows the distribution of these energy values grouped according to three types of impacts: merging events (grazing and non-grazing events that led to net growth), hit-and-run events, and head-on collisions. The majority (64%) of events result in a merger, while 31% are hit-and-run collisions. The remaining 5% of impacts are head-on collisions with high values of Q_s. While head-on collision events are uncommon, they account for all events that are capable of stripping a planet of its atmosphere. For collisions that led to the formation of Earth-analogs in the eight fragmentation runs performed by Chambers (2013), 70% resulted in net growth, 30% were hit-and-run events, and zero collisions led to net erosion. Our results are consistent, considering the relatively small number of simulations and the less massive initial disk used in Chambers (2013).

Zoom In Zoom Out Reset image size

We next computed the total number of giant impacts experienced by each Earth-analog during its formation. All but 7 of the 167 Earth-analogs experienced at least 1 giant impact. A histogram of these results is shown in Figure 10, where the bins with the highest occurrence are at two and three giant impacts. Owing to this being a discrete distribution with a number of events occurring within a fixed time window, we are able to describe the resulting probability mass function as a Poisson distribution. We fit these data of the number of giant impacts and found best fitting parameters of $\lambda =3.04$, which yields the red curve shown in Figure 10. This Poisson distribution function has a mean at 3.0 and a median of 3.4. While each of these giant impacts is likely to deplete a significant fraction of volatile content of the planet (c.f. Stewart et al. 2014), most of these giant impacts are unlikely to fully strip an atmosphere and ocean.

Zoom In Zoom Out Reset image size

The final giant impact received by an Earth-analog largely influences its dynamical state and bulk composition. With the relatively large number of simulations performed in this work, we can make probabilistic inferences from the distributions of the collision times of giant impacts onto Earth-analogs.

In Figure 11 we show the distribution of times of the final giant impact onto an Earth-analog in all of our fragmentation simulations. About 56% of our Earth-analog planets experience their last giant impact within the first 50 Myr of the simulation and 75% experience their final giant impact within the first 100 Myr. This mirrors our understanding of the range of feasible times for the Earth's Moon-forming impact (Figure 1). In terms of final giant impact times, we find that the Earth (assuming the Moon-forming impact was its final giant impact) is relatively typical.

Zoom In Zoom Out Reset image size

The time of the final giant impact onto an Earth-analog is correlated with the abundance of late (post-impact) accreted mass by that Earth-analog, primarily because the amount of solid material available in the disk decreases with time as planets are growing. Jacobson et al. (2014) illustrated this correlation by performing a large set of N-body perfect-accretion simulations using different models for the orbits of the giant planets. For simulations that included giant planets near their current orbits (comparable to our simulations), they found the fraction of mass accreted by an Earth-analog after its final giant impact ranged from about 10%–20% for impacts occurring at 10 Myr, to about 1% for final impacts at around 100 Myr.

From geochemical constraints, primarily the abundance of highly siderophile elements (HSEs) in the Earth's mantle compared to abundances found in chondritic meteorites, the fraction of mass that the Earth accreted after the Moon-forming impact has been inferred to be less than 1% (Chou 1978; Walker 2009). Assuming an estimate of the Earth's late accreted mass of 0.003–0.006 M${}_{\oplus }$, Jacobson et al. (2014) used their correlation to place constraints on the time of the Moon-forming impact to 95 ± 32 Myr (for simulations that included giant planet migration) to a median time of about 110 Myr (for giant planets near their present orbits) after the start of the solar system.

In Figure 12 we show, for all of our fragmentation simulations, a similar correlation between the final giant impact time for an Earth-analog and the fraction of mass the planet accreted after the impact. Compared to Jacobson et al. (2014), our results appear biased toward higher fractions of late mass accreted for a given final giant impact time. A direct comparison is difficult due to the differences in initial conditions and definitions. Jacobson et al. (2014) began with half of the disk mass in about 100 bodies that are Moon-size to Mars-size and half of the disk in thousands of non-mutually interacting Ceres-sized (about 470 km) planetesimals. Our disk begins with half of the mass in 26 Mars-size embryos and the other half in 260 Moon-size planetesimals, and fragments are allowed during the simulation with a minimum resolution of about half a Moon mass. Additionally, our definition for a giant impact relies on a minimum specific impact energy that is physically motivated and more stringent, whereas Jacobson et al. (2014) define a giant impact as a collision between an Earth-analog and a surviving Moon-size or Mars-size embryo (comparable to our initial planetesimals and embryos). Regardless of these differences, Jacobson et al. (2014) assume that all of the late mass accreted has a chondritic composition, which would allow for a direct correlation between the mantle concentration of HSEs and the mass accreted after the Moon-forming impact. Our simulations show that some of the late accreted mass is fragments. Studies of earlier stages of planetesimal formation have shown that planetesimals are differentiated and exhibit a wide range of iron core fractions (Carter et al. 2015). Fragments created in the final stages of planet formation are therefore likely dominated by rocky mantle and should be depleted in HSEs compared to chondrites. As a result, more mass can be accreted post-giant impact than implied by the mantle's HSE abundance, so our results are consistent with the observations and uncertainty in the timing of the Moon-forming impact.

Zoom In Zoom Out Reset image size

While our motivation in this work is not to model and constrain the timing of the Earth's Moon-forming impact, the distribution of mass ratios and impact velocities of the final giant impacts in our simulations may shed some light on the likelihood of a given scenario and may be useful for comparing Moon-formation models. In Figure 13 we show the ratio of the impactor mass to the target mass as a function of the impact velocity (scaled by the mutual escape velocity) for all final giant impacts onto the Earth-analogs that formed in our fragmentation simulations. Most impact events that involve bodies of comparable mass have impact velocities that are 1–2 times the mutual escape velocity, whereas smaller relative impactors reach impact velocities up to 5 times the mutual escape velocity.

Zoom In Zoom Out Reset image size

Here we briefly address some caveats that should be considered when interpreting these (and other similar N-body) results. First, our initial conditions begin with objects that are roughly Moon-sized to Mars-sized. In reality the initial disk would contain these in addition to a significantly larger number of small bodies. In most N-body accretion simulations a much higher fidelity is usually desired but must be balanced with available computational resources.

The choice of the minimum fragmentation mass in our model constrains the total number of objects. We selected a minimum fragmentation mass of a little under 0.5 a lunar-mass, the same value used in Chambers (2013), so we could compare results and also perform a statistically large sample of simulations in a reasonable time. It is expected that the final planets that form, and their collision histories will be different with the selection of a less massive (and more physically justified) minimum fragmentation mass.

Of particular concern in this work was whether our conclusions on the numbers and times of giant impacts are biased by our minimum fragment mass. In Figure 14 we show the specific impact energies of all impacts onto Earth-analogs split into two populations, (a) impactors more massive than one lunar mass and (b) impactors less massive than one lunar mass, where one lunar mass is approximately equal to the initial mass of the planetesimals in our simulations. We found that while only a third of total impactors have a mass of at least one lunar mass, 93% of giant impacts are from bodies at least as large as our initial planetesimal mass. Therefore, our results are unlikely to depend strongly on the initial fragment mass since very few of the giant impacts are from fragments.

Zoom In Zoom Out Reset image size

The initial conditions of the giant planets should be taken into consideration when interpreting the results presented herein. Classical models of the solar system have assumed that the giant planets formed in situ (Chambers 2001, 2013; Quintana & Lissauer 2014). All of our simulations presented herein have adopted this model and include analogs to Jupiter and Saturn near their present orbits. Recent models have suggested that the giant planets may have originated on more eccentric orbits or may have undergone significant migration (Raymond et al. 2009; Walsh et al. 2012; Raymond & Morbidelli 2014). Because these giant planets are the largest perturbers (not counting the central star), their configuration can significantly alter the growth of the terrestrial planets. For this article we focused on the classical theory of planet formation, which includes Jupiter and Saturn at their present orbits. Exploring the effects of various giant planet configurations on the collision history of Earth-analogs would benefit solar system models as well as planetary systems formed around other Sun-like stars in our galaxy.

Finally, we want to stress the importance of considering these processes in a statistical fashion and inferring results from probability distributions. Although planet formation in general is a chaotic process, accretion models have been fine-tuned for decades to reproduce and understand our solar system. Instead, the current architecture of the planets in our solar system should be thought of as one draw from the probability distribution of planetary systems that can be formed from the same protoplanetary disk. It is inference of this "generalized solar system" from which our own planetary system is drawn that simulations such as those performed here are sampling. Caution should be used if drawing strong conclusions based on a small number of simulations.