© ESO, 2017

1. Introduction

When the Nice model for the long-term dynamical evolution of the planetary system was first proposed (Tsiganis et al. 2005), an important virtue was found to be the possibility to explain a late lunar impact cataclysm, commonly referred to as the late heavy bombardment (LHB; Gomes et al. 2005). Since then, the dynamics that constitutes the core of the model has undergone several changes (Morbidelli et al. 2010; Levison et al. 2011), and the nature and timing of the LHB has also been further discussed (Morbidelli et al. 2012; Bottke et al. 2012).

Even so, some fundamental features have survived. The giant planets started out in a much tighter configuration than they exhibit at present. Beyond their orbits there was a planetesimal disk with a total mass of order 20−50 Earth masses. As the instability of the planetary orbits set in, this disk was dispersed across the whole solar system following the outward migration of Uranus and Neptune. An inevitable consequence is the cratering of planetary surfaces in the inner solar system (including that of the Moon).

However, the dynamical instability and rapid migration of the giant planets would also have affected other small body reservoirs closer to these target objects. One is the asteroid main belt (Gomes et al. 2005; Morbidelli et al. 2010), and another is the so-called E-belt (Bottke et al. 2012) related to the present Hungaria asteroid population. The question hence arises, whether the lunar and planetary cratering was dominated by trans-planetary (in fact, cometary) or asteroidal projectiles.

Gråe Jørgensen et al. (2009) argued for cometary bombardment based on Ir abundances in the Isua Greenstone belt (Greenland) rocks and lunar Apollo samples, while Strom et al. (2005) found the size distribution of lunar craters to be consistent with asteroidal projectiles and different from the likely cometary crater size distributions on the giant planet satellites. Moreover, Bottke et al. (2012) estimated that practically all the mare basins could have an E-belt origin, leaving nearly nothing for the comets.

This issue is not yet resolved. In addition, there is an open question about the amount of volatile delivery to Mars during the LHB (see Lammer et al. 2013). In the current paper we aim to address these issues by a quantitative estimate of the time-integrated lunar and planetary impact rates for a projectile population of the cometary type. The thrust of this paper is our new method to derive the impact probabilities by Monte Carlo simulation from a massive set of long-term orbit integrations for the projectiles in question.

In an accompanying paper (Rickman et al. 2016, hereinafter “Paper I”) we describe in detail how we construct our basic sample of orbital evolutions for the cometary projectiles in terms of initial orbits, dynamical model and integrations. A brief summary is provided in Sect. 2, and Sect. 3 describes the principles of using the MOID (minimum orbit intersection distance) to identify cases of potential collisions in the orbital evolution sample. In Sect. 4 we introduce the concept of physical evolution and lifetimes for the cometary projectiles and describe how we model this. In Sect. 5 we discuss the modelling of a secondary projectile population that arises from the deflections of projectile orbits at close encounters with the terrestrial planets. Section 6 provides the details of how the lunar and planetary impact probabilities are found using MOID methods (Rickman et al. 2014). In particular, a special technique for the lunar case is introduced here. In Sect. 7 we derive the actual impact rates for the Moon and planets, and we investigate the effect of different assumptions for the total mass and size distribution of the trans-planetary disk as well as different physical evolution models. Finally, a discussion of our results and a summary of our conclusions are provided in Sect. 8.

2. Orbital evolutions

We used initial orbits emanating from the calculations by Brož et al. (2013). These represent objects reaching for the first time into the inner solar system with perihelion distances less than 3.9 au. After down-selection and cloning, our sample consisted of 100 000 initial orbits. This down-selection means that we disregard a small part of the initial projectile population, mainly with P > 20 yr. We will derive impact rates per million objects, hence referring specifically to initial orbits with P < 20 yr. When these are combined with absolute numbers of objects transferred from the external disk, the disregarded objects are assumed to contribute the same specific impact rate as the retained objects. Our results are not severely biassed by this approximation.

The integrations were performed with the classic Radau integrator RA15 (Everhart 1985), using as dynamical model the gravitational six-body problem of the Sun, the four giant planets and the massless projectile. We use the current orbits of the giant planets, meaning that we do not attempt to treat the planetary instability or migration, which in the Nice Model is the reason for the delivery of the cometary objects into the orbits that we integrate. All these events occurred on a very short time scale (Morbidelli et al. 2010), and we assume that the LHB cometary bombardment took place soon afterwards, when the planets had settled into their new orbits.

In order to get good estimates of the impact probabilities associated with all orbits of the test objects, we had to use the best suited heliocentric osculating elements for these orbits. We deemed the best osculation epoch to be the time of perihelion and thus went into some effort to find this with a good accuracy. Our chosen epochs are typically within one week of the actual perihelion times, and these deviations do not matter for our purposes.

Each test object performed a pre-defined maximum of 100 000 revolutions around the Sun, but the integrations were ended prematurely in case the test object reached a semi-major axis exceeding 1000 au, reached an eccentricity exceeding 0.9999, or collided with the Sun or a planet. As shown in Paper I, this resulted in a median lifetime of ~6400 orbits, and the end state of ejection or orbital diffusion was the dominant one, comprising 93% of the cases. About 2.6% underwent collisions, and 4.4% survived after 100 000 orbits.

3. Identifying potentially collisional orbits

In Fig. 5 of Paper I we presented a graphical illustration of the minimum perihelion distances reached by the 100 000 test objects, and from the same data we have extracted the numbers of objects that reached perihelion distances smaller than the aphelion distances of the terrestrial planets. We present these numbers in Table 1. They provide a first indication of how the relative impact probabilities for the different planets are conditioned by the general nature of the dynamical evolutions that we have found.

We now aim to use the orbital integrations to estimate the number of collisions of the simulated cometary objects with the Moon and the terrestrial planets. The basic tool of our method is the calculation of MOIDs (see Wiśniowski & Rickman 2013). In particular, we associate each orbital revolution in our integration output with a set of target orbits to be described in Sect. 3.1, and we calculate the MOID for each target. If this is small enough, we deem that there is a chance for a collision during the orbital revolution in question. Our way to quantitatively determine this probability is described in Sect. 6.

3.1. Lunar and planetary orbits

Since we are dealing with a time period about 4 Gyr ago, we should assume a geocentric orbit for the Moon that is consistent with the current understanding of the long-term evolution of the lunar orbit. Of course, the tidal interactions in the Earth−Moon system have caused a gradual expansion of this orbit, but there is a large uncertainty over the integrated amount of this acceleration over the previous 4 Gyr. We base our estimate on results by Bills & Ray (1999) and use a semi-major axis for the Moon that is one half of the present value, that is, about 182 000 km. For simplicity, we take this orbit to be circular and situated in the ecliptic plane. Further, we treat the Earth−Moon system with the approximation that the Moon is effectively massless, and the system’s centre of mass resides in the Earth’s centre, moving on its elliptical orbit around the Sun.

We use a statistical approach to represent the orbits of the terrestrial planets. The semi-major axis is assumed to be the average one characterizing the present motion. The eccentricity and inclination are taken to be constant during the whole time interval covered by our integrations. This is often long enough to cover the secular periods of the planets, and therefore we use rough estimates of the averages plotted by Laskar (1994) for epochs about 4 Gyr ago. We list the values of the heliocentric a, e and i for all four terrestrial planets, and the geocentric values for the Moon, in Table 2.

It is true that the real motions of the terrestrial planets are characterized by quasi-periodic oscillations of e and i coupled to the circulations of ω and Ω, respectively. But since we neglect these oscillations, we are free to model the evolution of the angles at will. We have chosen to treat them as stochastic variables with values picked at random on the interval [0,2π] for each projectile orbital revolution under consideration.

3.2. MOID-based screening

Using all the elements described, we compute the MOID (ℳ) of the projectile and each planetary target for each orbit of the former, using the method of Wiśniowski & Rickman (2013). While the angular elements of the planets are taken at random, those of the projectiles are given by our integrations. Thus, the calculated MOID values reflect the varying collisional probabilities imposed by the real trends in the cometary orbital evolutions as governed by the perturbing, giant planets. This is important, because it is well known (e.g. Di Sisto et al. 2009) that the ω distribution of Jupiter family comets in particular is far from uniform, and it would hence be inaccurate to use the average impact probabilities for uniform angular elements presented by Rickman et al. (2014) or computed by the Öpik (1951) or Wetherill (1967) formulae.

We compare each ℳ with the collisional radius (R_coll) of the planet. This is done using a specific value of R_coll, calculated from the standard gravitational focussing formula with the actual approach speed of the projectile under the relevant geometric circumstances. If ℳ <R_coll, the collision probability is non-zero. We list the physical radii and escape speeds that we use for the calculation of R_coll in Table 3.

In Table 4 we list the numbers of ℳ <R_coll orbits – called collisional orbits – found for all the projectile orbital evolutions. The trend among the different planets is partly due to their different sizes, and partly due to the fact that planet-crossing orbits result from temporary visits into small perihelion distances. As shown in Table 1, the number of such visits is significantly smaller for Mercury-crossing than for Mars-crossing conditions.

The case of lunar impacts is special. Here we introduce a special limit ℳ₀ = 184 000 km for the terrestrial MOID, because when ℳ_E > ℳ₀, the lunar impact probability is always nil, while for ℳ_E < ℳ₀ it is often non-zero. We hence select the orbits with ℳ_E < ℳ₀ for calculating the respective lunar impact probabilities.

4. Physical evolution

We now consider the effects of physical evolution affecting the sizes of the projectiles and causing their eventual destruction. For this purpose, and for lack of a priori information about the LHB projectiles, we shall assume that they behaved like the current Jupiter family comets (JFCs). We consider three different options. In our reference scenario, there is no physical evolution at all, so all effects are simply ignored. In the other two scenarios, we assume a size distribution for the initial nuclei, and we let the radii of the individual nuclei decrease, orbit by orbit, due to the mass loss that the comets suffer. In one case, we consider the erosion caused by near-surface ice sublimation and associated dust outflow. In the other case, we investigate the consequences of the evolutionary scenario presented by Di Sisto et al. (2009), using the case that they found to yield the best fit to the orbital distribution of JFCs.

Specifically, in the two scenarios of physically evolving nuclei, we initiate a sample of 100 000 comets with given sizes, and we ascribe different orbital evolutions to all of them. For each subsequent orbit, we calculate the total amount of previous erosion in terms of the decrease of the radius, and we hence keep track of the current radius. When this is smaller than 0.1 km, we consider the comet to be destroyed. The initial size distribution is accounted for by weighting the objects according to their initial radii. Thus, our test objects may be considered to represent a real population with some given total number and the chosen distribution of initial radii.

4.1. Erosion model

In Paper I we present our model of surface erosion, which is based on calculations of the H₂O free sublimation flux as a function of heliocentric distance using a simple thermal model. This flux is scaled by a factor 0.05 to account for a reduced activity of JFCs, estimated by Di Sisto et al. (2009) to be 15% on the average, and the apparent existence of twice as many dormant as active JFCs (Paper I), indicating that JFCs experience no mass loss for 2/3 of the time. These results are applied to all the orbits of our sample by integrating the flux with respect to true anomaly. The integrated H₂O mass loss per orbit is augmented by an equal amount of mass assumed to represent the dust outflow. Finally, the depth of the eroded layer is found by using an assumed JFC density to be discussed in Sect. 7.1.

In this model we do not account for mass loss by splitting or fragmentation of the comet nuclei. We admit that this is a questionable assumption, but in our opinion the importance of splitting for the destruction of comet nuclei in general is not proven beyond doubt. However, in the next considered model splitting is the major contributor to the loss of comets.

4.2. Di Sisto model

In this model, the comet nucleus loses mass due to surface sublimation of H₂O ice, and thermal modelling of this process has yielded a polynomial fit to the amount of mass loss per unit area per orbital revolution as a function of the perihelion distance (Eq. (3) of Di Sisto et al. 2009). Dividing this by the bulk density of the nucleus, one gets an erosion depth, which represents the idealized case of a freely sublimating surface. This is multiplied by a factor less than unity, representing the “active fraction” of the comet. Following Di Sisto et al. (2009), we have used a constant active fraction of 0.15. Unfortunately, they did not mention which bulk density they used, so we take 500 kg m^-3 as a representative value from recent estimates of JFC densities.

We note for completeness that we modified their equation for the mass loss slightly. The formula is intended for use only up to perihelion distance q = 2.5 au, and we naturally adhere to this. However, close to the upper limit of the range of q, the polynomial function has a minimum value close to zero, after which the curve turns upward again. We have removed this increase by putting the mass loss equal to zero beyond the minimum. It is also worth emphasizing that our own erosion model differs in several respects from the one of the Di Sisto model. First, the sublimation flux that we use is generally larger than that of Di Sisto et al. Second, our active fraction is smaller, as explained above. Thirdly, in our model we consider sublimation up to a maximum of 6 au heliocentric distance, and erosion is hence modelled for all perihelion distances below this limit.

In Di Sisto et al. (2009), splitting events are modelled as stochastic with a probability f of happening per orbital revolution, given as a function of the perihelion distance by $\mathit{f}\mathrm{=}{\mathit{f}}_{\mathrm{0}}\mathrm{(}\mathit{q}\mathit{/}{\mathit{q}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{-}\mathit{\beta }}\mathit{.}$(1)Each splitting event involves the loss of a fraction s of the current mass of the nucleus, where s is given by $\mathit{s}\mathrm{=}\frac{{\mathit{s}}_{\mathrm{0}}}{\mathit{R}\mathit{/}{\mathit{R}}_{\mathrm{0}}}$(2)as a function of the radius R of the nucleus. In summary, the frequency of splitting events decreases with increasing perihelion distance, and the relative mass loss per event decreases with increasing size of the nucleus.

As parameters for this splitting model, we took those of Di Sisto et al.’s best fit case (their Model 4):

• β = 0.5

• f₀ = 1

• q₀ = 0.5 au

• s₀ = 0.001

• R₀ = 10 km.

This means that comets are assumed to split about once per orbit, when the perihelion distance is less than 1 au, losing about 1% of the mass each time, when the radius is about 1 km. We modified the Di Sisto et al. approach, as far as they described it, in order to model correctly the situation when f> 1, since this is of significant occurrence. Thus, we used a Poisson distribution with expectance f to draw a random integer, and we let the nucleus split this number of times during the orbit in question. After adding all the contributions to the mass loss for each orbit, we calculated the total erosion depth and applied this as a decrease of the nucleus radius.

The Di Sisto et al. (2009) investigation differed somewhat from ours regarding the dynamical model. Their initial conditions were different, the integrator was different, and they included explicitly Venus, Earth and Mars in addition to the giant planets as perturbing bodies. In addition, in contrast to us, they added a slight non-gravitational force to the equations of motion. Nonetheless, we have applied their physical evolution model without corrections to our sample of orbital evolutions.

However, the use of a physical model obviously requires the assumption of a size distribution for the comet nuclei. For this, Di Sisto et al. (2009) used a current JFC cumulative size distribution (CSD) based on Tancredi et al. (2006) and Fernàndez & Morbidelli (2006), assumed to be valid for q< 2.5 au: $\begin{array}{c}\\ {\mathit{N}}_{\mathrm{1}}\mathrm{(}\mathit{>}\mathit{R}\mathrm{)}\mathrm{=}\mathit{C}{\mathit{R}}^{\mathrm{-}{\mathit{s}}_{\mathrm{1}}}\hspace{0.17em}& \mathrm{0.1}\mathit{<}\mathit{R}\mathrm{\le }\mathrm{1}\\ {\mathit{N}}_{\mathrm{2}}\mathrm{(}\mathit{>}\mathit{R}\mathrm{)}\mathrm{=}\mathit{C}{\mathit{R}}^{\mathrm{-}{\mathit{s}}_{\mathrm{2}}}\hspace{0.17em}& \mathrm{1}\mathit{<}\mathit{R}\mathrm{\le }\mathrm{10}\end{array}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}$(3)with s₁ = 1.3, s₂ = 2.7, and R is given in km. They arbitrarily took a few values for the initial radius and applied these to equal numbers of initial comets. These evolved independently using the physico-dynamical model and as a result, a sample of comet orbits with q < 2.5 au was obtained, where each orbit was associated with a specific, current value of the nucleus radius. If the radius of a nucleus decreased below 0.1 km, the comet was considered lost. Weights were applied to different radius intervals in order to scale the relative numbers of orbits in the simulated sample into the observed size distribution as described above. From the weighted numbers, distributions of orbital elements could be constructed for comparison with observations. The preferred model of physical evolution was chosen based on this comparison.

Our problem is different. We do not need to test different physical evolution models, since we can use the one that Di Sisto et al. (2009) found to be the best. On the other hand, while these authors did not opt to solve for the best initial size distribution in order to get as close as possible to the observed size distribution without weighting, we deem this to be an essential check, and we will present our results in the following section.

4.3. Influence on size distribution

In both our models for the physical evolution of comet nuclei, we expect that this evolution may affect the slope of the comet size distribution, so that the power law index of the current CSD differs from that of the source population, from which the comets are captured. We have thus performed special simulations to measure the extent of this effect. The method is to launch a set of 100 000 comets with randomly chosen initial radii, uniformly distributed between 1 and 10 km, associate them in a random way to the 100 000 orbital evolutions, and run each physical evolution model on the whole set. We may simulate any initial CSD by means of weight factors coupled to the comets. In both models we consider comets to be alive only until the radius decreases below 1 km. From the entire pool of orbits of such comets with perihelion distances less than 2.5 au including their initial weights, we can assemble the current CSD and measure its slope index.

We have thus simulated test samples with many different values of the initial CSD index for both physical evolution models, and in Table 5 we present the comparison with the resulting current indices. For the erosion model, we note that the current distribution is always slightly shallower than the initial one. In the range of slopes commonly found for observed JFCs, the difference amounts to about 0.2−0.3 in agreement with the results by Weissman & Lowry (2003). This effect may have a rather trivial explanation. In the erosion model, the depth typically eroded away is independent of the size of the comet. Therefore, if the ratio between two initial radii is x_i > 1, the ratio between the current radii of the same comets is x_c > x_i, but the ratio between the numbers of comets remains the same.

The result for the Di Sisto model is more remarkable. While an initially very shallow CSD steepens up, the effect is the opposite for an initially very steep CSD. In particular, the current CSD of the Di Sisto model apparently cannot be steeper than a value around −2.3 for the CSD index. Physically, this is reasonable, since the smallest comets suffer an extremely high death rate owing to their efficient mass loss by splitting. Hence, if we feed a very large number of small comets into the JFC population by means of a steep initial CSD, this has no effect on the current CSD, since these comets die too fast.

A consequence of this is that the Di Sisto model is not able to reproduce the steep CSD of observed JFCs, which was assumed during the construction of the model (this index was −2.7, as mentioned above). Due to this inconsistency, we prefer to use the Di Sisto model only to illustrate the effects of assuming a furious splitting rate. We do not consider it a realistic model for the actual evolution of JFCs.

In view of the above, we have opted to use the same two CSD indices (−1.5 and −2.5) as in Paper I. They yield current CSD indices that are more or less consistent with observations of JFCs (Morbidelli & Rickman 2015), independent of which physical evolution model we employ.

5. Decoupling from Jupiter

As regards dynamical evolution, the cometary LHB projectiles have a present-day analogue in the Jupiter Family comets (JFCs). Thus, the time-integrated impact rates per object that we calculate should be similar to those characterizing the JFCs. However, one feature of the JFC dynamics is not properly included into our six-body gravitational model. This is the evolution into orbits decoupled from encounters with Jupiter due to small aphelion distances. Such objects have very long dynamical lifetimes and may be inner planet crossers for a long time, increasing their chance to hit the targets.

Comet 2P/Encke is a case in point. Its current orbit has perihelion and aphelion distances q = 0.33 au and Q = 4.11 au, respectively. Its association with the Taurid meteoroid complex and possibly, the Tunguska impactor (Kresák 1978) shows its relevance for discussions of the current Earth impact risk.

Concerning mechanisms for the transfer of comets into Encke-like orbits, Levison et al. (2006) gave strong preference to the perturbations by terrestrial planets including orbital deflections experienced at close encounters. In Paper I we described our investigation of such transfer, where we used an approximation that involved a strong focus on the closest possible encounters, and we found by comparison with Levison et al. (2006) that this method is severely underefficient in producing Encke-like comets. As a result, some insight into the dynamics of transfer was gained. However, it follows that we are unable to use the deflected comets of our investigation for realistically evaluating their contribution to the impact rates.

Instead, we make a crude estimate based on the results of Levison et al. (2006). They estimated three active, Encke-like comets¹. This was based on an estimated total of 540 usual JFCs with the same sizes (nuclear diameter D > 1 km) and q < 2.5 au. According to Table 3 of Paper I, more than 60 of these would hence have q < 1 au, and about 14 would have q < 0.5 au. Since Fig. 1b of Levison et al. (2006) shows that most Encke-like comets would have perihelion distances between 0.5 and 1 au, we conclude that the steady state number of Encke-like comets is ~5% of the number of potential impactors in the usual Jupiter Family.

The contribution of Encke-like comets to the impact rates is larger, since their orbital periods are shorter and they cross the planet orbits more frequently. However, the effect of this cannot be much larger than a factor two, so we conservatively estimate that the Encke-like comets may yield impact rates ~5–10% of those of usual JFCs. In the following, we will estimate the latter, and there will be considerable error bars due to uncertainties over the mass of the primordial disk and the transfer efficiency from the disk into the inner solar system. Therefore, it seems relevant to account for the Encke-like comets by simply adding 10% to the estimated impact rates.

We note that the production of Encke-like comets according to Levison et al. (2006) necessitates a rethinking about the physical evolution of comet nuclei. Since the time scales needed are much longer than the estimated active lifetimes of JFCs, these authors speculated that comets may spend most of their time in a dormant state, being reactivated upon large reductions of the perihelion distance (Rickman et al. 1991). Such a picture deviates fundamentally from the ones we apply in this paper (Sect. 4). So far, it remains speculative, but if it proves realistic, the confidence in our reference model is considerably increased.

6. Impact probabilities

6.1. The planetary case

Whenever ℳ < R_coll so that the impact probability is non-zero, we estimate this probability by the MOID method as described by Rickman et al. (2014). In almost all cases, the variant called MOID-chord is accurate enough (see Fig. 1). It means that we consider the circular contour of the planet’s collisional cross-section as imaged on the b-plane (Greenberg et al. 1988). The locus of all crossing points of the projectile for the different timings of the encounter that yield a collision is a chord within this circle, and the full timing range is found, dividing the length of the chord by the projection of the planet’s orbital velocity onto the b-plane. Dividing this timing range by the orbital period of the planet yields the collision probability.

Fig. 1 Illustration of the MOID-chord method of calculating the random probability of collision with a target planet, given a known value of the minimum orbit intersection distance (MOID). The target’s collisional radius defines a circle in the b-plane, and the time it takes for the target’s projection to move the distance EE′ divided by its orbital period gives the probability. From Rickman et al. (2014) with permission.

Only in special cases, when the projectile orbit has i < 0.1° or 0.999 < q/q_pl < 1, do we resort to the accurate MOID-track method (Rickman et al. 2014). These cases are of extremely rare occurrence in our dynamical sample.

Our method does not involve any super-sizing of the planets. Collisional probabilities are only calculated using the actual collisional radii. It is the use of a very rich sample of encounters, resulting from our cloning and integrations, that guarantees a good statistical accuracy in our derivation of a global, time-integrated impact rate for the population of projectiles. Another virtue of our technique is that each non-zero probability refers to an impact with a given velocity, and hence we can also find the distribution of impact velocities for each planet as well as the distribution of crater sizes – see Sect. 7.4.

Figure 2 illustrates some results for 4000 randomly chosen collisional orbits (1000 per target planet) in the reference model with infinite physical lifetimes. The encounters with ℳ <R_coll are ordered chronologically for each planet, and the respective collision probabilities are marked by symbols that identify the planets. The scatter seen among the individual collision probabilities for any planet is caused by two phenomena. On the one hand, the actual MOID values differ, and on the other hand, the encounter velocity vectors make different angles with the b-plane. The trend between the different planets is explained by the fact that, as seen from the discussion in Rickman et al. (2014), the plotted probability is proportional to the ratio between the planet’s radius and its orbital semi-major axis.

6.2. The lunar case

In order to compute the probability of an impact on the Moon for a projectile with ℳ_E < ℳ₀, we have to deal with a double timing problem. The first step is to consider the geocentric torus spanned by the Moon during its orbital motion as the target, and to find the range of the Earth’s mean anomaly at the epoch, which allows the projectile to pass through this torus. If there is no such range for the projectile in question, the lunar impact probability is obviously zero. In the complementary case of a non-vanishing range of relevant Earth timing, we have to identify the part of the lunar torus that the projectile crosses as a result of the combined heliocentric orbital motions of the Earth and the projectile. Then we consider all values of the Moon’s geocentric mean anomaly at the epoch and find the range, for which a lunar impact occurs. Finally, the lunar impact probability is found, in principle, as a double integral over the two mean anomalies.

Fig. 2 Individual collision probabilities for random samples of 1000 comet orbits per planet, all of which satisfy the condition ℳ < Rcoll for the planet in question. For each planet, the orbits are ordered chronologically. The coloured symbols denote: Mercury (black circles), Venus (red squares), Earth (green triangles), and Mars (blue diamonds).

We have developed a numerical code that performs these calculations, and we show some test results in Fig. 3. Here we compare the average impact probabilities of the Earth and Moon for objects that penetrate within the distance ℳ₀ from the Earth’s centre. These probabilities have been calculated for samples of 50 000 such objects with special combinations of q and i, always using a = 3.5 au. For these checks we took the Earth’s orbit as circular at 1 au in the ecliptic plane. We plot the ratio between the impact probabilities of the Earth and the Moon as a function of the perihelion distance of the projectiles for a constant inclination of 10°.

Fig. 3 The ratio (S) between the impact probabilities for the Earth and the Moon as found in test calculations (see the main text) is plotted versus the perihelion distance of the projectile for a constant inclination of 10°. Green dots: gravitational focussing neglected. Red and blue dots: gravitational focussing included, using the Moon’s ancient and contemporary orbital semi-major axis, respectively.

The figure illustrates several results. First, we note that replacing the Earth by the Moon without gravitational focussing would reduce the cross-section by a factor 13.45. The ratio S plotted in the diagram for the case when gravitational focussing is neglected (green dots) is always larger than 13.45, showing that in this case the fact of the Moon’s geocentric orbital motion has the effect of enhancing the Earth/Moon ratio. This is likely a geometrical effect rather than a consequence of the Moon’s velocity.

We also recognize that including the gravitational focussing has a large influence on the S ratio, when the projectiles have perihelion distances in the vicinity of 1 au, while the effect is only moderate for smaller perihelion distances. Of course, gravitational focussing occurs mainly for the Earth because of its much larger escape velocity. This is enhanced by the fact that the unperturbed approach velocity to the Earth−Moon system decreases, as q approaches 1 au. The red dots show this trend for the ancient lunar orbit. In fact, calculating S including gravitational focussing for the case where the Moon is at rest at the place of the Earth yields a value of 20.59 for q = 0.9 au, indicating that the Moon’s orbital motion has a similar effect with gravitational focussing as without it. We have also checked the influence of the inclination of the projectile orbit and found it to be very small for the typical values encountered. The explanation is that the approach velocity is given by the projectile’s Tisserand parameter relative to the Earth. This depends explicitly on cosi, which varies little over the range in question.

We have repeated the calculation of S(q) with gravitational focussing, using the contemporary semi-major axis of the lunar orbit, and the results are shown by the blue dots. The good agreement between S(q) for the two orbits shows that our choice of radius for the ancient orbit is not critical for the derived impact rates.

For completeness, we mention a peculiar behaviour that is not shown in Fig. 3 but appears for q values extremely close to unity. We then find that S(q) takes very large values (several hundred) due to the fact that the Earth’s impact probability exhibits a very high ridge at q = 1 au, when plotted versus q and i (Rickman et al. 2014), while the geometric effect of the Moon’s orbit is to wash out this ridge. However, the projectile orbits for which this effect would be appreciable are extremely few in our sample.

7. Impacts on the Moon and planets

7.1. Disk properties

The number of cometary impacts that we find will obviously depend on how many comets we feed into the system. This in turn comes from our assumptions about the primordial, trans-planetary disk – in particular, its total mass and size distribution. In addition, we need to specify the likelihood for disk objects to be transferred into the inner solar system (q< 3.9 au) during the dispersal caused by the outward migration of Uranus and Neptune.

Following Rickman et al. (2015), we will assume a special size distribution in the shape of a broken power law for diameters less than 1000 km supplemented by 1000 Pluto-sized objects. The latter play an important role in the Nice Model (Levison et al. 2011) by stirring the planetesimal orbits to significant eccentricities. The break in the power law occurs at a radius of R_o = 50 km. For the larger objects with radii between R_o and R_max = 500 km, we use a CSD index of −4, based on observational estimates for large TNOs. For the smaller objects with radii from R_min = 0.5 km to R_o, as mentioned above, we use two values (−1.5 and −2.5) bracketing the range of most observational estimates for JFCs.

For the total disk mass M_d, a likely range of 20−50 M_E (Earth masses) is often stated, and we will use two values probing the lower and upper ends of this interval. Rickman et al. (2015) considered a minimal disk with 18 M_E contained between R_min and R_max, and 2 M_E bound up in the 1000 Plutos, in order to minimize the collisional destruction of the smaller objects. This will be one of our two alternatives, and for the other one we double the mass for radii below R_max to 36 M_E. The Pluto-sized objects are irrelevant for our impact rates, since they are far too few to hit any of our targets, and we shall neglect them in the following.

We model the differential size distribution of the disk objects in terms of a power law frequency function ϕ(R′) given by $\mathit{\varphi }\mathrm{\left(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{R}}^{\mathrm{\prime }}\mathrm{=}\mathit{C}{\mathit{R}}^{\mathrm{\prime }\mathrm{-}\mathit{s}}\hspace{0.17em}\mathrm{d}{\mathit{R}}^{\mathrm{\prime }}\mathit{,}$(4)using the dimensionless variable R′ = R/R_o. In accordance with the above values for the CSD index, we put s₂ = 5 for R′ ≥ 1 and s₁ = 2.5 or s₁ = 3.5 for R′< 1. By introducing $\begin{array}{c}\\ {\mathrm{\Phi }}_{\mathrm{1}}\mathrm{\left(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\frac{\mathit{C}}{{\mathit{s}}_{\mathrm{1}}\mathrm{-}\mathrm{1}}\mathrm{(}{\mathit{R}}^{\mathrm{\prime }\mathrm{1}\mathrm{-}{\mathit{s}}_{\mathrm{1}}}\mathrm{-}\mathrm{1}\mathrm{)}\\ {\mathrm{\Phi }}_{\mathrm{2}}\mathrm{\left(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\frac{\mathit{C}}{\mathrm{4}}\mathrm{(}{\mathit{R}}^{\mathrm{\prime }\mathrm{-}\mathrm{4}}\mathrm{-}{\mathit{R}}_{\max }^{\mathrm{\prime }\mathrm{-}\mathrm{4}}\mathrm{)}\end{array}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}$(5)we can write the CSD in terms of the total number of objects as $\mathit{N}\mathrm{(}\mathit{>}{\mathit{R}}^{\mathrm{\prime }}\mathrm{)}\mathrm{=}\{\begin{array}{c}\\ {\mathrm{\Phi }}_{\mathrm{1}}\mathrm{\left(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{+}{\mathrm{\Phi }}_{\mathrm{2}}\mathrm{\left(}\mathrm{1}\mathrm{\right)}& {\mathit{R}}_{\min }^{\mathrm{\prime }}\mathrm{\le }{\mathit{R}}^{\mathrm{\prime }}\mathit{<}\mathrm{1}\\ {\mathrm{\Phi }}_{\mathrm{2}}\mathrm{\left(}{\mathit{R}}^{\mathrm{\prime }}\mathrm{\right)} & \mathrm{1}\mathrm{\le }{\mathit{R}}^{\mathrm{\prime }}\mathrm{\le }{\mathit{R}}_{\max }^{\mathrm{\prime }}\mathit{.}\end{array}$(6)We have: ${\mathit{R}}_{\min }^{\mathrm{\prime }}\mathrm{=}\mathrm{0.01}$ and ${\mathit{R}}_{\max }^{\mathrm{\prime }}\mathrm{=}\mathrm{10}$. With these values, the parameter C is very close to four times the number of disk objects in the larger category (with R′> 1).

There are four solutions for C, corresponding to the two values of s₁ and the low and high masses of the disk. When solving for C, we assume that all objects with R ≥ 50 km have a density of 1000 kg/m³, and for all the smaller objects we use half of this (500 kg/m³). We present the resulting numbers of disk objects with R ≥ 50 km and with all sizes in Table 6.

Our assumed mass distribution for the projectile population is top heavy: assuming the flatter of the two size distributions in the lower range, the mass ratio between the larger and smaller object populations is M₂/M₁ ≈ 2.7. For the steeper variant, we have M₂/M₁ = 1.0. However, due to the steep size distribution of the larger objects, relatively little mass is actually found in the very largest ones.

In addition to the above properties, we also need to assume an efficiency of transfer from the disk into orbits with q< 3.9 au, from which we pick the initial conditions for our integrations. Here we will consider two values. Quoting previous work by Levison & Duncan (1997) and Fernández et al. (2004), Volk & Malhotra (2008) used an efficiency of 30% for scattered disk objects to reach the Jupiter Family after encountering Neptune, and this can be translated into a crude value of f = 1/3 for our efficiency. However, we may also base our estimate on the results by Nesvorný et al. (2013), who simulated the capture of Jupiter Trojans following the migration of the ice giants. They found the likelihood for a primordial disk object to hit the Earth to be (2−4) × 10^-7 (Nesvorný, priv. comm.), while we will report a probability of about 3 × 10^-6 for objects of our initial population to hit the Earth. Both results are found by assuming infinite physical lifetimes. Hence, we may estimate an efficiency of f = 1/10 by comparing the two collision probabilities.

Both values of f are crude estimates. Using more accurate values is not warranted owing to the uncertainties over the actual disk properties. We note that a value intermediate between the two comes from comparing with the (Brož et al. 2013) results. The above-mentioned Earth impact probability of 3 × 10^-6 refers to a sample of 5000 objects coming from a total of 27 000 disk objects of Brož et al. (2013; see Paper I). This directly translates into a probability of about 5 × 10^-7 for the entire disk. However, about 1000 additional disk objects reached q< 3.9 au like our initial population though with longer orbital periods. Assuming these to have Earth impact probabilities similar to those we have calculated, we can roughly increase the above estimate to 6 × 10^-7, indicating that the Nesvorný et al. disk is more excited that the Brož et al. disk. Thus, the f value of the latter would be 1/4.5.

The cratering rates that we find depend on only one parameter, namely, the product fM_d. To effectively bracket the whole uncertainty range, we opt to combine on the one hand the larger values, yielding a total mass of 12 M_E for the objects populating our starting orbits, and on the other hand the smaller values, from which the corresponding mass is 1.8 M_E. For both of these choices, we use the two values of the size distribution index s₁.

7.2. Cumulative impact rates

As described above, we use three models for the physical evolution of the cometary projectiles. In the reference model there is no evolution at all, but in the other two models (the erosional and Di Sisto models) the radii of the projectiles decrease monotonously in manners that depend on the changing properties of the projectile orbits. However, we use the same procedure to simulate the impact rates in all three cases.

In this procedure, we consider 100 000 sample comets. Each one gets its own, individual orbital evolution picked at random from our integrations. As initial nucleus radii (R_in) we consider 1000 discrete values, equally spaced in lgR_in. We divide the interval from lgR_min to lgR_max into 1000 intervals of equal width and place 100 comets at the middle of each interval. In the reference model each such set of 100 comets keep identical radii, even though the orbits evolve differently. In the other two models, since the orbital evolutions differ, the 100 radii will soon diverge.

Using the number of disk objects as a function of radius and the transfer efficiency into our starting orbits as described in Sect. 7.1, we associate each sample comet with its specific number of projectiles. We use this throughout the simulation as a weight factor. The weight factor only depends on the initial radius and thus remains constant throughout the simulation.

Concentrating on one target at a time, we look for the collisional orbits performed by the sample comets. Many comets have no such orbits at all and will not contribute anything to the impact rate, but others may have dozens. In the physically evolving models, for each new orbit by a sample comet, its radius is updated by subtracting the proper amount of erosion. The comet stays in existence until its current radius decreases below 0.1 km, when it falls victim to the erosional death and ceases to exist. As long as this does not happen, each orbit is checked for a potential impact. If the orbit is found to be collisional, the necessary data about the impact is entered into an impact data base, where we list the identity of the target, the probability of the event, the weight factor, current radius and approach speed of the projectile.

First, we derive a cumulative impact rate for each target by simply adding the collision probabilities of all orbits for that target. This refers to 100 000 projectiles entering into the inner solar system, and thus we derive the cumulative impact rate per such initial projectile for each target. These rates are presented in Table 7 for the reference model as expected numbers of impacts for a total of one million initial projectiles of all sizes. We also show the expectances for the numbers of impacts by the subset of these with radii R> 50 km.

A few comments to this Table are in order. The impact rate listed for the Earth is the one referred to above in Sect. 7.1, when a value for the transfer efficiency from the primordial disk was derived. Comparing this to the lunar impact rate, we find the average Earth/Moon ratio to be 24.9. Such a comparison was presented in more detail in Sect. 6.2 as a function of orbital elements. In this sense the value just found is an average over the elements of the Earth−Moon collisional orbits. We further note that the Mercury/Mars ratio of impact rates is as large as 0.35. Thus, even though many more comets cross the orbit of Mars than that of Mercury, and Mars presents a considerably larger collisional cross section in absolute units, the chances of impact on Mercury are enhanced by the rapid orbital motion of the innermost planet, leading to a comparable average impact rate. This effect was illustrated in Fig. 2.

Since the lifetimes of the projectiles are size independent in the reference model, the total impact rates integrated over all sizes do not depend on the initial size distribution. However, when we turn to the models of shrinking nuclei, this independence no longer holds true. Thus, in Tables 8 and 9 we present the same quantities as in Table 7, but there are now two columns for the total impact rates corresponding to the two values of the slope index of smaller projectiles.

Table 8 shows several interesting results, when compared to Table 7. The total impact rate by projectiles of all sizes decreases by about an order of magnitude for Mercury when we take the erosion by H₂O sublimation into account – a little less for the shallower size distribution and a little more for the steeper one. But the targets moving further away from the Sun are considerably less affected. In the extreme case of Mars, the total impact rate decreases by a factor 3.4−4.1 depending on the slope index. Hence, the Mercury/Mars ratio has dropped from 0.35 to the range 0.09−0.14. While Mercury was hit more than twice as frequently as the Moon in the reference model, the lunar impact rate is equal to or larger than that of Mercury in the erosional model.

The reason for the decrease of the total impact rate is that many comets with small initial radii do not survive until they would enter their collisional orbits, when erosion is allowed for. We note that this effect is large, even though our erosional model appears quite conservative in its prediction of cometary mass loss. When it comes to impacts by large objects (R> 50 km), the situation is different. Most of the comets that start in this range do not have time to shrink below the 50 km limit prior to their collisional orbits, so the decrease in this part of the impact rate is much smaller than for the total impact rate. As a crude estimate, the relative decrease amounts to about 30% as an average over all the targets.

Interestingly, the trend among the different targets is not the same as for the total impact rates. This time, Mercury has the smallest decrease (26.5%), and the Earth and Moon have the largest (31.5%). The effect is likely caused by a difference in time scales. The small comets that erode away entirely do so faster than the large comets pass the 50 km limit. For the latter, it is conceivable that secular perturbations play an important role in decreasing the perihelion distances (see Paper I), so that Mercury and Venus are protected by the slowness of this dynamical evolution.

We now see a slight change of the Earth/Moon ratio for the total impact rate. This ratio is now 26.4 – a significant increase from the reference model. The most likely explanation is that the collisional orbits now correspond to a younger population that is less dynamically excited. The approach velocities to the Earth−Moon system are thus smaller, and hence the Earth’s gravitational focusing increases in importance. On the other hand, the Earth/Moon ratio of impact rates by R> 50 km objects stays the same as in the reference model. This is natural, since the dynamical heating precedes most of these impacts just like in the reference model.

Table 9 shows the consequences of applying a model that involves a very high splitting rate, perhaps in conflict with observational evidence. The possible virtue of this model is to imply very short physical lifetimes, especially for small comets with low perihelion distances. Let us first compare the total impact rates with those of the erosional model. We see an extreme decrease for all targets but especially for Mercury. This amounts to a factor of more than 30 for the shallower size distribution, and more than 100 for the steeper one. These factors become smaller in a monotonous way, when we go to more distant targets. For Mars, they are 11 and 16, respectively.

Comparing with the reference model, we see that less than one comet in 300 manages to hit Mercury before it disappears by the physical evolution of the Di Sisto model. For Mars, we get less than about one comet in 40. Of course, this dramatic destruction mainly concerns the small comets. The comets that start out with R > 50 km suffer much smaller losses. However, these losses are larger in the Di Sisto model than in the erosional model. Now the trend is monotonous – for Mars, 76.4% of the large initial comets pass the 50 km limit before the collisional orbits are reached, while for Mercury this increases to 86.3%.

This means that the time scale for the shrinkage of the nuclei is shorter with the Di Sisto model for comets of all sizes. It seems likely that this inhibits the dynamical protection of Mercury and Venus that operated in the erosional model. Concerning the Earth/Moon ratio for the total impact rate, this has now increased to 27.2, since the impactor population is now extremely young. Interestingly, for the R > 50 km objects the ratio has now dropped to 23.7 compared to 24.9 in the reference model. A possible explanation is that, as compared with the reference model, the early impacts by relatively small objects are now suppressed by shrinkage below the radius limit, while the larger, long-lived objects tend to approach the Earth only after the orbits become dynamically excited.

7.3. The largest cometary impactors

Combining the numbers in Tables 6 and 7, we see that the expected number of impacts by projectiles with R>R_o (diameters exceeding 100 km) is not very large even in the reference model. For instance, for the Moon this number is 3.15 for the disk with the largest number of large objects and the highest transfer efficiency into the inner solar system. Hence, for this size range, which contains a large fraction of the total mass of projectiles, we need to account for the probability distribution of the actual number of impacts and derive the relevant statistics of some quantities, including the largest impactor radius and the total mass of the large impactors.

To solve this problem, we take the Poisson distribution as relevant. This means that the occurrence of an impact by a large enough projectile is independent of how many such impacts have already occurred, which seems trivially guaranteed. Thus, if λ is the expectance for the number of large impacts, we model the actual number k of large impacts that occurred in the real solar system by means of a probability mass function $\mathit{f}\mathrm{\left(}\mathit{k}\mathrm{;}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\frac{{\mathit{\lambda }}^{\mathit{k}}{\mathit{e}}^{\mathrm{-}\mathit{\lambda }}}{\mathit{k}\mathrm{!}}\mathrm{·}$(7)

Now, consider the above-described impact data base for an arbitrary target, an arbitrary disk-transfer model and an arbitrary model of physical evolution. Each entry corresponds to one collisional orbit, and we use the current radius R of the projectile, its weight w, and the impact probability p. Since w is the number of initial comets represented by this sample object, the product wp means an expected number of actual impacts by projectiles with the current radius R.

Let us order the entries according to decreasing current radii. Going through this ordered list and tracing the increasing sum of wp, we eventually come to an entry, for which this sum exceeds the value 5 (this choice is arbitrary but seems natural). At this point, we call the current radius R₅. Thus, the expectance for the number of impacts by projectiles with radius R₅ or larger is λ ≈ 5. We then perform a Monte Carlo simulation drawing 1000 random values of k according to Eq. (7). As an example, for the Moon we find R₅ to range between 19 and 42 km for the four disk models, when the reference model is applied. In the extreme case of the Di Sisto model, these values are 4 and 22 km.

For each case we generate k random values of R, using the already determined set of values for the sum of wp at different radii, by picking a uniform random number x on the interval [0,λ]. The radius for which x is closest to the sum of wp is taken as a random value of R. We take the largest of these k random R values as one random realization of R_max, the radius of the largest projectile hitting the target. We finally use the set of 1000 values of R_max to determine the median and the lower and upper quartiles, assuming these to characterize the true parent distribution of R_max.

We also calculate the masses of all these simulated impactors using the densities mentioned above. For each of the 1000 realizations of the set of impactors with R ≥ R₅, we thus get a total mass M₅. We derive the median and quartiles of M₅ and assume these to characterize its true parent distribution. For all the smaller radii we just form the sum of all the current projectile masses, each mass being multiplied by wp, to get the total impactor mass M_small. This approximation is warranted, since the number of such projectiles is large enough for its statistical uncertainty to be neglected.

The median radius of the largest impactor is invariably found to be largest for the disk model with the largest value of fM_d and the more shallow size distribution of the smaller objects. Similarly, it is always smallest for the opposite combination of the smallest value of fM_d and the steepest size distribution of the smaller objects. In the following, we focus on these two extreme models, which we call the maximum and minimum models, respectively. In Fig. 4 we present the medians and quartiles of R_max for Mercury, the Moon and Mars in the reference, erosional and Di Sisto models, using the maximum and minimum disk models.

Fig. 4 Maximum cometary impactor radius for Mercury, the Moon and Mars. Each error bar extends from the lower to the upper quartile of the cumulative probability distribution with the median marked by a filled circle for the maximum model and a cross for the minimum model of the disk. The red colour is used for the reference model, the blue colour for the erosional model, and the magenta colour for the Di Sisto model. The shaded band for the Moon indicates an observational constraint based on the Jupiter Trojan population (see the main text).

It is evident that the results for the reference and erosional models are quite similar, whereas the Di Sisto model predicts much smaller impactors. Naturally, these trends are similar for all three targets. The difference between the maximum and minimum disk models is as large as the one between the physical evolution models. In this regard, it is interesting to use an observational constraint based on the population of Jupiter Trojans. The cumulative number of Trojans brighter than absolute magnitude H is well known for the larger objects and can thus be used. Nesvorný et al. (2013) found the trapping probability of Trojans during the time of planet migration to be p_t ~ 6−8 × 10^-7 per object in the primordial disk, while the collision probability with the Moon can be estimated at p_c ~ 1−2 × 10^-8 by combining the value listed in Table 7 with the estimates of transfer efficiency f in Sect. 7.1.

Multiplying the cumulative H distribution of Trojans by p_c/p_t, the number 1 is reached at H ≃ 9.1−9.8, which for an albedo of 0.07 means a radius range of 27.5−37.5 km. We consider this to be a realistic constraint on the size of the largest lunar, cometary LHB impactor, and the range is marked in Fig. 4 as a grey band in the strip of the Moon. Evidently, this gives preference for the minimum disk model over the maximum one, and we conclude that reality is likely closer to the former.

By adding the two above-mentioned mass contributions M₅ and M_small, we get our results for the total mass M_tot of cometary LHB projectiles hitting each target. In Fig. 5 we show these results for the same targets as above. The statistical error bars refer only to the contributions by the largest objects and are thus somewhat underestimated.

Fig. 5 Total mass of all the cometary impactors on Mercury, the Moon and Mars, plotted on a logarithmic scale. The error bars, and the colours and symbols, have the same meaning as in Fig. 4.

For Mars as well as the Earth, it is interesting to compare M_tot to their estimated water inventories. In the case of the Earth, we get a maximum total impactor mass that is less than 2 × 10²⁰ kg. This is far less than the ocean mass (1.35 × 10²¹ kg), which in turn is considerably less than the Earth’s total estimated water content. Thus, clearly, most of the water was picked up by the Earth during its early era, before the LHB, as predicted by Morbidelli et al. (2000). For Mars, Fig. 5 shows that the highest obtainable values of M_tot are less than 1 × 10²⁰ kg. If we take the minimum disk model to estimate M_tot ≃ 7 × 10¹⁸ kg and use a mass fraction of 20% for water in a comet nucleus, guided by Rosetta results for comet 67P (Rotundi et al. 2015), we obtain about 10 m Global Equivalent Layer, which is much less than the estimated total inventory of martian water (Lammer et al. 2013). We conclude that the presence of liquid water on the surface of Mars during the Noachian did not occur because the planet swept up its water at that time through the LHB. This water must have arrived much earlier – probably simultaneous with or even preceding the arrival of water onto the Earth (see Brasser 2013).

In the case of the Moon, the first thing to note is that M_tot is less than 10¹⁹ kg for all model combinations. This can be compared with the mass estimates for the lunar basin-forming impactors (Levison et al. 2001), from which a crude, total LHB impactor mass can be derived. Using a median approach speed of 18 km s^-1 for the lunar projectiles (see Fig. 6 in Sect. 7.4), we estimate this mass to be ≃10¹⁹ kg within a factor of a few. Thus, even if we had confidence in the maximum disk model and evolutionary models with relatively long lifetimes, we would just barely be able to match the observations in this sense. For the minimum model, there is no doubt that comets provided only a minor part of the total mass of lunar LHB impactors.

7.4. The largest impact basins

In order to calculate crater sizes, we need to know the speed of impact, V_imp. This is obtained from the asymptotic approach velocity U and the escape velocity V_e at the surface of the target by ${\mathit{V}}_{\mathrm{imp}}^{\mathrm{2}}\mathrm{=}{\mathit{U}}^{\mathrm{2}}\mathrm{+}{\mathit{V}}_{\mathrm{e}}^{\mathrm{2}}\mathit{.}$(8)Since we know U for each collisional orbit as well as the corresponding projectile weight and impact probability, we can easily derive the distributions of U for impacts on all the targets. We present these graphically for the reference model in the left panel of Fig. 6 by means of cumulative probability curves. A significant contribution by velocities larger than the circular speed (U_c) at the mean distance of the target from the Sun is a striking feature for all planets. For Mars and the Earth, about 10% of the encounters occur with U>U_c, for Venus almost 15% and for Mercury about 30%. In these cases, the Tisserand parameters T_t with respect to the target planets fall in the range T_t < 2.

In part, the far tails of the U distributions are due to the secular heating of the projectile population, which occurs in the reference model. We showed in Paper I that this leads to a significant fraction of high inclinations. The fact that Mercury is more affected than the other planets, followed by Venus, may partly reflect a difference in dynamical age. It takes longer time, on average, for comets to be able to encounter the innermost planets than the outermost ones.

Fig. 6 Left panel: cumulative distributions of approach velocities to different planets for collisional orbits in the reference model. The curve for Mercury shows some wiggles owing to relatively sparse material. The thin vertical lines indicate the circular speeds at the mean distance of each planet. Right panel: cumulative distributions of approach velocities to the Earth for collisional orbits in the reference, erosional and Di Sisto models. For the latter two, the slope index s1 = 2.5 has been used. The wiggles in the curve for the Di Sisto model are due to the sparse material.

In the right panel of Fig. 6 we show the changes of the U distribution for the Earth, as we move from the reference model to those of shrinking nuclei. For the latter we have chosen the variant with the shallower size distribution of the smaller objects. Choosing the steeper size distribution, the results are similar, but the difference from the reference model is even more pronounced. The obvious trend is for the velocities to decrease, when the erosional effects are stronger and the lifetimes of the nuclei are reduced. The long-term dynamical heating of the Jupiter Family is thereby avoided.

Since the size of the craters (or basins) formed by impacts on the targets increases with velocity as well as radius of the projectile, this should mean that, in a statistical sense, the largest basin diameter on any target is quite sensitive to the efficiency of erosion. For a higher erosion rate and shorter lifetimes, not only does the maximum projectile size decrease, but the likelihood of a high velocity also decreases.

To calculate the probability distribution of the largest basin diameter, we use the same procedure as for the maximum projectile radius. For each entry in the impact data base we calculate the crater diameter (D) as described in the following paragraphs, using the speed of impact as calculated by Eq. (8). We then order the entries according to decreasing D, and we identify the entry with crater diameter D₅, for which the sum of wp reaches a value λ ≈ 5. The Monte Carlo simulation is made in the same way as for the maximum projectile radius, and we now determine the median and quartiles of D_max.

We use the general formula proposed by Holsapple & Housen (2007) for the radius ℛ_tr of the transient crater formed by a point source impactor with radius R, hitting the target surface with a velocity whose vertical component is v: ${\mathrm{ℛ}}_{\mathrm{tr}}\mathrm{=}{\mathit{K}}_{\mathrm{1}}\mathit{R}[\frac{\mathit{gR}}{{\mathit{v}}^{\mathrm{2}}}{\left(\frac{{\mathit{\rho }}_{\mathrm{t}}}{\mathit{\rho }}\right)}^{\mathrm{2}\mathit{\nu }\mathit{/}\mathit{\mu }}\mathrm{+}{\left(\frac{{\mathit{Y}}_{\mathrm{t}}}{{\mathit{\rho }}_{\mathrm{t}}{\mathit{v}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathrm{+}\mathit{\mu }\mathit{/}\mathrm{2}}{\left(\frac{{\mathit{\rho }}_{\mathrm{t}}}{\mathit{\rho }}\right)}^{\mathit{\nu }\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{2}\mathit{/}\mathit{\mu }\mathrm{)}}{]}^{\mathrm{-}\mathit{\mu }\mathit{/}\mathrm{(}\mathrm{2}\mathrm{+}\mathit{\mu }\mathrm{)}}\mathit{,}$(9)where ρ_t and Y_t are the density and strength of the target material, g is the surface gravity of the target, and ρ is the density of the projectile. We use ρ_t = 2700 kg/m³ and Y_t = 6 × 10⁷ dyne/cm² for all the targets, the latter being an average of the range 2 × 10⁷−2 × 10⁸ dyne/cm² (Asphaug et al. 1996). For the numerical parameters we use μ = 0.55, ν = 0.4 and K₁ = 0.93 as appropriate for rocky targets (Holsapple & Housen 2007). We assume the impacts to occur at random on the spherical target surface, yielding an average of 2/3 for the projection factor transforming V_imp into v.

The transient crater evolves by slumping into the simple crater, whose diameter is ${\mathrm{𝒟}}_{\mathrm{spl}}\mathrm{=}\mathrm{2}{\mathrm{ℛ}}_{\mathrm{tr}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{c}}_{\mathrm{sl}}\mathrm{)}\mathit{,}$(10)where the slumping coefficient c_sl is taken as 0.3 in accordance with many previous papers. The diameter of the final, complex crater is in case , and ${\mathrm{𝒟}}_{\mathrm{cpx}}\mathrm{=}{\mathrm{𝒟}}^{\mathrm{\ast }}\mathrm{(}{\mathrm{𝒟}}_{\mathrm{spl}}\mathit{/}{\mathrm{𝒟}}^{\mathrm{\ast }}{\mathrm{)}}^{\mathrm{1.18}}$(11)in case . For the transition diameter we use  km for the Moon (Melosh 1989), and for Mars and Mercury we use  km and 10 km, respectively (Kring 2006). The resulting is taken as the crater diameter D.

Fig. 7 Maximum cometary impact basin diameter for Mercury, the Moon and Mars. The meaning of the “error bars”, markers and colours is the same as in Fig. 4. The grey lines indicate the diameters of the largest basins considered to be LHB-related on the three targets. The commonly used names of these basins are written above the respective lines.

In Fig. 7 we present the results for the maximum basin diameters using the same two disk models as for the maximum projectile radii and total impactor masses. If we had confidence in the maximum model, we would say that both the Caloris and Imbrium basins could reasonably be formed by LHB comets, unless the Di Sisto model is the most reasonable picture of physical evolution. However, the above calibration using the Trojans gives stronger support for the minimum disk model. Hence, we conclude that the majority of large LHB-related impact basins on the Moon, Mercury and Mars were formed by the asteroidal component. For the Moon, the largest cometary basin would rather be similar to Mare Crisium (D ≃ 555 km) in size. This corroborates the conclusion by Bottke et al. (2012) and supports their hypothesis that the dispersal of the asteroid E-belt was the main reason for the observed impact evidence.

There is a difference between Mars and the two other targets in that none of our cometary models is even close to explaining structures of sizes comparable to the Hellas basin. In fact, our models predict basins on Mercury of at least the same size as on Mars in stark contrast to the difference between the diameters of Hellas and Caloris. A speculative suggestion that would need to be confirmed is that the E-belt source of impactors is more likely to favour Mars as target due to its closeness to the martian orbit.

7.5. The areal crater density

We have also used the crater diameters resulting from the impact data bases to form the sum of wp for all entries with D> 20 km. This represents the expectance for the total number of such craters formed by the LHB comets on the target in question. Assuming these to cover the surface uniformly, we divide by the total surface area of the target and calculate the areal density of D > 20 km craters per million km². Of course, this density is model dependent, and the resulting values span a wide range. However, guided by the preference for the minimum disk model, and noting that the model with the lower disk mass and the shallow size distribution yields similar results for the largest impactor radius, we now concentrate on these two models only.

Fig. 8 Areal density of D > 20 km craters for Mercury, the Moon and Mars. For each target the three physical evolution models are plotted using the same colours as in Fig. 4. The upper cross is for the steeper size distribution of the small disk objects, and the lower cross is for the shallower one. The low mass disk is used in both cases. The lines at the common level of 100 craters per million km2 mark the observational reference data.

Figure 8 shows the results of this calculation for the same three targets as in the previous illustrations. However, the results are now plotted on a log scale. For each physical evolution model, we mark the crater density for each of the two disk models. We also include an observational background, obtained as follows. For the Nectarian strata of the lunar highlands, there is an estimate of 100 D > 20 km craters per million km² by Marchi et al. (2012). For the highlands of Mars and Mercury, similar numbers are reported by Werner (2014). We thus use this crude number as a common crater density for all three targets – something that may seem like a remarkable coincidence.

If we had confidence in the reference model, we would conclude that the size distribution is quite far from the steep variant, because otherwise too many craters would be produced – especially on Mercury. However, we deem the reference model to be unlikely in view of the observational evidence for mass loss and destruction of comets. Hence, we should rather use the data from the evolving models, and these have one feature is common: the comets provide only a minority of the craters observed on all targets. The cometary contribution is minute in the Di Sisto model, while in the erosional model it may be significant though far from dominant.

8. Discussion and conclusions

We have simulated the impact rates on the Moon and terrestrial planets due to objects originating in the primordial, transplanetary disk, assuming that this disk was dispersed by the outward migration of the ice giants as predicted by the Nice Model. Our calculations thus refer to the cometary component of the LHB in case the timing of the planetary migration is relevant for offering an explanation of the LHB (Gomes et al. 2005). In case the migration occurred at a much earlier time (Kaib & Chambers 2016; Weaver et al. 2016), our results would refer to a much earlier cometary bombardment. This would be unrelated to the LHB, and the LHB would need to be explained in a different way. We note that the dispersal of the hypothetical E-belt (Bottke et al. 2012) would then also be unrelated to the LHB.

Thus, finding an alternative explanation for the LHB would remain as a major challenge. It is far from obvious that any small body reservoir in the solar system would then be rich enough to provide the required flux of projectiles. The scattered disk would already be severely diluted (Nesvorný et al. 2013), and the possible role of the Oort Cloud is questionable. Moreover, the geochemical evidence for an important asteroidal component of the LHB (Kring & Cohen 2002) would essentially leave the main belt as a source, but this would already have been diluted more or less to its present state, which appears insufficient to sustain the LHB.

We hence argue that our simulations most probably deal with the actual Late Heavy Bombardment. As initial conditions we use a set of 5000 orbits emanating from dynamical simulations within the Nice Model by Brož et al. (2013). These represent objects originating in the primordial disk and reaching for the first time perihelia within 3.9 au. Compared to the disk models generally used in earlier papers, the present one is relatively excited due to viscous stirring (Levison et al. 2011). As we have mentioned, this may lead to a decreased efficiency of transfer from the disk into the inner solar system from the earlier standard value of 1/3 into our new estimate of 1/10, which we find by comparing our result for the Earth impact rate to an analogous number by Nesvorný et al. (2013). Since the present scattered disk is a remnant of the primordial disk, we argue that this reduced transfer efficiency will also hold for the capture of Jupiter Family comets in contrast to most earlier estimates. This would tend to reduce the discrepancy between the estimated populations of the scattered disk and the Oort Cloud, which was found by Brasser & Morbidelli (2013).

The dynamical model in which we trace the orbital evolutions of our comet clones is a self-consistent six-body model with the Sun, the four giant planets and the massless object. We thus neglect all other influences, but in Sect. 5 we discussed the statistical chances for our comets to become deflected into Encke-type orbits by interactions with terrestrial planets. Based on (Levison et al. 2006), we estimated that the steady-state number of Encke-like comets may be ~5% of the number of Jupiter Family comets. Since the Encke-like comets have larger impact probabilities due to their shorter orbital periods, we also estimated that they would cause an additional impact rate amounting to about 5−10% on top of the one caused by the modelled comets. Since this is just a small increase in comparison with other uncertainties, we did not include it into our calculations.

We explore the possibilities for the LHB cometary bombardment as far as we can by considering on the one hand four different models of the initial population of projectiles. These are based on a high-mass and a low-mass primordial disk, associated with a high and a low transfer efficiency, respectively, and a high and a low index for the size distribution of the smaller disk objects. On the other hand, we use three different models for the physical evolution of the comets. This way we believe that we cover the whole range of possible outcomes. The physical evolutions include two extreme cases. One is to neglect any physical destruction of the comets and assume that the nuclei survive intact for an arbitrary time, and the other is based on the paper by Di Sisto et al. (2009), which argued for a model with a very large splitting rate of Jupiter Family comets that makes the lifetimes of small comets very short. Intermediate to these, we consider an erosional model accounting for mass loss by ice sublimation by using thermal model results for each particular orbit in our simulations and assuming a crude value for the activity level of the average comet nuclei.

Realistically, the first model – our reference model – must be discarded as inconsistent with the observed mass losses of comets. However, we are not able to judge, how far it is from the truth. Both the other models must be deemed primitive as well in view of the observed behaviour of comets, which often exhibit rejuvenation episodes following major orbital perturbations (Rickman et al. 1991). Moreover, for the erosional model we may have overestimated the average activity level of active comets, possibly because of an observational bias against measurements on low-active comets. Finally, the splitting rate predicted by the Di Sisto et al. (2009) model appears difficult to reconcile with observational evidence. Therefore, we prefer to keep an open mind and to use the reference model for an upper limit to the cometary impact rate with no prejudice about how large an overestimate it represents.

We now consider the statistical accuracy of our results. The best guarantee is the sheer size of our sample of orbital evolutions. We start with 100 000 objects, and in Paper I we show that 50 000 of these remain after ≃6400 orbits and 25 000 after ≃16 000 orbits. From Paper I it is also obvious that the typical ages of the projectiles are low enough to be well sampled by our dynamical material – especially when physical evolution is included. After 100 000 orbits we still have more than 4400 objects surviving, which shows that long-term dynamical effects are also reasonably well sampled. One issue concerns the identification of collisional orbits, which uses one random choice of the target angular elements per projectile orbital revolution. Could the resulting impact rates have been significantly different for a different, random choice? We have checked this by performing one extra choice, and the result is that the numbers of collisional orbits and impact rates are very similar, even though the individual encounters are of course different.

The thrust of this paper is the Monte Carlo simulation approach that we use to get full information about the collisional histories of the comets in a statistical sense. This means that for each comet orbit we calculate MOIDs with respect to all the terrestrial planets. Collisional orbits are identified by MOID values less than the planet’s collisional radius, and random impact probabilities are calculated individually for these. For the Moon we use a terrestrial MOID large enough to encompass the lunar orbit (assumed to be half its current size), and the random impact probability is computed by solving a two-dimensional timing problem involving both the Earth and the Moon. Armed with these results, we follow the dynamical and physical evolutions of 100 000 projectiles entering into the inner solar system for each of the above-mentioned model combinations, and we register all the non-zero impact probabilities along with other relevant data, as long as the projectile survives against erosion.

We let the 100 000 sample comets represent different numbers of real comets according to their initial sizes in order to simulate the size distribution and impactor population of the relevant model. These numbers are used as weight factors along with the data on impact probabilities, so that we obtain the actual, expected numbers of impacts for the different targets together with the sizes and velocities of the respective projectiles. From this we use Poisson statistics to estimate the probability distribution for the radius of the biggest projectile hitting each target. In particular, for the Moon the result is typically ~30–40 km or ~20 km depending on physical lifetimes, using the minimal model of the projectile population. If we instead consider the maximal model, the corresponding estimates are ~60–90 km or ~40–60 km. From a comparison with the capture efficiency of Jupiter Trojans by Nesvorný et al. (2013), we expect the actual value to be 27.5−37.5 km, and we conclude that the minimal model is to be preferred.

This means a primordial, trans-planetary disk of relatively low mass, that is, about 20 Earth masses or less. We note that, from Table 6, our minimal model of this disk contains roughly 1−2 × 10¹¹ objects with diameters D> 2 km. According to Brasser & Morbidelli (2013) this would imply a current scattered disk with about 1−2 × 10⁹ such members. In Paper I we came to a similar conclusion based on the steady-state capture rate of JFCs that we estimated using the same dynamical simulations as in this paper. There is hence consistency between our two independent estimates of the same quantity, which lends some credence to the result.

By a similar statistical analysis we also estimate the largest impact basin diameter for those targets, where the scars would have survived. Trusting the minimal model of the projectile population, we conclude that the largest lunar, cometary impact basin is only about 500 km in diameter – comparable to Mare Crisium but far from the size of Mare Imbrium. This means that the LHB was mostly asteroidal as far as the Moon is concerned. Similarly, for Mercury, the largest cometary impact basin is much smaller than Caloris. For Mars, the Hellas basin is much larger than any cometary basin even for the maximal model. Furthermore, we have calculated the areal density of cometary craters with D> 20 km on all these three targets for the minimal model. We conclude that for infinite physical lifetimes, this might be close to or even exceed the observed densities depending on the size distribution of the projectiles. However, for lifetimes that seem more reasonable to expect, the comets make only a minor contribution to the small highland craters on all the targets.

Our work hence supports the conclusion by Bottke et al. (2012) of an asteroidal predominance of the LHB cratering. Regarding the size distribution of small disk objects, including the precursors of today’s JFCs, we note the recent finding from the New Horizons mission (Singer et al. 2015, 2016) that the impact craters on Pluto and Charon indicate a projectile CSD with a power-law slope of −2 – right in the middle of our assumed range.

We have thus established that there is no reason to cast doubt on a late-occurring planet instability in the Nice Model based on concerns that this would cause a too heavy cometary bombardment. One rather has to turn the question around and ask, if there is solid evidence that comets did still play an important role in the LHB, because this is clearly predicted by our results. We note in this regard that Marty et al. (2016) found the Ar abundance in the coma of comet 67P/Churyumov-Gerasimenko, combined with the assumption of a substantial cometary component of the LHB, would be consistent with the amount of Ar in the Earth’s atmosphere. Possibly, the water abundance measured in lunar mare basalts (Greenwood et al. 2011) could also derive from cometary LHB impacts.

Finally, we have estimated the total mass of all the cometary projectiles hitting the targets during the course of the LHB. For the Earth we get much less than the ocean mass. For Mars we also get much less than the estimated bulk water inventory (Lammer et al. 2013). The total impactor mass can be estimated at ~7 × 10¹⁸ kg, and using a H₂O abundance of ~20% in comets from Rosetta measurements (Rotundi et al. 2015), the Global Equivalent Layer is only about 10 m. For ~30% CO₂/H₂O by mass (Bockelée-Morvan et al. 2016), a CO₂ surface pressure of 90 mb might have been built up, which is insufficient to support liquid surface water.

We tentatively conclude that the water that flowed on Mars during the Noachian cannot owe its existence to the comets that impacted Mars during the LHB. If impacts were crucial to release subsurface H₂O reservoirs, these were mostly caused by asteroids. The water on both the Earth and Mars must in any case have been delivered long before the LHB – probably during the build-up phase of both planets (Morbidelli et al. 2000; Brasser 2013).

Acknowledgments

This work was supported by the Polish National Science Centre under Grant No. 2011/01/B/ST9/05442 as well as Grant No. 74/10:2 of the Swedish National Space Board. We gratefully acknowledge advice from Simone Marchi regarding the calculation of crater sizes and from Stephanie Werner concerning the observed crater areal densities. We thank the referee, David O’Brien, for valuable comments.

References

All Tables

All Figures

Fig. 1 Illustration of the MOID-chord method of calculating the random probability of collision with a target planet, given a known value of the minimum orbit intersection distance (MOID). The target’s collisional radius defines a circle in the b-plane, and the time it takes for the target’s projection to move the distance EE′ divided by its orbital period gives the probability. From Rickman et al. (2014) with permission.

In the text

	Fig. 2 Individual collision probabilities for random samples of 1000 comet orbits per planet, all of which satisfy the condition ℳ < Rcoll for the planet in question. For each planet, the orbits are ordered chronologically. The coloured symbols denote: Mercury (black circles), Venus (red squares), Earth (green triangles), and Mars (blue diamonds).
In the text

Fig. 3 The ratio (S) between the impact probabilities for the Earth and the Moon as found in test calculations (see the main text) is plotted versus the perihelion distance of the projectile for a constant inclination of 10°. Green dots: gravitational focussing neglected. Red and blue dots: gravitational focussing included, using the Moon’s ancient and contemporary orbital semi-major axis, respectively.

In the text

Fig. 4 Maximum cometary impactor radius for Mercury, the Moon and Mars. Each error bar extends from the lower to the upper quartile of the cumulative probability distribution with the median marked by a filled circle for the maximum model and a cross for the minimum model of the disk. The red colour is used for the reference model, the blue colour for the erosional model, and the magenta colour for the Di Sisto model. The shaded band for the Moon indicates an observational constraint based on the Jupiter Trojan population (see the main text).

In the text

	Fig. 5 Total mass of all the cometary impactors on Mercury, the Moon and Mars, plotted on a logarithmic scale. The error bars, and the colours and symbols, have the same meaning as in Fig. 4.
In the text

Fig. 6 Left panel: cumulative distributions of approach velocities to different planets for collisional orbits in the reference model. The curve for Mercury shows some wiggles owing to relatively sparse material. The thin vertical lines indicate the circular speeds at the mean distance of each planet. Right panel: cumulative distributions of approach velocities to the Earth for collisional orbits in the reference, erosional and Di Sisto models. For the latter two, the slope index s1 = 2.5 has been used. The wiggles in the curve for the Di Sisto model are due to the sparse material.

In the text

	Fig. 7 Maximum cometary impact basin diameter for Mercury, the Moon and Mars. The meaning of the “error bars”, markers and colours is the same as in Fig. 4. The grey lines indicate the diameters of the largest basins considered to be LHB-related on the three targets. The commonly used names of these basins are written above the respective lines.
In the text

Fig. 8 Areal density of D > 20 km craters for Mercury, the Moon and Mars. For each target the three physical evolution models are plotted using the same colours as in Fig. 4. The upper cross is for the steeper size distribution of the small disk objects, and the lower cross is for the shallower one. The low mass disk is used in both cases. The lines at the common level of 100 craters per million km2 mark the observational reference data.

In the text