Onset of Giant Planet Migration before 4480 Million Years Ago

The nature of early bombardment to the inner solar system is a long-standing problem in planetary science. We use the term late accretion to refer to protracted and ongoing bombardment wherein leftover material composed of comets, planetesimals, and asteroids continues to strike the planets after their formation (Table 1). Although a matter of debate, a signal of this bombardment is traced by the chondritic relative abundances of the highly siderophile elements (HSEs, which include the platinum group elements: Os, Ir, Ru, Rh, Pt, and Pd), regarded as pollutants by the earliest phase of late accretion to a planet's mantle (Walker 2009); these elemental abundances are in turn used to infer an exogenous mass augmentation to the terrestrial planets after separation of their cores and mantles (e.g., Bottke et al. 2010).

There exist two competing paradigms to describe the time-dependent flux of late accretion. The first argues for a burst of relatively short duration in the intensity of impact bombardment hailed as a "lunar cataclysm" or, more commonly, "late heavy bombardment" (LHB); this is assumed to have affected the whole inner solar system (Morbidelli et al. 2012, 2018 and references therein; see Ryder 2002) at around 3950 Ma, or perhaps earlier. The LHB is largely implied by interpretations of uncorrected ^40–39Ar chronology from lunar samples (e.g., Fernandes et al. 2013) and has been promulgated by several popular dynamical models involving late-stage giant planet migration (e.g., Gomes et al. 2005). Alternatively, the second concept proposes an overall monotonic decline in the bombardment rate, which is an interpretation most broadly consistent with crater counting studies (Hartmann 1970; Neukum et al. 1975; see Werner & Ivanov 2015 for a recent review).

To transcend the inconsistency between the two different models for late accretion cited above—which also has a strong bearing on the Hadean Earth's proclivity to host a nascent biosphere (e.g., Maher & Stevenson 1988; Fritz et al. 2014)—requires a comprehensive review of geochronological data for the inner solar system and a new investigation of planetary dynamics for the outer solar system.

Direct records of the earliest phases of solar system evolution that could shed light on the inception of late accretion in the first ∼100 Myr are either rare or absent from the crusts of the Moon, Earth, and Mars. On Earth, this deficit has been ascribed to endogenous crust–mantle interactions that recycled the primary crust (e.g., Hopkins et al. 2008, 2010), or less popularly, to the destructive thermal and mechanical effects of impact bombardments that accompanied the earliest stages of late accretion. If intense enough (such as in the event at ca. 4500 Ma attributed to the "late veneer"; Genda et al. 2017), late accretion can result in melting and isotopic rehomogenization of the crust and mantle that compose the planet's silicate reservoirs (Allègre et al. 2008). As a matter of fact, several lines of evidence exist that Earth's terrestrial silicate reservoirs experienced such a re-equilibration event well after cessation of primary accretion (Roth et al. 2014). This is attested by analysis of terrestrial samples using the time-sensitive radiogenic tracers of lithosphere production and destruction such as ¹⁸²Hf–¹⁸²W, ^235,238U–^207,206Pb, ¹²⁹I–¹²⁹Xe, and ¹⁴⁶Sm–¹⁴²Nd (e.g., Albarède & Juteau 1984; Allègre et al. 1995, 2008; Zhang 2002). For the Moon, the mere existence of ancient lunar zircons (zirconium orthosilicate, or Zr(SiO₄)) proves that its crust was not completely melted by endogenous processes such as impacts, at least since ca. 4417 Ma (Nemchin et al. 2009) and likely earlier. Thus far, direct records (as opposed to model ages) of pre-4.42 Ga lunar materials appear to have been destroyed even as evidence accumulates that the Moon experienced ruinous basin-forming events before 4100 Ma (Norman & Nemchin 2014). Indirect evidence from coupled U–Pb and Hf-isotopes also shows that significant differentiation and crust formation on the Moon commenced at ca. 4510 Ma, commensurate with the proposed time of the Moon's own formation (Barboni et al. 2017). Remote sensing studies reveal that Venus's crust may be no older than about 1000 Ma (McKinnon et al. 1997; see Rolf et al. 2018), and the oldest features on Mercury may be comparable in age to those of the Moon and Mars (Strom & Neukum 1988; Marchi et al. 2013). The firmest minimum date for the perpetuation of the crusts of Earth and Mars comes from preserved ages for the oldest terrestrial zircons at ca. 4400 Ma (Valley et al. 2014) and the oldest comparable Martian zircons, which have been dated at 4428–4470 Ma (Humayun et al. 2013; Bouvier et al. 2018). That these minerals survive at all demonstrates that the crusts of the sampled terrestrial worlds did not experience wholesale melting since that time (Figure 2 of Hopkins et al. 2015; Heisinger & Head 2006). Thus, available geological and cosmochemical data from the terrestrial planets, which could directly bear on the temporal flux profile of late accretion, are confined to the last ca. 4450 Myr.

Earth's Moon shows an unusual history: impact-modified lunar samples returned by the Apollo and Luna missions from the lunar near side yield ^40–39Ar, U–Pb, and Rb–Sr ages that seem to cluster around 3950 Ma (Turner et al. 1973; Tera et al. 1974); this is far younger than the oldest documented ages of the lunar crust cited above (Nemchin et al. 2009; Barboni et al. 2017). Such a grouping of ages has long been used to argue that an LHB-type event affected the inner solar system as a rapid rise, and ensuing demise, in the impactor flux (Ryder 1990, 2002). Intriguingly, neither lunar, nor Mercurian, nor Martian crater chronologies show evidence for such an LHB event (Werner 2014), and it has been suggested that the 3950 Ma ages may instead be from a biased sampling compromised from widespread contamination of debris associated with the Moon's Imbrium basin (e.g., Shearer & Borg 2006). Statistical analyses of crater distributions point instead to a monotonic decline in impactor flux from at least 4400 Ma (Figure 1). It is therefore equivocal whether there was an LHB at all, or that if there was, it occurred with as yet no clear temporal constraints (e.g., Hartmann 2003).

Zoom In Zoom Out Reset image size

Beyond Earth, the Moon, and Mars, highly retentive U–Pb and Pb–Pb ages in eucrite meteorites attributed to asteroid 4 Vesta (dubbed HEDs, for the howardite-eucrite-diogenite clan of achondrite meteorites collectively termed vestoids) are known to record a solid crust at 4563 Ma (Ireland & Wlotzka 1992) and subsequent thermal modifications until ca. 4450 Ma (Zhou et al. 2013; Hopkins et al. 2015). However, practically no such ages younger than ca. 4450 Ma exist for any known meteorite class, including the vestoids. Younger ages for meteorite classes are instead largely confined to less retentive ^40–39Ar geochronology data sets (Bogard 1995), which generally display age continua from ∼4450 Ma extending to about 3000 Ma, with occasional resetting events up to the present (e.g., Fernandes et al. 2013).

Despite these discrepancies, the lunar crater record continues its historical role as the standard reference for the nature of late accretion experienced by all solar system bodies (Öpik 1960; Shoemaker et al. 1962). Recent reanalyses, however, of purported ^40–39Ar age clusters at ca. 3950 Ma from lunar rocks are inconclusive as to whether they were caused by an LHB or are the result of a monotonic decline (Boehnke & Harrison 2016). As such, the problem of late accretion reduces to deciding for which  period can the lunar crater record be used to interpret what  primordial impact history connected with the ongoing accretion process after planet formation. Furthermore, how  translatable is the lunar record to other solid surfaces of the solar system? Progress in deciphering this record would greatly enhance our understanding of the formation and evolution of the solar system as a whole and the thermal histories of the crusts of the terrestrial planets; the latter includes Earth's capacity to host a persistent biosphere from the earliest times (Mojzsis et al. 1996; Abramov & Mojzsis 2009; Abramov et al. 2013).

Notwithstanding debate over the existence of an LHB, a suite of dynamical models has been constructed to explain how a discrete and temporary surge in impactors could have transpired at some time in the formative stages of the early solar system (e.g., Gomes et al. 2005; Bottke et al. 2012; Marchi et al. 2012; Morbidelli et al. 2012, 2018; see Nesvorný et al. 2018, Clement et al. 2019). The common theme to all of these models is that they rely on a dynamical instability in the outer solar system. Indeed, support for orbital migration of the giant planets (Jupiter, Saturn, Uranus, and Neptune) in the first few hundred million years comes from several independent lines of evidence stemming from the orbital properties of both these planets themselves and analysis of small-body populations beyond Neptune (see Nesvorný 2018 for a recent review). For instance, Clement et al. (2019) present dynamical simulations that invoke a very early instability shortly after the dispersal of the protoplanetary disk. However, it is important to note that without temporal limits placed by geochronological data, caution is warranted in accepting any such model at face value. Giant planet migration is often, and erroneously, conflated with Grand Tack (Table 1). An early Grand Tack instability naturally reproduces the small mass of Mars (even if there are other models that can accomplish this strict constraint; e.g., Hansen 2009; Walsh et al. 2011; Izidoro et al. 2014; Brasser et al. 2017). Grand Tack is limited to the first million years or so, whereas giant planet migration (Nice model) occurs after that.

Subsequent to primary accretion, the giant planets are thought to have radially migrated because of their gravitational interactions with a massive external population of planetesimals (Fernandez & Ip 1984). The migration of Neptune is chronicled in the orbital distribution of Pluto and the Kuiper Belt objects (Malhotra 1993, 1995). Dynamical models also show that giant planet migration requires planet–planet scattering (Thommes et al. 1999; Brasser et al. 2009; Morbidelli et al. 2010; Levison et al. 2011; Nesvorný & Morbidelli 2012). Giant planet migration scenarios mutually predict that the inner solar system experienced a short but intense surge of cratering by comets sometime in the first 500 Myr (Gomes et al. 2005).

Two points must be emphasized about the nature of giant planet migration:

• (1)  
  the specific trigger for migration—and ensuing bombardment—is still mostly unknown; and
• (2)  
  its timing is inexplicit and could occur after a delay of a few million years to upward of several hundred million years (Levison et al. 2011).

Although debated, clear evidence for the role of outer solar system materials (i.e., comets) in this migration-induced bombardment is lacking on the Moon (Kring & Cohen 2002; Strom et al. 2005; see Greenwood et al. 2011). Dynamical models of terrestrial planet formation all envisage a steady decline in bombardment rate once migration ceases (Neukum et al. 2001; Brasser et al. 2016). This expectation is entirely consistent with the crater chronologies cited above but is at odds with an LHB-like "spike" at any time after about 4400 Ma (Hartmann et al. 2000). For instance, if giant planet migration commenced sometime after 4400 Ma to give rise to a discrete LHB-type event, such an episode should be obvious as a distinct cluster of younger reset ages shared by all inner solar system bodies—including the Moon—in radiogenic systems relatively sensitive to thermal resetting, such as in ^40–39Ar and Rb–Sr (e.g., Bogard 1995; Fernandes et al. 2013). If, however, the migration commenced before 4400 Ma, endogenous crustal processes on geologically active worlds such as Earth and Mars are expected to have erased such evidence. This is not necessarily the case for a small, airless, old and cold body like asteroid 4 Vesta (Zhou et al. 2013), a fact that we use later in this analysis. The classical LHB case for a spike affecting the whole solar system after about 4100 Ma (Ryder 2002) is further undermined by numerous documented pre-4100 Ma (and older) lunar ages (e.g., Petro & Pieters 2006; Zellner 2017). These ages are surmised to reflect partial or complete resetting of high closure temperature mineral geochronological systems (e.g., Norman & Nemchin 2014; Hopkins & Mojzsis 2015; Kelly et al. 2018). A useful tool to interrogate early solar system bombardments is "closure temperature" (hereafter T_c), defined as the temperature of a system at the time given by its apparent age (Dodson 1973; Reiners et al. 2005). We return to this point later.

Taken together, it should be evident that at least from the lunar record no firm conclusions can be drawn about the proposed mechanics of the early (before 4400 Ma) bombardment process, nor of the precise timing of giant planet migration that caused it. A new approach is warranted, and for this we turn to the asteroids.

The diverse members of the asteroid belt (DeMeo & Carry 2014) completed their formation before the inner planets and are thus witnesses to the solar system's earliest post-primary accretion history. We use compiled chronological data for asteroidal meteorites and compare different radiogenic systems with different T_c in different minerals and whole-rock samples to analyze this history. Our approach circumvents the debate between the monotonic decline and LHB "cataclysm" scenarios cited above because asteroidal meteorite ages are not anchored to the lunar record, and there is no dispute that they preserve a history that predates the Moon's own formation. Robust chronometers like U–Pb in zircon from meteorites such as the brecciated eucrites—a subset of vestoids that have undergone substantial reworking through impacts—document thermal reset ages (Hopkins et al. 2015) for a system with substantially higher closure temperature (T_c ∼ 1000°C; e.g., Reiners et al. 2015) than, for example, ^40–39Ar (250–450°C). Minerals that host radiogenic systems with intermediate T_c, such as Rb–Sr and U–Pb in phosphate (e.g., basic Ca-phosphate, or apatite) as well as whole-rock Sm–Nd geochronology (Borg et al. 2015), are expected to record ages of intermediate bombardment intensity. This approach bridges the gap between the highly retentive U–Pb in zircon and the low retentivity of ^40–39Ar (Figure 2). Accordingly, T_c provides a clear theoretical basis for understanding meteorite mineral and whole-rock ages as cooling ages arising from the interplay between the kinetics of diffusion (or annealing), diffusant properties, crystal dimensions, and accumulation rates in cooling radio-isotopic systems from processes such as thermal effects of impacts.

Zoom In Zoom Out Reset image size

Our thesis is simple: high-T_c highly retentive radiogenic systems will chiefly record older ages as opposed to the less retentive radiogenic systems. The most retentive systems such as U–Pb in zircon (when available) will only ever record their own formation in a melt, or (rarely) the most energetic events caused either by higher impact velocities or by impacts from larger bodies, or both. The prediction is that U–Pb in zircon registers ages of crust formation, or reset ages when the impact intensity is so high that the crust is molten and zircon crystallizes in melt sheets (see Wielicki et al. 2012). The ^40–39Ar system will only retain ages when the impact flux was low(er) and the solar system had dynamically relaxed.

An LHB-style surge in the impact spike from giant planet migration at any time in solar system history is expected to reset ages for the asteroid belt and appear as a cluster of ages in a variety of radiogenic systems, and not just in ^40–39Ar.

To understand this further and explain the data in Figure 2, we now use our compendium from meteorites in Supporting Information (Table S2) with dynamical models to better constrain when the giant planets could have commenced their migration. This task requires that we model the temporal impact bombardment history onto the terrestrial planets by accounting for the major asteroids, leftover planetesimals of terrestrial planet formation, and comets arriving from the outer solar system. The outcomes of these models for late accretion, linked as they are to the geochronological constraints imposed herein, are then used to revise previous analyses of the thermal consequences of impact bombardments to the crusts of the planetary bodies of the inner solar system. Special focus is placed on the proclivity of the early Hadean Earth to support an emergent biosphere.

The bombardment history of the terrestrial planets consists of three components: (i) dynamical evacuation of the (E-belt) asteroids, (ii) a contribution from leftover planetesimal material from terrestrial planet formation, and (iii) an influx of comets. In the following sections, we discuss the dynamical setup that determines the amount of late accretion onto the Moon and the terrestrial planets from E-belt asteroids and leftover planetesimals. The contribution from comets is discussed separately below owing to the fact that the dynamical model is much more involved.

5.1. The E-belt

The E-belt is a now-empty implied extension of the asteroid belt all the way to Mars that was evacuated by giant planet migration (Bottke et al. 2012). We compute the contribution to late accretion from the E-belt by running 32 dynamical N-body simulations of the evolution of this reservoir in the presence of the gravitational influence of the Sun and the eight planets on their current orbits. We also added asteroids 4 Vesta and 1 Ceres as massive bodies. Each simulation also has 512 massless asteroids. Our model employed the following initial conditions: the asteroids have a uniform distribution in semimajor axes a from 1.7 to 2.1 au and a main-asteroid-belt-like Gaussian distribution in eccentricity e and inclination i, with mean values of 015 and 85, respectively, and standard deviations of 007 and 7°, respectively. The other three angles (longitude of the ascending node, argument of perihelion, and mean anomaly) are chosen uniformly randomly between 0° and 360°. The differences between each simulation are the different initial conditions of the planetesimals, even though these are statistically identical. We also ensured that initially none of the asteroids were on Mars-crossing orbits.

The simulations were run for 1 Gyr of model time with SWIFT RMVS3 (Levison & Duncan 1994) with a time step of 3.65 days. Planetesimals were removed once they were farther than 100 au from the Sun (whether bound or unbound) or when they collided with a planet or with the Sun. The impact probabilities were independently obtained from counting direct impacts onto the planets and from employing a post-processing procedure (Wetherill 1967) that computes the average impact probability. The latter is the only source to compute the impact probability with the Moon. We have previously tested this computational method against counting direct impacts with the planets, and we find that the two results match those in Bottke et al. (2012). The impact probabilities and average impact velocities for the terrestrial planets, 4 Vesta and 1 Ceres, are listed in Table 2. Analysis of the dynamics of late accretion to the terrestrial planets is forthcoming in future work.

Table 2.  Impact Probabilities, Average Impact Velocities, and Mass Accreted from E-belt Material with the Terrestrial Planets, the Moon, 4 Vesta, and 1 Ceres

Planet	$\langle {p}_{\mathrm{imp}}\rangle$ (%)	$\langle {v}_{\mathrm{imp}}\rangle$ (km s−1)	Expected Mass Striking the Surface (ppm)
Mercury	${0.90}_{-0.03}^{+0.15}$	${41}_{-24}^{+27}$	${71}_{-56}^{+128}$
Venus	${5.4}_{-1.0}^{+0.7}$	${26}_{-14}^{+19}$	${28}_{-17}^{+39}$
Earth	${5.5}_{-1.0}^{+0.7}$	${21}_{-9}^{+14}$	${24}_{-15}^{+35}$
Moon	${0.25}_{-0.04}^{+0.03}$	${18}_{-12}^{+16}$	${93}_{-80}^{+154}$
Mars	${2.2}_{-0.2}^{+0.1}$	${11}_{-5}^{+11}$	${94}_{-67}^{+155}$
4 Vesta	0.016	9	1500
1 Ceres	0.0022	11	56

Download table as:  ASCIITypeset image

We caution that the average impact velocity increases with time and that the values listed here are averaged over 1000 Myr. Consequently, a more refined temporal analysis of the flux exclusive to a particular epoch (e.g., limited to only 4450–4440 Ma) requires that we run many simulations to establish convergence. We also tabulate the expected amount of mass striking the surfaces of Earth, the Moon, and Vesta from Monte Carlo impact experiments (Brasser et al. 2016) and their 2σ uncertainties (see below). The uncertainties obtained from the assumed shallow size–frequency distribution are much larger than that in the impact probabilities. We further computed that the impact probability of objects from the main asteroid belt with the Moon is at least a factor of 30 lower than that of the E-belt, consistent with earlier estimates (Morbidelli et al. 2010). Since the two reservoirs are thought to initially have a similar mass (Bottke et al. 2005), the contribution from the main asteroid belt is thus insignificant in the first 1000 Myr and henceforth neglected in our analysis.

5.2. Leftover Planetesimals

We have previously constrained the mass in leftover planetesimals from terrestrial planet formation (Brasser et al. 2016) based on the HSE abundance in the lunar mantle (Walker 2009). We then used that information to show that the purported "late veneer" (an early, albeit significant component of late accretion subsequent to core formation) to Earth and Mars consisted mostly of single, large impactors of dominantly enstatite or ordinary chondrite composition (Genda et al. 2017; Brasser et al. 2018; Woo et al. 2018, 2019). Because our analysis of leftover material requires that the Moon was present, we start our dynamical simulations at the time of the Moon-forming event, assigned here at 4500 Ma (Barboni et al. 2017). We also possess a database of terrestrial planet formation simulations in the framework of the Grand Tack model (Walsh et al. 2011). Snapshots of the planetesimal population were taken 60 Myr after the start of the simulations, which we approximately assign at 4500 Ma. We split the total number of planetesimals into 32 cases of 1074 massless planetesimals and added all the major planets, plus 4 Vesta and 1 Ceres, on their current orbits. Simulations were run with SWIFT RMVS3 for 500 Myr with a time step of 3.65 days with the same removal conditions described above for the E-belt simulations. Results for Earth, the Moon, and Vesta are provided in Table 3. Uncertainties are 2σ.

Table 3.  Average Impact Probabilities and Impact Velocities of Leftover Planetesimals with the Terrestrial Planets, the Moon, 4 Vesta, and 1 Ceres

Planet	$\langle {p}_{\mathrm{imp}}\rangle$ (%)	$\langle {v}_{\mathrm{imp}}\rangle$ (km s−1)	Expected Mass Striking the Surface (ppm)
Mercury	${5.2}_{-1.5}^{+1.4}$	${36}_{-20}^{+24}$	${919}_{-734}^{+1280}$
Venus	${17}_{-3.1}^{+3.7}$	${26}_{-12}^{+15}$	${170}_{-130}^{+177}$
Earth	${13}_{-2.0}^{+2.5}$	${22}_{-9}^{+13}$	${115}_{-89}^{+130}$
Moon	${0.59}_{-0.20}^{+0.12}$	${19}_{-12}^{+15}$	250 (fixed)
Mars	${0.85}_{-0.17}^{+0.13}$	${15}_{-9}^{+10}$	${79}_{-69}^{+139}$
4 Vesta	5.4 × 10−4	13	534
1 Ceres	7.6 × 10−4	13	710

Download table as:  ASCIITypeset image

No direct record exists of how much mass was in these leftover planetesimals. Lunar HSE abundances show that the Moon accreted a further 0.025 wt.% after its formation while it still had a magma ocean. Combined with the average impact probability with the Moon and employing a size–frequency distribution akin to that of the main asteroid belt suggests that the mass in leftover planetesimals at 4500 Ma is approximately ${4.7}_{-3.7}^{+4.4}\times {10}^{-3}$ Earth masses. To study the effect of the size–frequency distribution on the mass accreted by the planets, we rely on Monte Carlo impact experiments.

We recently combined N-body and Monte Carlo simulations to establish the nature of late accretion onto the terrestrial planets and the likely total mass accreted by Earth and Mars when calibrated to lunar HSE abundances (Brasser & Mojzsis 2017). Due to the assumed shallow size–frequency distribution of the impactors, the ratio of the total mass accreted by Earth versus that of the Moon is generally much higher than the ratio of their gravitational cross sections (Bottke et al. 2010). The accreted mass is dominated by stochastic accretion from a few large impactors (Tremaine & Dones 1993). We presupposed that the bodies in the simulations had minimum diameters of 1 km, a maximum diameter of 2000 km, and a bulk density of 2500 kg m⁻³ for all objects regardless of their size. Two cases are tested: we either presume the same shallow size–frequency distribution as that of the main asteroid belt, or with a changeover to a steeper cumulative slope of 4.5 at a projectile diameter of 50 km (Johnson et al. 2016). For further details of the method we refer the reader to Brasser et al. (2016). The expected values of material accreted by Earth and Mars are much lower than their actual values of 0.7 wt% and 0.8 wt%, respectively (Day et al. 2016); these outcomes are discussed extensively in Hartmann et al. (2000). We did not include 4 Vesta and 1 Ceres in the Monte Carlo code because they are not expected to be struck by large bodies, so the listed values are estimates obtained from the dynamical simulations.

5.3. Comets

The cometary impact flux onto the terrestrial planets is caused by the late migration of the giant planets. By late we mean well after the gas of the protosolar nebula has dissipated (Table 1); we assume that the gas disk is dispersed by ca. 4563 Ma based on the youngest chondrule ages (e.g., Morris et al. 2015; Bollard et al. 2017). Late giant planet migration is caused by the scattering of distant planetesimals beyond Neptune and by mutual scattering among the giant planets (e.g., Tsiganis et al. 2005).

We computed the cometary flux onto the terrestrial planets in three sequential steps. Step 1: Find the best initial conditions for the giant planets and the planetesimal disk that ultimately has the highest probability of reproducing the current architecture of the terrestrial planets. Step 2: Compute the fraction of comets from the trans-Neptunian disk that venture closer than 1.7 au from the Sun as a proxy for objects that reach the terrestrial planets. Step 3: Compute the impact probability of these comets with the terrestrial planets. Step 4: Apply Monte Carlo impact simulations to the comets to compute the amount of material that strikes Earth and the Moon.

These steps are now explained in detail.

Step 1: Initial conditions. The evolution of the giant planets during late migration is chaotic (Tsiganis et al. 2005), and the probability of the planets ending up near their current orbits is low: ∼5% (Brasser & Lee 2015; Nesvorný 2015a). The aim is to increase that probability by varying the initial conditions and decide which combination of three input parameters best reproduces the current configuration of the giant planets. Here we follow the procedure outlined in Wong et al. (2019).

The initial conditions employed here are identical to a loose five-planet configuration (Nesvorný 2011) because this configuration best reproduces the current orbital structure of the Kuiper Belt (Nesvorný 2015b). To start, we assume that there are five giant planets in the quadruple resonance 3:2, 3:2, 2:1, 3:2. The initial conditions for this configuration were provided by D. Nesvorný (2019, personal communication), but we allowed the initial semimajor axis of Jupiter to vary (first free parameter); the initial semimajor axes of the other planets were computed from their resonant locations. Variations in the initial conditions within each set were obtained by giving a uniformly random deviation of 10⁻⁶ au to the position vector of each planetesimal.

In our simulations, the planetesimal disk consisted of 3000 objects. The surface density of the planetesimal disk scales with heliocentric distance as Σ ∝ r⁻¹. The total disk mass is the second free parameter. The inner edge of the disk is located 1 au outward from the outermost ice giant; the outer edge of the disk is our third free parameter.

We run many sets of 128 simulations each wherein the three free parameters are permutated. We make use of the symplectic integrator SyMBA (Duncan et al. 1998). Planets and planetesimals were removed once they were farther than 1000 au from the Sun (whether bound or unbound) or when they collided with a planet or ventured closer than 0.5 au to the Sun. We evolve every simulation using a time step of 0.4 yr to 500 Myr or until fewer than four planets remain, whichever occurs first. The migration evolution and final semimajor axes and eccentricities for the giant planets are provided in Figure 3. We show cases for a planetesimal disk of 18 (red symbols) and 19 Earth masses (blue symbols). The black circles and error bars denote the current semimajor axes and eccentricities and the extent of their secular variation. The arrows indicate the direction and amount of migration. We conclude that the combination of parameters that yields the highest probability of reproducing the current solar system according to the criteria of Nesvorný (2015b) has Jupiter initially at 5.6 au, a disk mass of 18 Earth masses, and an outer edge at 27 au. These parameters are used as input for Step 2.

Zoom In Zoom Out Reset image size

Step 2. Compute the fraction of comets that enter the inner solar system.  In order to calculate the impact probability of comets onto the Moon and the terrestrial planets, we use a piecewise approach and split the work across two sets of simulations.

The model used the best initial conditions determined in Step 1 to run one set of 512 simulations with the same initial configuration of the planets, disk mass, and disk outer edge, but now comets were removed from the simulation when they ventured closer than 1.7 au from the Sun. The vectors and the time of crossing 1.7 au are stored for the next step.

We find that 22% of all comets from the trans-Neptunian disk crossed the 1.7 au barrier, or the equivalent of ca. 4 M_⊕. The outcome of this experiment is used in step 3.

Step 3. Calculate the impact probability of the comets with the terrestrial planets.  All the comets that crossed the 1.7 au barrier were then integrated for 10 Myr in simulations that included all eight planets on their current orbits, as well as 4 Vesta and 1 Ceres. This duration is about two orders of magnitude longer than the average dynamical lifetime of a Jupiter-family comet (Brasser & Morbidelli 2013). The comets were introduced into the simulation with the position and velocity when they crossed the barrier from Step 2. All runs were performed with SWIFT RMVS3 with the same time step and removal criteria as for the leftover planetesimals.

The comets are dynamically controlled by Jupiter, and consequently they only ever spend little time in the inner solar system (Levison & Duncan 1997; Di Sisto et al. 2009). We therefore do not record any physical impacts on the terrestrial planets and resort to an averaging procedure (Wetherill 1967) that computes the impact probability and therefore the total accreted mass. Table 4 shows the cometary impact probabilities with Earth, the Moon, and 4 Vesta when taking into account the 22% probability of the comets venturing into the inner solar system from Step 2 herein. We also list the average impact velocity and the expected amount of mass in comets that strikes the surface of each of these bodies.

Table 4.  Impact Probabilities and Average Impact Velocities of Cometary Material with the Terrestrial Planets, Moon, 4 Vesta, and 1 Ceres

Planet	$\langle {p}_{\mathrm{imp}}\rangle$ (×10−6)	$\langle {v}_{\mathrm{imp}}\rangle$ (km s−1)	Expected Mass Striking the Surface (ppm)
Mercury	${0.13}_{-0.08}^{+0.22}$	${34}_{-5}^{+10}$	${42}_{-13}^{+14}$
Venus	${1.9}_{-1.3}^{+2.9}$	${26}_{-2.5}^{+3.8}$	${42}_{-7}^{+7}$
Earth	${3.1}_{-2.8}^{+3.3}$	${22}_{-1.7}^{+2.9}$	${57}_{-8}^{+8}$
Moon	${0.16}_{-0.08}^{+0.20}$	${20}_{-1.9}^{+3.7}$	${239}_{-70}^{+64}$
Mars	${2.5}_{-1.2}^{+3.3}$	${13}_{-0.6}^{+1.9}$	${409}_{-42}^{+59}$
4 Vesta	0.0093	11	3900
1 Ceres	0.025	11	2900

Note. Monte Carlo simulations were only run for Earth, the Moon, and Mars.

Download table as:  ASCIITypeset image

Step 4: Monte Carlo impact experiments.  We apply the same Monte Carlo experiments to the comets to compute the amount of material that impacts Earth and the Moon. The simplest size–frequency distribution that matches the Kuiper Belt is used, with a change in slope at a diameter of ca. 60 km and cumulative slope indices of 4.8 and 1.9, respectively, at the high and low ends (Fraser & Kavelaars 2009). Since the distribution exhibits a knee at diameters of 60 km, most of the mass in the belt resides in objects of approximately this size, and thus we expect that the ratio of accreted masses between the planets is approximately equal to their gravitational cross sections. We assumed a bulk density of 1400 kg m⁻³ (∼2× computed average comet density; Peale 1989). Outcomes of these analyses are shown in Figures 4–6.

Zoom In Zoom Out Reset image size

To synopsize this section, we describe how bombardment to the terrestrial planets consists of comets, leftovers of accretion, and E-Belt asteroids. The results of our analyses for the rates of impacts from these three sources to Earth and the Moon are reported in Figure 5.

Zoom In Zoom Out Reset image size

Bombardment that accompanied late accretion strongly modulated the initial physical, chemical, and (pre-)biological state of the Hadean Earth. Here we assess how our new mass-production functions for late accretion affected the young Earth's crust using improved impact parameters for comets, leftovers of accretion, and E-belt asteroids.

Zoom In Zoom Out Reset image size

Our new global thermal cratering analysis consists of (i) a stochastic cratering model that populates a cuboid representing the Hadean Earth's surface with craters using the mass-production functions derived from the results of our new dynamical models; (ii) analytical expressions that calculate a temperature field for each crater (Abramov & Mojzsis 2009; Abramov et al. 2013); (iii) a three-dimensional thermal model of the terrestrial lithosphere, where craters are allowed to cool by conduction and radiation; and (iv) instantaneous plus cumulative melt production, as well as defined thermal volumes of the crust as a function of time (Kieffer & Simonds 1980; Pierazzo & Melosh 2000). The model output reports the amount of melting in the lithosphere, both instantaneous and cumulative, as well as thermal volumes in the top 4 km of the lithosphere (surface to depth) at any point in time corresponding to temperature fields known for various microbial communities in the crust (mesophile, thermophile, hyperthermophile).

6.1. Late Accretion to the Hadean Earth Surface

The outputs from simulations of impact bombardment experienced by Earth between 4500 and 4000 Myr are shown in Figure 7. In our baseline model, the lithosphere thickness does not change from 4500 to 4000 Ma and remains 24 km thick (solid lines in Figure 7(a)). We also ran a model where we start with a magma ocean at 4.5 Ga and have lithospheric thickness increase with time due to radiative cooling (dashed line in Figure 7(a)). Mean surface temperature and geothermal gradient were assumed to be 20°C and 70°C km⁻¹, respectively, and a basaltic composition was assigned for the initial lithosphere. Total delivered mass was 7.76 × 10²¹ kg, equivalent to 0.13% of Earth's mass (Supplementary video file 1), slightly more than the baseline model, although 0.065 wt% delivered mass was also tested. The size–frequency distribution of these late accretion impacts was approximated using the size–frequency distribution of the main asteroid belt. The velocity distribution of the impactors, as well as the decline of impact flux with time, was derived from dynamical modeling as described herein. Impactor density of these rocky bodies was assumed to be 3000 kg m⁻³, and impact angle of each impactor was stochastically generated from a Gaussian distribution centered at 45°. Some of our simulations began with a magma ocean scenario (Supplementary video file 2). We explored this scenario owing to the fact that we deem it to be a likely consequence of the largest post-Moon formation impactors in late accretion implicated, for example, in Earth's postulated late veneer (e.g., Genda et al. 2017). Based on previous work by us (Abramov & Mojzsis 2009) and others (Chyba 1990; Marchi et al. 2014), we estimate that a minimum total delivered mass of 9.31 × 10²¹ kg of accreted material is required to melt the terrestrial lithosphere at any particular point in time. Assuming asteroid 1 Ceres's density, this corresponds to an object of approximately 1200 km in diameter (about the size of Pluto's moon, Charon) striking Earth to melt the crust. This condition is also satisfied by a multitude of smaller objects, but we consider that scenario to be less likely than the single-strike case (see Brasser et al. 2016 for this argument). At a 13% impact probability with Earth, the total mass in leftover planetesimals required to account for melting the terrestrial lithosphere is 0.011 Earth masses at 4500 Ma; this is well within acceptable parameter space (Brasser et al. 2016).

Zoom In Zoom Out Reset image size

6.2. Asteroid 4 Vesta as an Archive

Our compilation of asteroidal meteorite age data does not show a common LHB-like spike pattern near 3950 Ma, and we can rule out a specific solar-system-wide bombardment at that time caused by giant planet migration. Anticipated LHB patterns do not exist at 4100–4200 Ma either, nor are they present back to ca. 4400 Ma. Highly retentive radiogenic U–Pb and Pb–Pb ages older than 4530 Ma—variably preserved in ordinary chondrites, vestoids, and other meteorite groups in Figure 2—are attributable to parent body formation ages (Ireland & Wlotzka 1992; Bogard 1995). Crustal processes inherent to these small bodies that could reset ages after 4530 Ma are precluded owing to the demise of the principle short-lived radioactive heat source to the early solar system ²⁶Al (e.g., Lichtenberg et al. 2016) and the dearth of other radioactive heat sources with sufficient power output (e.g., ⁶⁰Fe, ²⁴⁴Pu). Hence, our preferred interpretation of ages substantially younger than 4530 Ma for any asteroidal meteorite is that they document impact-induced thermal resets past a particular T_c of a radiogenic system from a protracted decline of late accretion.

Intriguingly, a number of meteorite classes display patterns of U–Pb and Pb–Pb in zircon and apatite age resets that end at about 4450 Ma. The whole-rock Pb–Pb and Sm–Nd systems, with intermediate T_c of 600°C (below that of U–Pb in zircon, but above ^40–39Ar), terminate later, at around 4400 Ma. We find that ^40–39Ar ages tend to pick up near 4480 Ma, roughly where U–Pb, Pb–Pb, and Sm–Nd ages leave off, with little overlap between the groups. Argon data also show a continuum in ages extending from about 4480 to 3000 Ma, and younger. We view this behavior of sequential termination of ages in minerals with decreasing T_c as the expected result of a long monotonic decline of leftover debris of planet formation and the product of occasional inter-asteroid collisions (Bottke et al. 2005). The result is in agreement with crater counting statistics for the Moon and generally comports with the meta-analysis of compiled lunar ^40–39Ar data by Boehnke & Harrison (2016). It is decidedly at odds, however, with any bombardment model that calls for cataclysmic spikes in impact activity at any time in solar system history since about 4450 Ma.

It may also be the case that the relative rarity of resolved ^40–39Ar ages older than ca. 4450 Ma for the asteroidal meteorites results from the "Stonewall Effect" (Hartmann et al. 2000; Hartmann 2003), in which a history of earlier thermal events becomes masked. In the Stonewall scenario, the last impacts energetic enough to reset ages are recorded, and earlier incidents tend to become so severely overprinted that they become effectively obliterated. The near absence of reset ages younger than 4450 Ma in, for example, robust U–Pb and Pb–Pb chronologies for asteroidal meteorites also agrees with trends of declining bombardment intensity in line with our new dynamical models. Collectively, meteorite age distributions viewed through the lens of dynamical studies presented herein show that giant planet migration commenced around 4480 Ma, since the ensuing cometary bombardment records no cluster of ages at any time after.

7.1. Giant Planet Migration Occurred around 4480 Ma

7.2. Life Could Have Emerged on Earth Any Time after 4400 Ma

The giant planets radially migrated after dissipation of the protoplanetary nebula around the young Sun, yet the timing of this migration is unconstrained. Here we show that radiometric ages from various meteorite classes can be reconciled with giant planet migration if it occurred before 4480 Ma. Our result resolves crater chronologies with late accretion timescales and dynamic models. It also implies that conditions on the Hadean Earth during late accretion do not preclude the emergence of an enduring biosphere as early as ∼170 million years after solar system formation.

This work was funded by the Collaborative for Research in Origins (CRiO), which is supported by The John Templeton Foundation (principal investigator: Steven Benner/FfAME); the opinions expressed in this publication are those of the authors and do not necessarily reflect the views of the John Templeton Foundation. S.J.M. and N.M.K. were supported by the NASA Cosmochemistry program (grant No. NNX14AG31G) in part of this work. S.J.M. thanks the Earth-Life Science Institute (ELSI) at the Tokyo Institute of Technology and the Geophysical Fluid Dynamics Group at the Swiss Federal Institute of Technology (ETH) in Zürich, on visiting professorial appointments during significant phases of this research. R.B. is grateful for financial support from JSPS Kakenhi (17KK0089). O.A. and S.J.M. thank the NASA Solar System Workings Program, grant No. 80NSSC17K0732. S.C.W. appreciates the support by the Research Council of Norway via the Centre of Excellence grant to the Centre of Earth Evolution and Dynamics (CEED 223272). We thank Steven Benner, Thomas Carell, Guillaume Caro, Nicolas Dauphas, Vera Fernandes, Bernard Marty, Michael Russell, and Grenville Turner for discussions. Editorial handling by Michael Endl and reviews by an anonymous reviewer helped to improve this work.

S.J.M. formulated the study, supervised all work, interpreted data, and wrote the manuscript. R.B. developed dynamical models of accretion, helped design the study, analyzed and interpreted the results from the dynamical models, and cowrote the manuscript. N.M.K. compiled the thermochronologies of solar system objects. O.A. wrote the revised analytical thermal bombardment codes and ran simulations of the crustal evolution of the Hadean Earth using our new impact bombardment production functions. S.C.W. performed supportive crater statistical analyses and shared her crater chronology data.

The authors declare no conflict of interest.