Formation of Lunar Basins from Impacts of Leftover Planetesimals

In the standard model of terrestrial planet formation (Wetherill 1990), accretional collisions between 1 and 1000 km planetesimals lead to gradual build up of lunar-to-Mars-size protoplanets that gravitationally interact and further grow during the late stage of giant impacts (Chambers & Wetherill 1998). Hafnium–tungsten (Hf–W) isotopic system analyses indicate that the Moon-forming impact on proto-Earth happened t = 30–150 Myr after the appearance of the first solar system solids T ≃ 4.57 Ga (Kleine & Walker 2017). The newborn Moon was molten, gradually cooled down, and was eventually able to support impact structures on its surface. This is time zero for the lunar crater record. The lunar surface recorded ≃50 basin-scale impacts of d > 20 km bodies (Miljković et al. 2016) since time zero, at least some of which formed relatively late. The radiometric dating of Apollo impact melts indicates that the Imbrium basin, which dominates the lunar near side, formed at t ≃ 650 Myr (Zhang et al. 2019). From all other lunar basins only Orientale and Schrödinger have lower accumulated density of superposed craters than Imbrium (Fassett et al. 2012), and must therefore be younger.

Having (at least) three basin-forming impacts happening at t ≳ 650 Myr is unexpected from the planet formation perspective because in the inner solar system, where the accretion processes have relatively short timescales (<100 Myr), the impact flux should have rapidly declined over time. The ubiquity of ∼3.9 Gyr ages in the Apollo samples has therefore motivated the impact spike hypothesis where it was assumed that the Imbrium-era impacts mark an epoch of enhanced bombardment (Tera et al. 1974), and prompted a search for possible causes. For example, it has been suggested that a dynamical instability in the outer solar system—if it happened suitably late—could have destabilized asteroid and comet reservoirs, and produced an Imbrium-era spike (Gomes et al. 2005). The late-instability model, however, has fundamental problems with the survival of the terrestrial planets (Kaib & Chambers 2016; e.g., the fully formed Earth often collides with Venus) and the delay itself (Nesvorný et al. 2018), leaving the problem of the origin of lunar basins unresolved.

To settle this matter, we construct a physical model for three key populations of impactors in the inner solar system: (i) leftover planetesimals in the terrestrial planet zone (0.3–1.75 au), (ii) main-belt asteroids (1.75–4 au), and (iii) comets. The model is based on N-body simulations of planets and small bodies over the age of the solar system (Section 2). The simulations follow the growth of terrestrial planets, collisional and dynamical decline of populations, and record impacts of (i)–(iii) on the terrestrial worlds. For consistency, they include the outer planet instability (Tsiganis et al. 2005)—a brief period in the early solar system when the giant planets were scattered to their current orbits. The instability is assumed to happen within t ∼ 10 Myr after the protoplanetary gas disk dispersal (Clement et al. 2018). We note that the nature and timing of the instability are inconsequential for the main results presented here.

2.1. Leftover Planetesimals

2.2. Asteroids

A dynamical model for asteroid impactors was published in Nesvorný et al. (2017a). The model used the same setup for the planetesimal-driven migration and instability of the outer planets as the terrestrial planet formation model described above. To start with, the terrestrial planets were placed on the low-eccentricity and low-inclination orbits. The surface density profile of asteroids was assumed to follow Σ(r) ∝ r⁻¹. The radial profile was smoothly extended from r > 2 au, where the model can be calibrated on observations of main-belt asteroids (see below), to r = 1.75–2 au. This fixed the initial number of bodies in the now largely extinct E-belt (Bottke et al. 2012). The results had large statistics (50,000 model asteroids) and full temporal coverage (4.57 Gyr).

The flux of asteroid impactors was calibrated from today's asteroid belt. We showed that the model distribution at the simulated time t = 4.57 Gyr (i.e., at the present epoch) was a good match to the orbital distribution of asteroids. The number of model d > 10 km asteroids at t = 4.57 Gyr was set to be equal to the number of d > 10 km main-belt asteroids (≃8200), as measured by the Wide-field Infrared Survey Explorer (Mainzer et al. 2019). When propagated backward in time—using the simulation results—this provided the number of asteroids and asteroid impactors over the whole solar system history.

We empirically approximated the asteroid impact flux. In terms of the cumulative number of Earth impactors with diameters >d at times >t, the best fit yields

$$\begin{eqnarray}&&F{(d,t)={F}_{1}(d)\exp [-(t/\tau )}^{\alpha }]+{F}_{2}(d)T,\end{eqnarray} \tag{ 1 }$$

with τ = 65 Myr, α = 0.6, T = 4570 − t, and t in Myr. The first term in Equation (1) accounts for the decline of asteroid impactors during early epochs. The second term represents the constant impact flux in the last 3 Gyr. There are two size-dependent factors in Equation (1). F₁(d) is assumed to follow the size distribution of main-belt asteroids; the fit to simulation results gives F₁(10 km) = 225. Also, from the main-belt size distribution, F₁(1 km) = 3.0 × 10⁴. F₂(d) is calibrated on modern near-Earth asteroids (NEAs). Nesvorný et al. (2021a) estimated ∼3 impacts of d > 10 km NEAs on the Earth per Gyr; we thus have F₂(10 km) = 3 × 10⁻³ Myr⁻¹.

2.3. Comets

A dynamical model for comets was developed in Nesvorný et al. (2017b). To start with, a million cometesimals were distributed in a disk beyond Neptune, with Neptune on an initial orbit at 23 au. The bodies were given low orbital eccentricities, low inclinations, and the surface density Σ(r) ∝ r⁻¹. The disk was truncated at 30 au to assure that Neptune stopped migrating near its current orbital radius at ≃30 au. The simulations were run from the time of the gas disk dispersal to the present epoch. The effects of outer planet (early) migration/instability, galactic tides, and perturbations from passing stars were accounted for in the model. The results were shown to be consistent with the orbital distribution of modern comets, Centaurs, and the Kuiper Belt.

The size distribution of outer disk cometesimals is calibrated from the number of large comets and Centaurs observed today, the size distributions of Jupiter Trojans and Kuiper-belt objects, and from the general condition that the initial setup leads to plausible migration/instability histories of the outer planets. The calibration gives ∼6 × 10⁹ d > 10 km and ∼5 × 10⁷ d > 100 km cometesimals in the original disk. The size distribution is expected to closely follow a power law with the cumulative index γ ≃ −2.1 for 10 < d < 100 km, and have a transition to a much steeper slope for d > 100 km. The distribution is a product of the initial size distribution that was modified by collisional grinding. We account for the physical lifetime of comets following the method described in Nesvorný et al. (2017b).

The impact flux of comets on the terrestrial worlds is computed with the Öpik algorithm (Bottke et al. 1994). The results of the Öpik code are normalized to the initial number of comets in the original disk (see above). The calibrated model gives us the flux of cometary impactors over the whole history of the Solar System. An excellent approximation of the cumulative impact flux of comets on the Earth is

$$\begin{eqnarray}\begin{array}{rcl}F(d,t) & = & {C}_{{\rm{s}}}(d)\left\{{F}_{1}\exp \left[-{\left(t/{\tau }_{1}\right)}^{{\alpha }_{1}}\right]\right.\\ & & +\,\left.{F}_{2}\exp \left[-{\left(t/{\tau }_{2}\right)}^{{\alpha }_{2}}\right]+{F}_{3}(4570-t)\right\},\end{array}\end{eqnarray} \tag{ 2 }$$

with F₁ = F₂ = 6.5 × 10³, τ₁ = 7 Myr, α₁ = 1, τ₂ = 13 Myr, α₂ = 0.44, F₃ = 4 × 10⁻³ Myr⁻¹, C_s(d) = 1 for d = 10 km, and t in Myr.

The integrated history of lunar impacts (Figure 1) shows that leftover planetesimals dominated the early impact flux (t < 1.1 Gyr or T > 3.5 Ga; T is measured looking backward from today). Asteroids took over and produced the most impacts in the last ≃3.5 Gyr. This may explain why the size distribution of modern lunar impactors is similar to NEAs, but that of early impactors was not (Strom et al. 2005; Minton et al. 2015). The cometary flux was never large enough, in the whole history of the inner solar system, to be competitive. Indeed, the isotopic signatures of comets are difficult to find in lunar samples (Joy et al. 2012). The overwhelming majority of craters observed on the lunar surface must date back to T > 3.5 Ga, when the impact flux was orders of magnitude higher than it is today. The model predicts ≃100–500 d > 20 km impacts for t > 30–150 Myr (T < 4.42–4.54 Ga), which can be compared to only ≃50 known basins (Neumann et al. 2015). We therefore see that the number of impacts suggested by the model would be excessive if the Moon formed at t = 30–150 Myr (Kleine & Walker 2017) and recorded all large impacts since its formation.

Zoom In Zoom Out Reset image size

This suggests that the lunar record is incomplete. The Moon was fully molten when it accreted from the debris disk created by the giant impact on proto-Earth (Canup et al. 2022). The subsequent evolution and solidification of the lunar magma ocean (LMO) was controlled by a number of geophysical processes, including tidal heating, formation of an insulating flotation crust and crust overturn (Meyer et al. 2010; Elkins-Tanton et al. 2011). The basins that formed when the LMO was still present would have extremely reduced topographic and crustal thickness signatures (Miljković et al. 2021); they may be unidentifiable today. Suppose, for example, that the lunar surface started recording basin-scale impacts at t ≃ 190 Myr (T ≃ 4.38 Ga)—the oldest known basins (e.g., South Pole–Aitken) would date back to this time. If so, the model implies that ≃50 d > 20 km impacts should be recorded (Figure 1), in close agreement with the number of known lunar basins. The estimated time of LMO solidification, t = 160–220 Myr or T=4.35–4.41 Ga, is consistent with constraints from the radiogenic crustal ages and highly siderophile elements (HSEs) in the lunar mantle (Elkins-Tanton et al. 2011; Borg et al. 2015; Morbidelli et al. 2018; Zhu et al. 2019). A detailed analysis of HSEs will be published elsewhere (Nesvorný et al. 2022, in preparation).

Our impact model gives general support to empirically derived lunar chronologies (Neukum et al. 2001; Marchi et al. 2009; Robbins 2014), but differs in specifics (Figure 2). Classically, the lunar N₁(T) chronology function—the number of accumulated D > 1 km craters in km² of the lunar surface since T—was obtained by fitting the measured crater densities on terrains with known radiometric ages. Neukum et al. (2001), for example, suggested ${N}_{1}(T)=a[\exp ({bT})-1]+{cT}$ with a = 5.44 × 10⁻¹⁴ km⁻², b = 6.93 Gyr⁻¹, and c = 8.38 × 10⁻⁴ Gyr⁻¹ km⁻². The exponential term is identified here with the declining impact flux of leftover planetesimals (Figure 1). A stretched exponential function $\exp [-{(t/\tau )}^{\alpha }]$, however, more accurately approximates the declining flux in our simulations. The half-life of impact decline is δ t_half = τ(t/τ)^1−α /(2α), where τ = 1/b ≃ 144 Myr and α = 1 in Neukum et al. (2001), and τ = 6 Myr and α = 0.45 in this work. This gives fixed δ t_half ≃ 72 Myr in the classical chronology, but changing half-life in our model (e.g., t_half = 46 Myr for t = 200 Myr and t_half = 88 Myr for t = 650 Myr). ⁵

Zoom In Zoom Out Reset image size

The Imbrium-era basins represent a crucial constraint on any impact chronology. In the model, ≃0.31 Imbrium-scale lunar impacts (d ≥ 100 km; Miljković et al. 2013; Schultz & Crawford 2016; Miljković et al. 2021) happen for t > 600 Myr. From the standard Poisson statistics, and folding in a generous ∼50% uncertainty in the model flux calibration (Section 2), we estimate that the Imbrium formation was a 15%–35% probability event. This is low, but not an alarmingly low probability, especially because there is a very good reason for that. There are only two basins, Orientale and Schrödinger that formed after Imbrium. They were produced by d ≃ 50–64 km and d ≃ 20 km impactors, respectively (Johnson et al. 2016a; Miljković et al. 2016), which is consistent with the model expectation of ∼2 basin-scale lunar impacts for t > 600 Myr (Figure 1). Having only two smaller basins with post-Imbrium formation ages, however, is surprising. For example, if the current asteroid-belt size distribution is adopted for reference, there should (on average) be ∼7.4 d > 20 km impacts for every d > 100 km impact. This suggests that the Imbrium basin formed unusually late, by chance, to have only two smaller and younger basins than Imbrium (there would be many more younger/smaller basins otherwise), and justifies the 15%–35% probability of Imbrium formation quoted above. In fact, our impact chronology is (nearly) optimal to satisfy the Imbrium-era constraints (Figure 3).

Zoom In Zoom Out Reset image size

The Earth receives ≃20 times more impacts than the Moon. When a large impactor strikes the Earth, it produces a vapor-rich ejecta plume. As the plume cools down, glassy spherules form and fall back, producing a global layer that can be several millimeters to many centimeters thick. Some ∼16 spherule beds have been found in the late Archean period (T = 2.5–3.5 Ga; Marchi et al. 2021), although preservation biases and incomplete sampling may be an issue. At least some of these layers may have been produced by very large, d ∼ 50 km impactors, but most are thought to record d > 10 km impacts (Marchi et al. 2010, Johnson et al. 2016b). This can be compared with the model predictions. We estimate ≃20 d > 10 km impacts on the Earth for T = 2.5–3.5 Ga (Figure 4), of which ∼2 should be d > 50 km. The leftover planetesimals and main-belt asteroids contribute equally to impacts in this time interval (∼10 impacts each). Whereas the asteroid impacts should be more uniformly spread over the late Archean, nearly all planetesimal impacts happen before ∼3 Ga. The model gives ∼10 and ∼2 d > 10 km asteroid impacts on the Earth in the past 2.5 and 0.6 Gyr, respectively, which is consistent with the current impact flux of large NEAs (Nesvorný et al. 2021a).

Zoom In Zoom Out Reset image size

Even though comets were dwarfed by planetesimals and asteroids in terms of the overall bombardment, evidence for cometary impacts can be found in the composition of Earth's atmosphere (Marty et al. 2016). Comets start to be released from the trans-Neptunian region near the onset of the outer planet instability. The cometary impact profile is found here to be more extended in time than previously thought (Morbidelli et al. 2018), with 10% of cometary impacts happening at >55 Myr and 1% at >370 Myr after the instability. Adopting a case with the outer planet instability at t < 10 Myr and Moon formation at t ∼ 50 Myr, we estimate that the Earth would have accreted ∼2 × 10²² g of cometary material after the Moon-forming impact. This is consistent with comets being the main source of noble gases in the Earth's atmosphere (Marty et al. 2016; but a negligible source of Earth's water— the mass of Earth's oceans is 1.4 × 10²⁴ g). In fact, the atmospheric abundance of noble gases can be used to constrain the delay between the outer planet instability and Moon formation, Δt. We estimate 20 < Δt < 60 Myr, with the exact value depending on the physical lifetime of comets (Section 2). This suggests, if the instability happened very early (t < 10 Myr; Clement et al. 2018; Liu et al. 2022), that the Moon must have formed early as well (t < 70 Myr or T > 4.5 Gyr; Thiemens et al. 2019).

In a broader context, our work highlights the importance of impacts for the early Earth. The bulk of accreted mass was brought by 80–400 d > 100 km (Imbrium-scale) impacts and the stochastic accretion of very large planetesimals (d ≳ 1,000 km), as needed to explain Earth's HSEs (Bottke et al. 2010). If the surface of Earth's Hadean crust was widely reprocessed by these impacts, this could explain the age distribution of Hadean zircons and absence of terrestrial rocks older than 4.3 Gyr (Marchi et al. 2014). There is roughly a 50% chance that the last d > 100 km impact on the Earth happened as late as T ≲ 3.5 Ga. In total, the early Earth received ∼2,000–10,000 d > 10 km impacts, each representing the magnitude of the Cretaceous–Paleogene (K/Pg) extinction-scale event (Alvarez et al. 1980). By T ∼ 3.7 Ga, when we have the earliest firm evidence for biotic life (Nutman et al. 2016), the impact flux has already declined ∼1000 times, and the mean interval between K/Pg-scale impacts stretched to >1 Myr.

The main results of this work can be summarized as follows.

• 1.  
  The leftover planetesimals produce most lunar impacts in the first 1.1 Gyr (t < 1.1 Gyr or T > 3.5 Ga). Asteroids produce most impacts in the last 3.5 Gyr. The transition from leftover planetesimals to asteroids has been imprinted in the crater size distributions (Strom et al. 2005, Head et al. 2010, Orgel et al. 2018). The comet contribution to the crater record is found to be insignificant (for the early instability case adopted here).
• 2.  
  Some 500 d > 20 km planetesimals from the terrestrial planet zone (0.5–1.5 au) are expected to impact the Moon since its formation. The early crater record must have been erased because the lunar surface was unable to support basin-scale impact structures. The ∼50 known lunar basins formed after t ≃160–215 Myr (T ≲ 4.36–4.41 Ga). This is consistent with the long-lived LMO (Morbidelli et al. 2018, Zhu et al. 2019). The South Pole–Aitken basin should date back to T = 4.36–4.41 Ga.
• 3.  
  About two lunar basins are expected to form for t ≳ 650 Myr (T ≲ 3.92 Ga). The Imbrium basin formation (T ≃ 3.92 Ga, d ≳ 100 km impactor) is estimated to happen with a 15%–35% probability in our model. Imbrium should have formed unusually late, relative to the expectations from the lunar impact chronology, to have only two smaller/younger basins than Imbrium (Orientale and Schrödinger); there would be many more younger/smaller basins otherwise.
• 4.  
  The lunar chronology can be given as a sum of two terms: the stretched exponential function (the leftover planetesimal branch) and a constant (the asteroid or NEA branch). This is similar to the classical (empirical) crater chronologies (Neukum et al. 2001; Marchi et al. 2009; Robbins 2014), except that the cratering rate profile in the first ∼1 Gyr had a longer tail than the exact exponential. This can lead to modest, ∼50 Myr differences in the estimates of lunar basin ages.
• 5.  
  Our model predicts ≃20 d > 10 km impacts on the Earth for T = 2.5–3.5 Ga. This is similar to the number of known spherule beds in the late Archean (Bottke et al. 2012; Johnson et al. 2016b; Marchi et al. 2021). Both the leftover planetesimals and main-belt asteroids contribute to impacts in this time interval. Whereas the asteroid impacts were more uniformly spread over the late Archean, nearly all planetesimal impacts should have happened before 3 Ga.

Bottke et al. (2007) studied lunar impacts of leftover planetesimals and concluded that the leftover planetesimals cannot produce basin-scale impacts on the Moon during the Imbrium era, because the collisional and dynamical decline of planetesimal impactors was presumably too quick. Comparing their impact profile with the ones obtained here, we identify the main reason behind this: the planetesimal population in Bottke et al. (2007) dynamically decayed by 2 orders of magnitude in the first 100 Myr. We do not find any such strong initial trend in our work. The problem in question is most likely related to the approximate nature of initial conditions in Bottke et al. (2007), where the present-day NEAs were used as a proxy for terrestrial planetesimals. We note that modern NEAs have relatively short dynamical lifetimes.

Brasser et al. (2020) developed a dynamical model for lunar impactors with four different components: E-belt, asteroid belt, comets, and leftover planetesimals. The main difference relative to our work is that their model was not absolutely calibrated from independent means. To match the lunar impact constraints in their model, Brasser et al. (2020) increased ∼10 times the population of the E-belt (see Bottke et al. 2012), and suggested that the E-belt asteroids were the main source of basin-forming impacts. The best-fit contribution of planetesimals was found to be negligible; leftovers were therefore concluded to be only a minor source of lunar impacts. This can be compared to our work where we find the dominant role of leftover planetesimals.

The work of D.N. was supported by the NASA Emerging Worlds program. F.R. acknowledges support from the Brazilian National Council of Research—CNPq. The work of D.V. was supported by the Czech Science Foundation (grant No. 21–11058S). A.M. received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant No. 101019380 HolyEarth). R.D. acknowledges support from the NASA Emerging Worlds program, grant 80NSSC21K0387. The simulations were performed on NASA's Pleiades Supercomputer. We greatly appreciate the support of the NASA Advanced Supercomputing Division. The SWIFT code is available on www.boulder.swri.edu/~hal/swift.html. The data that support the plots within this Letter and other findings of this study are available from the corresponding author upon reasonable request.