© ESO, 2013

1. Introduction

The late heavy bombardment (LHB) is an important period in the history of the solar system. It is often defined as the process that made the huge but relatively young impact basins (a 300 km or larger diameter crater) on the Moon like Imbrium and Orientale. The sources and extent of the LHB, however, has been undergoing recent revisions. In the past, there were two end-member schools of thought describing the LHB. The first school argued that nearly all lunar basins, including the young ones, were made by impacting planetesimals left over from terrestrial planet formation (Neukum et al. 2001; Hartmann et al. 2000, 2007; see Chapman et al. 2007, for a review). The second school argued that most lunar basins were made during a spike of impacts that took place near 3.9 Ga (e.g., Tera et al. 1974; Ryder et al. 2000).

Recent studies, however, suggest that a compromise scenario may be the best solution: the oldest basins were mainly made by leftover planetesimals, while the last 12–15 or so lunar basins were created by asteroids driven out of the primordial main belt by the effects of late giant-planet migration (Tsiganis et al. 2005; Gomes et al. 2005; Minton & Malhotra 2009; Morbidelli et al. 2010; Marchi et al. 2012; Bottke et al. 2012). This would mean the LHB is limited in extent and does not encompass all lunar basins. If this view is correct, we can use studies of lunar and asteroid samples heated by impact events, together with dynamical modelling work, to suggest that the basin-forming portion of the LHB lasted from approximately 4.1–4.2 to 3.7–3.8 billion years ago on the Moon (Bogard 1995, 2011; Swindle et al. 2009; Bottke et al. 2012; Norman & Nemchin 2012).

The so-called “Nice model” provides a coherent explanation of the origin of the LHB as an impact spike or rather a “sawtooth” (Morbidelli et al. 2012). According to this model, the bombardment was triggered by a late dynamical orbital instability of the giant planets, in turn driven by the gravitational interactions between said planets and a massive trans-Neptunian disk of planetesimals (see Morbidelli 2010, for a review). In this scenario, three projectile populations contributed to the LHB: the comets from the original trans-Neptunian disk (Gomes et al. 2005), the asteroids from the main belt (Morbidelli et al. 2010) and those from a putative extension of the main belt towards Mars, inwards of its current inner edge (Bottke et al. 2012). The last could have been enough of a source for the LHB, as recorded in the lunar crater record (Bottke et al. 2012), while the asteroids from the current main belt boundaries would have only been a minor contributor (Morbidelli et al. 2010).

The Nice model, however, predicts a very intense cometary bombardment of which there seems to be no obvious traces on the Moon. In fact, given the expected total mass in the original trans-Neptunian disk (Gomes et al. 2005) and the size distribution of objects in this disk (Morbidelli et al. 2009), the Nice model predicts that about 5 × 10⁴ km-size comets should have hit the Moon during the LHB. This would have formed 20 km craters with a surface density of 1.7 × 10^-3 craters per km². But the highest crater densities of 20 km craters on the lunar highlands is less than 2 × 10^-4 (Strom et al. 2005). This discrepancy might be explained by a gross overestimate of the number of small bodies in the original trans-Neptunian disk in Morbidelli et al. (2009). However, all impact clast analyses of samples associated to major LHB basins (Kring & Cohen 2002; Tagle 2005) show that also the major projectiles were not carbonaceous chondrites or similar primitive, comet-like objects.

The lack of evidence of a cometary bombardment of the Moon can be considered as a fatal flaw in the Nice model. Curiously, however, in the outer solar system we see evidence of the cometary flux predicted by the Nice model. Such a flux is consistent with the number of impact basins on Iapetus (Charnoz et al. 2009), with the number and the size distribution of the irregular satellites of the giant planets (Nesvorný et al. 2007; Bottke et al. 2010) and of the Trojans of Jupiter (Morbidelli et al. 2005), as well as with the capture of D-type asteroids in the outer asteroid belt (Levison et al. 2009). Moreover, the Nice model cometary flux is required to explain the origin of the collisional break-up of the asteroid (153) Hilda in the 3/2 resonance with Jupiter (located at ≃4 AU, i.e. beyond the nominal outer border of the asteroid belt at ≃3.2 AU; Brož et al. 2011).

Missing signs of an intense cometary bombardment on the Moon and the evidence for a large cometary flux in the outer solar system suggest that the Nice model may be correct in its basic features, but most comets disintegrated as they penetrated deep into the inner solar system.

To support or reject this possibility, this paper focusses on the main asteroid belt, looking for constraints on the flux of comets through this region at the time of the LHB. In particular we focus on old asteroid families, produced by the collisional break-up of large asteroids, which may date back at the LHB time. We provide a census of these families in Sect. 2.

In Sect. 3, we construct a collisional model of the main belt population. We show that, on average, this population alone could not have produced the observed number of families with D_PB = 200–400 km. Instead, the required number of families with large parent bodies is systematically produced if the asteroid belt was crossed by a large number of comets during the LHB, as expected in the Nice model (see Sect. 4). However, for any reasonable size distribution of the cometary population, the same cometary flux that would produce the correct number of families with D_PB = 200–400 km would produce too many families with D_PB ≃ 100 km relative to what is observed. Therefore, in the subsequent sections we look for mechanisms that might prevent detection of most of these families.

More specifically, in Sect. 5 we discuss the possibility that families with D_PB ≃ 100 km are so numerous that they cannot be identified because they overlap with each other. In Sect. 6 we investigate their possible dispersal below detectability due to the Yarkovsky effect and chaotic diffusion. In Sect. 7 we discuss the role of the physical lifetime of comets. In Sect. 8 we analyse the dispersal of families due to the changes in the orbits of the giant planets expected in the Nice model. In Sect. 9 we consider the subsequent collisional comminution of the families. Of all investigated processes, the last one seems to be the most promising for reducing the number of visible families with D_PB ≃ 100 km while not affecting the detectability of old families with D_PB = 200–400 km.

Finally, in Sect. 10 we analyse a curious portion of the main belt, located in a narrow semi-major axis zone bounded by the 5:2 and 7:3 resonances with Jupiter. This zone is severely deficient in small asteroids compared to the other zones of the main belt. For the reasons explained in the section, we think that this zone best preserves the initial asteroid belt population, and therefore we call it the “pristine zone”. We checked the number of families in the pristine zone, their sizes, and ages and we found that they are consistent with the number expected in our model invoking a cometary bombardment at the LHB time and a subsequent collisional comminution and dispersion of the family members. The conclusions follow in Sect. 11.

2. A list of known families

Although several lists of families exist in the literature (Zappalá et al. 1995; Nesvorný et al. 2005; Parker et al. 2008; Nesvorný 2010), we are going to identify the families once again. The reason is that we seek an upper limit for the number of old families that may be significantly dispersed and depleted, while the previous works often focussed on well-defined families. Moreover, we need to calculate several physical parameters of the families (such as the parent-body size, slopes of the size-frequency distribution (SFD), a dynamical age estimate if not available in the literature) which are crucial for further modelling. Last but not least, we use more precise synthetic proper elements from the AstDyS database (Knežević & Milani 2003, version Aug. 2010) instead of semi-analytic ones.

We employed a hierarchical clustering method (HCM, Zappalá et al. 1995) for the initial identification of families in the proper element space (a_p,e_p,sinI_p), but then we had to perform a lot of manual operations, because i) we had to select a reasonable cut-off velocity v_cutoff, usually such that the number of members N(v_cutoff) increases relatively slowly with increasing v_cutoff. ii) The resulting family should also have a “reasonable” shape in the space of proper elements, which should somehow correspond to the local dynamical features¹. iii) We checked taxonomic types (colour indices from the Sloan DSS MOC catalogue version 4, Parker et al. 2008), which should be consistent among family members. We can recognise interlopers or overlapping families this way. iv) Finally, the SFD should exhibit one or two well-defined slopes, otherwise the cluster is considered uncertain.

Our results are summarised in online Table 1 and the positions of families within the main belt are plotted in Fig. 1. Our list is “optimistic”, so that even not very prominent families are included here².

There are, however, several potential problems we are aware of:

• 1.

  There may be inconsistencies among different lists of families. For example, sometimes a clump may be regarded as a single family or as two separate families. This may be the case of: Padua and Lydia, Rafita and Cameron.

• 2.

  To identify families we used synthetic proper elements, which are more precise than the semi-analytic ones. Sometimes the families look more regular (e.g., Teutonia) or more tightly clustered (Beagle) when we use the synthetic elements. This very choice may, however, affect results substantially! A clear example is the Teutonia family, which also contains the big asteroid (5) Astraea if the synthetic proper elements are used, but not if the semi-analytic proper elements are used. This is due to the large differences between the semi-analytic and synthetic proper elements of (5) Astraea. Consequently, the physical properties of the two families differ considerably. We believe that the family defined from the synthetic elements is more reliable.

• 3.

  Durda et al. (2007) often claim a larger size for the parent body (e.g., Themis, Meliboea, Maria, Eos, Gefion), because they try to match the SFD of larger bodies and the results of SPH experiments. This way they also account for small bodies that existed at the time of the disruption, but which do not exist today since they were lost due to collisional grinding and the Yarkovsky effect. We prefer to use D_Durda instead of the value D_PB estimated from the currently observed SFD. The geometric method of Tanga et al. (1999), which uses the sum of the diameters of the first and third largest family members as a first guess of the parent-body size, is essentially similar to our approach³.

2.1. A definition of the production function

To compare observed families to simulations, we define a “production function” as the cumulative number N(>D) of families with parent-body size D_PB larger than a given D. The observed production function is shown in Fig. 2, and it is worth noting that it is very shallow. The number of families with D_PB ≃ 100 km is comparable to the number of families in the D_PB = 200–400 km range.

It is important to note that the observed production function is likely to be affected by biases (the family sample may not be complete, especially below D_PB ≲ 100 km) and also by long-term collisional/dynamical evolution which may prevent a detection of old comminutioned/dispersed families today (Marzari et al. 1999).

From the theoretical point of view, the slope q of the production function N(>D) ∝ D^q should correspond to the cumulative slopes of the SFDs of the target and projectile populations. It is easy to show⁴ that the relation is $\mathit{q}\mathrm{=}\mathrm{2}\mathrm{+}{\mathit{q}}_{\mathrm{target}}\mathrm{+}\frac{\mathrm{5}}{\mathrm{3}}{\mathit{q}}_{\mathrm{project}}\mathit{.}$(1)Of course, real populations may have complicated SFDs, with different slopes in different ranges. Nevertheless, any populations that have a steep SFD (e.g. q_target = q_project = −2.5) would inevitably produce a steep production function (q ≐ −4.7).

In the following analysis, we drop cratering events and discuss catastrophic disruptions only, i.e. families which have largest remnant/parent body mass ratio less than 0.5. The reason is that the same criterion LR/PB < 0.5 is used in collisional models. Moreover, cratering events were not yet systematically explored by SPH simulations due to insufficient resolution (Durda et al. 2007).

2.2. Methods for family age determination

If there is no previous estimate of the age of a family, we used one of the following three dynamical methods to determine it: i) a simple (a_p,H) analysis as in Nesvorný et al. (2005); ii) a C-parameter distribution fitting as introduced by Vokrouhlický et al. (2006); iii) a full N-body simulation described e.g. in Brož et al. (2011).

In the first approach, we assume zero initial velocities, and the current extent of the family is explained by the size-dependent Yarkovsky semimajor axis drift. This way we can obtain only an upper limit for the dynamical age, of course. We show an example for the Eos family in Fig. 3. The extent of the family in the proper semimajor axis vs the absolute magnitude (a_p,H) plane can be described by the parametric relation $\mathrm{0.2}\mathit{H}\mathrm{=}{\log }_{\mathrm{10}}\frac{\mathrm{|}{\mathit{a}}_{\mathrm{p}}\mathrm{-}{\mathit{a}}_{\mathrm{c}}\mathrm{|}}{\mathit{C}}\mathit{,}$(2)where a_c denotes the centre of the family, and C is the parameter. Such relation can be naturally expected when the semimajor-axis drift rate is inversely proportional to the size, da/dt ∝ 1/D, and the size is related to the absolute magnitude via the Pogson equation $\mathit{H}\mathrm{=}\mathrm{-}\mathrm{2.5}{\log }_{\mathrm{10}}\mathrm{(}{\mathit{p}}_{\mathit{V}}{\mathit{D}}^{\mathrm{2}}\mathit{/}{\mathit{D}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{)}$, where D₀ denotes the reference diameter and p_V the geometric albedo (see Vokrouhlický et al. 2006 for a detailed discussion). The limiting value, for which all Eos family members (except interlopers) are above the corresponding curve, is C = 1.5to2.0 × 10^-4 AU. Assuming reasonable thermal parameters (summarised in Table 2), we calculate the expected Yarkovsky drift rates da/dt (using the theory from Brož 2006) and consequently can determine the age to be t < 1.5to2.0 Gyr.

The second method uses a histogram N(C,C + ΔC) of the number of asteroids with respect to the C parameter defined above, which is fitted by a dynamical model of the initial velocity field and the Yarkovsky/YORP evolution. This enables us to determine the lower limit for the age too (so the resulting age estimate is $\mathit{t}\mathrm{=}\mathrm{1.}{\mathrm{3}}_{-0.2}^{\mathrm{+}\mathrm{0.15}}\hspace{0.17em}\mathrm{Gyr}$ for the Eos family).

In the third case, we start an N-body simulation using a modified SWIFT integrator (Levison & Duncan 1994), with the Yarkovsky/YORP acceleration included, and evolve a synthetic family up to 4 Gyr. We try to match the shape of the observed family in all three proper orbital elements (a_p, e_p, sin I_p). In principle, this method may provide a somewhat independent estimate of the age. For example, there is a “halo” of asteroids in the surroundings of the nominal Eos family, which are of the same taxonomic type K, and we may fit the ratio N_halo/N_core of the number of objects in the “halo” and in the family “core” (Brož et al., in prep.).

The major source of uncertainty in all methods are unknown bulk densities of asteroids (although we use the most likely values for the S or C/X taxonomic classes, Carry 2012). The age scales approximately as t ∝ ρ_bulk. Nevertheless, we are still able to distinguish families that are young from those that are old, because the allowed ranges of densities for S-types (2 to 3 g/cm³) and C/X-types (1 to 2 g/cm³) are limited (Carry 2012) and so are the allowed ages of families.

2.3. Which families can be of LHB origin?

The ages of the observed families and their parent-body sizes are shown in Fig. 4. Because the ages are generally very uncertain, we consider that any family whose nominal age is older than 2 Gyr is potentially a family formed ~4 Gyr ago, i.e. at the LHB time. If we compare the number of “young” (<2 Gyr) and old families (>2 Gyr) with D_PB = 200–400 km, we cannot see a significant over-abundance of old family formation events. On the other hand, we almost do not find any small old families.

Only 12 families from the whole list may be possibly dated back to the LHB, because their dynamical ages approach ~3.8 Gyr (including the relatively large uncertainties; see Table 3, which is an excerpt from Table 1).

If we drop cratering events and the families of Camilla and Hermione, which do not exist any more today (their existence was inferred from the satellite systems, Vokrouhlický et al. 2010), we end up with only five families created by catastrophic disruptions that may potentially date from the LHB time (i.e. their nominal age is more than 2 Gy). As we shall see in Sect. 4, this is an unexpectedly low number.

Moreover, it is really intriguing that most “possibly-LHB” families are larger than D_PB ≃ 200 km. It seems that old families with D_PB ≃ 100 km are missing in the observed sample. This is an important aspect that we have to explain, because it contradicts our expectation of a steep production function.

3. Collisions in the main belt alone

Before we proceed to scenarios involving the LHB, we try to explain the observed families with ages spanning 0–4 Gyr as a result of collisions only among main belt bodies. To this purpose, we used the collisional code called Boulder (Morbidelli et al. 2009) with the following setup: the intrinsic probabilities P_i = 3.1 × 10^-18 km^-2 yr^-1, and the mutual velocities V_imp = 5.28 km s^-1 for the MB vs. MB collisions (both were taken from the work of Dahlgren 1998). The assumption of a single V_imp value is a simplification, but about 90% collisions have mutual velocities between 2 and 8 km s^-1 (Dahlgren 1998), which assures a similar collisional regime.

The scaling law is described by the polynomial relation (r denotes radius in cm) ${\mathit{Q}}_{\mathit{D}}^{\mathit{\star }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{{\mathit{q}}_{\mathrm{fact}}}\mathrm{(}{\mathit{Q}}_{\mathrm{0}}{\mathit{r}}^{\mathit{a}}\mathrm{+}\mathit{B\rho }{\mathit{r}}^{\mathit{b}}\mathrm{)}$(3)with the parameters corresponding to basaltic material at 5 km s^-1 (Benz & Asphaug 1999, Table 4):

Even though not all asteroids are basaltic, we use the scaling law above as a mean one for the main belt population. Below, we discuss also the case of significantly lower strengths (i.e. higher q_fact values).

We selected the time span of the simulation 4 Gyr (not 4.5 Gyr) since we are interested in this last evolutionary phase of the main belt, when its population and collisional activity is nearly same as today (Bottke et al. 2005). The outcome of a single simulation also depends on the “seed” value of the random-number generator that is used in the Boulder code to decide whether a collision with a fractional probability actually occurs or not in a given time step. We thus have to run multiple simulations (usually 100) to obtain information on this stochasticity of the collisional evolution process.

The initial SFD of the main belt population conditions was approximated by a three-segment power law (see also thin grey line in Fig. 5, 1st row) with differential slopes q_a = −4.3 (for D > D₁), q_b = −2.2, q_c = −3.5 (for D < D₂) where the size ranges were delimited by D₁ = 80 km and D₂ = 16 km. We also added a few big bodies to reflect the observed shape of the SFD at large sizes (D > 400 km). The normalisation was N_norm(D > D₁) = 350 bodies in this case.

We used the observed SFD of the main belt as the first constraint for our collisional model. We verified that the outcome our model after 4 Gyr is not sensitive to the value of q_c. Namely, a change of q_c by as much as ± 1 does not affect the final SFD in any significant way. On the other hand, the values of the remaining parameters (q_a, q_b, D₁, D₂, N_norm) are enforced by the observed SFD. To obtain a reasonable fit, they cannot differ much (by more than 5–10%) from the values presented above.

We do not use only a single number to describe the number of observed families (e.g. N = 20 for D_PB ≥ 100 km), but we discuss a complete production function instead. The results in terms of the production function are shown in Fig. 5 (left column, 2nd row). On average, the synthetic production function is steeper and below the observed one, even though there is approximately a 5% chance that a single realization of the computer model will resemble the observations quite well. This also holds for the distribution of D_PB = 200–400 km families in the course of time (age).

In this case, the synthetic production function of D_PB ≳ 100 km families is not significantly affected by comminution. According to Bottke et al. (2005), most of D > 10 km fragments survive intact and a D_PB ≳ 100 km family should be recognisable today. This is also confirmed by calculations with Boulder (see Fig. 5, left column, 3rd row).

To improve the match between the synthetic and the observed production function, we can do the following: i) modify the scaling law, or ii) account for a dynamical decay of the MB population. Using a substantially lower strength (q_fact = 5 in Eq. (3), which is not likely, though) one can obtain a synthetic production function which is on average consistent with the observations in the D_PB = 200–400 km range.

Regarding the dynamical decay, Minton & Malhotra (2010) suggest that initially the MB was three times more populous than today while the decay timescale was very short: after 100 Myr of evolution the number of bodies is almost at the current level. In this brief period of time, about 50% more families will be created, but all of them will be old, of course. For the remaining  ~3.9 Gyr, the above model (without any dynamical decay) is valid.

To conclude, it is possible – though not very likely – that the observed families were produced by the collisional activity in the main belt alone. A dynamical decay of the MB population would create more families that are old, but technically speaking, this cannot be distinguished from the LHB scenario, which is discussed next.

4. Collisions between a “classical” cometary disk and the main belt

In this section, we construct a collisional model and estimate an expected number of families created during the LHB due to collisions between cometary-disk bodies and main belt asteroids. We start with a simple stationary model and we confirm the results using a more sophisticated Boulder code (Morbidelli et al. 2009).

Using the data from Vokrouhlický et al. (2008) for a “classical” cometary disk, we can estimate the intrinsic collisional probability and the collisional velocity between comets and asteroids. A typical time-dependent evolution of P_i and V_imp is shown in Fig. 6. The probabilities increase at first, as the trans-Neptunian cometary disk starts to decay, reaching up to 6 × 10^-21 km^-2 yr^-1, and after 100 Myr they decrease to zero. These results do not differ significantly from run to run.

4.1. Simple stationary model

In a stationary collisional model, we choose an SFD for the cometary disk, we assume a current population of the main belt; estimate the projectile size needed to disrupt a given target according to (Bottke et al. 2005) ${\mathit{d}}_{\mathrm{disrupt}}\mathrm{=}{\mathrm{(}{\mathrm{2}{\mathit{Q}}_{\mathit{D}}^{\mathit{\star }}\mathit{/}{\mathit{V}}_{\mathrm{imp}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{D}}_{\mathrm{target}}\mathit{,}$(4)where ${\mathit{Q}}_{\mathit{D}}^{\mathit{\star }}$ denotes the specific energy for disruption and dispersion of the target (Benz & Asphaug 1999); and finally calculate the number of events during the LHB as ${\mathit{n}}_{\mathrm{events}}\mathrm{=}\frac{{\mathit{D}}_{\mathrm{target}}^{\mathrm{2}}}{\mathrm{4}}\hspace{0.17em}{\mathit{n}}_{\mathrm{target}}\mathrm{\int }{\mathit{P}}_{\mathrm{i}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}{\mathit{n}}_{\mathrm{project}}\mathrm{\left(}\mathit{t}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{t,}$(5)where n_target and n_project are the number of targets (i.e. main belt asteroids) and the number of projectiles (comets), respectively. The actual number of bodies (27 000) in the dynamical simulation of Vokrouhlický et al. (2008) changes in the course of time, and it was scaled such that it was initially equal to the number of projectiles N(>d_disrupt) inferred from the SFD of the disk. This is clearly a lower limit for the number of families created, since the main belt was definitely more populous in the past.

The average impact velocity is V_imp ≃ 10 km s^-1, so we need the projectile sizes to disrupt given target sizes listed in Table 5.

We tried to use various SFDs for the cometary disk (i.e., with various differential slopes q₁ for D > D₀ and q₂ for D < D₀, the elbow diameter D₀ and total mass M_disk), including rather extreme cases (see Fig. 7). The resulting numbers of LHB families are summarised in Table 6. Usually, we obtain several families with D_PB ≃ 200 km and about 100 families with D_PB ≃ 100 km. This result is robust with respect to the slope q₂, because even very shallow SFDs should produce a lot of these families⁵. The only way to decrease the number of families significantly is to assume the elbow at a larger diameter D₀ ≃ 150 km.

It is thus no problem to explain the existence of approximately five large families with D_PB = 200–400 km, which are indeed observed, since they can be readily produced during the LHB. On the other hand, the high number of D_PB ≃ 100 km families clearly contradicts the observations, since we observe almost no LHB families of this size.

4.2. Constraints from (4) Vesta

The asteroid (4) Vesta presents a significant constraint for collisional models, being a differentiated body with a preserved basaltic crust (Keil 2002) and a 500 km large basin on its surface (a feature indicated by the photometric analysis of Cellino et al. 1987), which is significantly younger than 4 Gyr (Marchi et al. 2012). It is highly unlikely that Vesta experienced a catastrophic disruption in the past, and even large cratering events were limited. We thus have to check the number of collisions between one D = 530 km target and D ≃ 35 km projectiles, which are capable of producing the basin and the Vesta family (Thomas et al. 1997). According to Table 6, the predicted number of such events does not exceed ~2, so given the stochasticity of the results there is a significant chance that Vesta indeed experienced zero such impacts during the LHB.

4.3. Simulations with the Boulder code

To confirm results of the simple stationary model, we also performed simulations with the Boulder code. We modified the code to include a time-dependent collisional probability P_i(t) and impact velocity V_imp(t) of the cometary-disk population.

We started a simulation with a setup for the cometary disk resembling the nominal case from Table 6. The scaling law is described by Eq. (3) with the parameters given in Table 7, suitable for asteroids (basalt) and comets (water ice).

The intrinsic probabilities P_i = 3.1 × 10^-18 km^-2 yr^-1 and velocities V_imp = 5.28 km s^-1 for the MB vs MB collisions were again taken from the work of Dahlgren (1998). We do not account for comet-comet collisions since their evolution is dominated by the dynamical decay. The initial SFD of the main belt was similar to the one in Sect. 3, q_a = −4.2, q_b = −2.2, q_c = −3.5, D₁ = 80 km, D₂ = 14 km, and only the normalisation was increased up to N_norm(D > D₁) = 560 in this case.

The resulting SFDs of 100 independent simulations with different random seeds are shown in Fig. 5 (middle column). The number of LHB families (approximately 10 with D_PB ≃ 200 km and 200 with D_PB ≃ 100 km) is even larger compared to the stationary model, as expected, because we had to start with a larger main belt to get a good fit of the currently observed MB after 4 Gyr of collisional evolution.

To conclude, the stationary model and the Boulder code give results that are compatible with each other, but that clearly contradict the observed production function of families. In particular, they predict far too many families with D = 100 km parent bodies. At first sight, this may be interpreted as proof that there was no cometary LHB on the asteroids. Before jumping to this conclusion, however, one has to investigate whether there are biases against identifying of D_PB = 100 km families. In Sects. 5–9 we discuss several mechanisms that all contribute, at some level, to reducing the number of observable D_PB = 100 km families over time. They are addressed in order of relevance, from the least to the most effective.

5. Families overlap

Because the number of expected D_PB ≥ 100 km LHB families is very high (of the order of 100) we now want to verify if these families can overlap in such a way that they cannot be distinguished from each other and from the background. We thus took 192 main belt bodies with D ≥ 100 km and selected randomly 100 of them that will break apart. For each one we created an artificial family with 10² members, assume a size-dependent ejection velocity V ∝ 1/D (with V = 50 m/s for D = 5 km) and the size distribution resembling that of the Koronis family. The choice of the true anomaly and the argument of perihelion at the instant of the break-up event was random. We then calculated proper elements (a_p,e_p,sinI_p) for all bodies. This type of analysis is in some respects similar to the work of Bendjoya et al. (1993).

According to the resulting Fig. 8 the answer to the question is simple: the families do not overlap sufficiently, and they cannot be hidden that way. Moreover, if we take only bigger bodies (D > 10 km), these would be clustered even more tightly. The same is true for proper inclinations, which are usually more clustered than eccentricities, so families could be more easily recognised.

6. Dispersion of families by the Yarkovsky drift

In this section, we model long-term evolution of synthetic families driven by the Yarkovsky effect and chaotic diffusion. For one synthetic family located in the outer belt, we have performed a full N-body integration with the SWIFT package (Levison & Duncan 1994), which includes also an implementation of the Yarkovsky/YORP effect (Brož 2006) and second-order integrator by Laskar & Robutel (2001). We included 4 giant planets in this simulation. To speed-up the integration, we used ten times smaller sizes of the test particles and thus a ten times shorter time span (400 Myr instead of 4 Gyr). The selected time step is Δt = 91 d. We computed proper elements, namely their differences Δa_p,Δe_p,ΔsinI_p between the initial and final positions.

Then we used a simple Monte-Carlo approach for the whole set of 100 synthetic families – we assigned a suitable drift Δa_p(D) in semimajor axis, and also drifts in eccentricity Δe_p and inclination ΔsinI_p to each member of 100 families, respecting asteroid sizes, of course. This way we account for the Yarkovsky semimajor axis drift and also for interactions with mean-motion and secular resonances. This Monte-Carlo method tends to smear all structures, so we can regard our results as the upper limits for dispersion of families.

While the eccentricities of small asteroids (down to D ≃ 2 km) seem to be dispersed enough to hide the families, there are still some persistent structures in inclinations, which would be observable today. Moreover, large asteroids (D ≥ 10 km) seem to be clustered even after 4 Gyr, so that more than 50% of families can be easily recognised against the background (see Fig. 9). We thus can conclude that it is not possible to disperse the families by the Yarkovsky effect alone.

7. Reduced physical lifetime of comets in the MB crossing zone

To illustrate the effects that the physical disruption of comets (due to volatile pressure build-up, amorphous/crystalline phase transitions, spin-up by jets, etc.) can have on the collisional evolution of the asteroid belt, we adopted here a simplistic assumption. We considered that no comet disrupt beyond 1.5 AU, whereas all comets disrupt the first time that they penetrate inside 1.5 AU. Both conditions are clearly not true in reality: some comets are observed to blow up beyond 1.5 AU, and others are seen to survive on an Earth-crossing orbit. Thus we adopted our disruption law just as an example of a drastic reduction of the number of comets with small perihelion distance, as required to explain the absence of evidence for a cometary bombardment on the Moon.

We then removed all those objects from output of comet evolution during the LHB that had a passage within 1.5 AU from the Sun, from the time of their first passage below this threshold. We then recomputed the mean intrinsic collision probability of a comet with the asteroid belt. The result is a factor ~3 smaller than when no physical disruption of comets is taken into account as in Fig. 6. The mean impact velocity with asteroids also decreases, from 12 km s^-1 to 8 km s^-1.

The resulting number of asteroid disruption events is thus decreased by a factor ~4.5, which can be also seen in the production function shown in Fig. 5 (right column). The production of families with D_PB = 200–400 km is consistent with observations, while the number of D_PB ≃ 100 km families is reduced to 30–70, but is still too high, by a factor 2–3. More importantly, the slope of the production function remains steeper than that of the observed function. Thus, our conclusion is that physical disruptions of comets alone cannot explain the observation, but may be an important factor to keep in mind for reconciling the model with the data.

8. Perturbation of families by migrating planets (a jumping-Jupiter scenario)

In principle, families created during the LHB may be perturbed by still-migrating planets. It is an open question what the exact orbital evolution of planets was at that time. Nevertheless, a plausible scenario called a “jumping Jupiter” was presented by Morbidelli et al. (2010). It explains major features of the main belt (namely the paucity of high-inclination asteroids above the ν₆ secular resonance), and is consistent with amplitudes of the secular frequencies of both giant and terrestrial planets and also with other features of the solar system. In this work, we thus investigated this particular migration scenario.

We used the data from Morbidelli et al. (2010) for the orbital evolution of giant planets. We then employed a modified SWIFT integrator, which read orbital elements for planets from an input file and calculated only the evolution of test particles. Four synthetic families located in the inner/middle/outer belt were integrated. We started the evolution of planets at various times, ranging from t₀ to (t₀ + 4 Myr) and stopped the integration at (t₀ + 4 My), in order to test the perturbation on families created in different phases of migration. Finally, we calculated proper elements of asteroids when the planets do not migrate anymore. (We also had to move planets smoothly to their exact current orbital positions.)

The results are shown in Fig. 10. While the proper eccentricities seem to be sufficiently perturbed and families are dispersed even when created at late phases of migration, the proper inclinations are not very dispersed, except for families in the outer asteroid belt that formed at the very beginning of the giant planet instability (which may be unlikely, as there must be a delay between the onset of planet instability and the beginning of the cometary flux through the asteroid belt). In most cases, the LHB families could still be identified as clumps in semi-major axis vs inclination space. We do not see any of such (a_p,sinI_p)-clumps, dispersed in eccentricity, in the asteroid belt⁶.

The conclusion is clear: it is not possible to destroy low-e and low-I families by perturbations arising from giant-planet migration, at least in the case of the “jumping-Jupiter” scenario⁷.

9. Collisional comminution of asteroid families

We have already mentioned that the comminution is not sufficient to destroy a D_PB = 100 km family in the current environment of the main belt (Bottke et al. 2005).

However, the situation in case of the LHB scenario is different. Both the large population of comets and the several-times larger main belt, which has to withstand the cometary bombardment, contribute to the enhanced comminution of the LHB families. To estimate the amount of comminution, we performed the following calculations: i) for a selected collisional simulation, whose production function is close to the average one, we recorded the SFDs of all synthetic families created in the course of time; ii) for each synthetic family, we restarted the simulation from the time t₀ when the family was crated until 4 Gyr and saved the final SFD, i.e. after the comminution. The results are shown in Fig. 11.

It is now important to discuss criteria, which enable us to decide if the comminutioned synthetic family would indeed be observable or not. We use the following set of conditions: D_PB ≥ 50 km, D_LF ≥ 10 km (largest fragment is the first or the second largest body, where the SFD becomes steep), LR/PB < 0.5 (i.e. a catastrophic disruption). Furthermore, we define N_members as the number of the remaining family members larger than observational limit D_limit ≃ 2 km and use a condition N_members ≥ 10. The latter number depends on the position of the family within the main belt, though. In the favourable “almost-empty” zone (between a_p = 2.825and2.955 AU), N_members ≥ 10 may be valid, but in a populated part of the MB one would need N_members ≳ 100 to detect the family. The size distributions of synthetic families selected this way resemble the observed SFDs of the main belt families.

According to Fig. 5 (3rd row), where we can see the production functions after comminution for increasing values of N_members, families with D_PB = 200–400 km remain more prominent than D_PB ≃ 100 km families simply because they contain much more members with D > 10 km that survive intact. Our conclusion is thus that comminution may explain the paucity of the observed D_PB ≃ 100 km families.

10. “Pristine zone” between the 5:2 and 7:3 resonances

We now focus on the zone between the 5:2 and 7:3 mean-motion resonances, with a_p = 2.825to2.955 AU, which is not as populated as the surrounding regions of the main belt (see Fig. 1). This is a unique situation, because both bounding resonances are strong enough to prevent any asteroids from outside to enter this zone owing the Yarkovsky semimajor axis drift. Any family formation event in the surroundings has only a minor influence on this narrow region. It thus can be called “pristine zone” because it may resemble the belt prior to creation of big asteroid families.

We identified nine previously unknown small families that are visible on the (e_p,sinI_p) plot (see Fig. 12). They are confirmed by the SDSS colours and WISE albedos, too. Nevertheless, there is only one big and old family in this zone (D_PB ≥ 100 km), i.e. Koronis.

That at most one LHB family (Koronis) is observed in the “pristine zone” can give us a simple probabilistic estimate for the maximum number of disruptions during the LHB. We take the 192 existing main belt bodies which have D ≥ 100 km and select randomly 100 of them that will break apart. We repeat this selection 1000 times and always count the number of families in the pristine zone. The resulting histogram is shown in Fig. 13. As we can see, there is very low (<0.001) probability that the number of families in the pristine zone is zero or one. On average we get eight families there, i.e. about half of the 16 asteroids with D ≥ 100 km present in this zone. It seems that either the number of disruptions should be substantially lower than 100 or we expect to find at least some “remnants” of the LHB families here.

It is interesting that the SFD of an old comminutioned family is very flat in the range D = 1to10 km (see Fig. 11) – similar to those of some of the “less certain” observed families! We may speculate that the families like (918) Itha, (5567) Durisen, (12573) 1999 NJ₅₃, or (15454) 1998 YB₃ (all from the pristine zone) are actually remnants of larger and older families, even though they are denoted as young. It may be that the age estimate based on the (a_p,H) analysis is incorrect since small bodies were destroyed by comminution and spread by the Yarkovsky effect too far away from the largest remnant, so they can no longer be identified with the family.

Finally, we have to ask an important question: what does an old/comminutioned family with D_PB ≃ 100 km look like in proper-element space? To this aim, we created a synthetic family in the “pristine zone”, and assumed the family has N_members ≃ 100 larger than D_limit ≃ 2 km and that the SFD is already flat in the D = 1to10 km range. We evolved the asteroids up to 4 Gyr due to the Yarkovsky effect and gravitational perturbations, using the N-body integrator as in Sect. 6. Most of the D ≃ 2 km bodies were lost in the course of the dynamical evolution, of course. The resulting family is shown in Fig. 14. We can also imagine that this family is placed in the pristine zone among other observed families, to get a feeling of whether it is easily observed or not (refer to Fig. 12).

It is clear that such family is hardly observable even in the almost empty zone of the main belt! Our conclusion is that the comminution (as given by the Boulder code) can explain the paucity of D_PB ≃ 100 km LHB families, since we can hardly distinguish old families from the background.

11. Conclusions

In this paper we investigated the cometary bombardment of the asteroid belt at the time of the LHB, in the framework of the Nice model. There is much evidence of a high cometary flux through the giant planet region, but no strong evidence of a cometary bombardment on the Moon. This suggests that many comets broke up on their way to the inner solar system. By investigating the collisional evolution of the asteroid belt and comparing the results to the collection of actual collisional families, our aim was to constrain whether the asteroid belt experienced an intense cometary bombardment at the time of the LHB and, if possible, constrain the intensity of this bombardment.

Observations suggest that the number of collisional families is a very shallow function of parent-body size (that we call in this paper the “production function”). We show that the collisional activity of the asteroid belt as a closed system, i.e. without any external cometary bombardment, in general does not produce such a shallow production function. Moreover, the number of families with parent bodies larger than 200 km in diameter is in general too small compared to the observations. However, there is a lot of stochasticity in the collisional evolution of the asteroid belt, and about 5% of our simulations actually fit the observational constraints (shallowness of the production function and number of large families) quite well. Thus, in principle, there is no need for a bombardment due to external agents (i.e. the comets) to explain the asteroid family collection, provided that the real collisional evolution of the main belt was a “lucky” one and not the “average” one.

If one accounts for the bombardment provided by the comets crossing the main belt at the LHB time, predicted by the Nice model, one can easily justify the number of observed families with parent bodies larger than 200 km. However, the resulting production function is steep, and the number of families produced by parent bodies of 100 km is almost an order of magnitude too large.

We have investigated several processes that may decimate the number of families identifiable today with 100 km parent bodies, without considerably affecting the survival of families formed from larger parent bodies. Of all these processes, the collisional comminution of the families and their dispersal by the Yarkovsky effect are the most effective ones. Provided that the physical disruption of comets due to activity reduced the effective cometary flux through the belt by a factor of ≈5, the resulting distribution of families (and consequently the Nice model) is consistent with observations.

To better quantify the effects of various cometary-disruption laws, we computed the numbers of asteroid families for different critical perihelion distances q_crit and for different disruption probabilities p_crit of comets during a given time step (Δt = 500 yr in our case). The results are summarised in Fig. 15. Provided that comets are disrupted frequently enough, namely the critical perihelion distance has to be at least q_crit ≳ 1 AU, while the probability of disruption is p_crit = 1, the number of D_PB ≥ 100 km families drops by the aforementioned factor of ≈5. Alternatively, q_crit may be larger, but then comets have to survive multiple perihelion passages (i.e. p_crit have to be lower than 1). It would be very useful to test these conditions by independent models of the evolution and physical disruptions of comets. Such additional constraints on cometary-disruption laws would then enable study of the original SFD of the cometary disk in more detail.

We can also think of two “alternative” explanations: i) physical lifetime of comets was strongly size-dependent so that smaller bodies break up easily compared to bigger ones; ii) high-velocity collisions between hard targets (asteroids) and very weak projectiles (comets) may result in different outcomes than in low-velocity regimes explored so far. Our work thus may also serve as a motivation for further SPH simulations.

We finally emphasize that any collisional/dynamical models of the main asteroid belt would benefit from the following advances:

• i)

  determination of reliable masses of asteroids of various classes. This may be at least partly achieved by the Gaia mission in the near future. Using up-to-date sizes and shape models (volumes) of asteroids one can then derive their densities, which are directly related to ages of asteroid families.

• ii)

  Development of methods for identifying asteroid families and possibly targeted observations of larger asteroids addressing their membership, which is sometimes critical for constructing SFDs and for estimating parent-body sizes.

• iii)

  An extension of the SHP simulations for both smaller and larger targets, to assure that the scaling we use now is valid. Studies and laboratory measurements of equations of states for different materials (e.g. cometary-like, porous) are closely related to this issue.

The topics outlined above seem to be the most urgent developments to be pursued in the future.

Online material

---

Acknowledgments

The work of M.B. and D.V. has been supported by the Grant Agency of the Czech Republic (grants no. 205/08/0064, and 13-01308S) and the Research Programme MSM0021620860 of the Czech Ministry of Education. The work of W.B. and D.N. was supported by NASA’s Lunar Science Institute (Center for Lunar Origin and Evolution, grant number NNA09DB32A). We also acknowledge the usage of computers of the Observatory and Planetarium in Hradec Králové. We thank Alberto Cellino for a careful review which helped to improve final version of the paper.

References

All Tables

All Figures