© ESO, 2015

1. Introduction

Comet nuclei are usually thought to be icy planetesimals, formed beyond the snow line in the nascent solar system. As such, they are naturally considered as precious targets of space missions – e.g., Rosetta. This concept is supported by the properties of comet nuclei derived from observations. The low bulk densities (e.g., Rickman 1989; Davidsson et al. 2007), the negligible tensile strength inferred for comet D/Shoemaker-Levy 9 (Asphaug & Benz 1996), and the low tensile strength of the surface layer required to explain their activity (Blum et al. 2014) are consistent with low-velocity accretion, in line with the general expectation for planetesimals formed at large distance from the Sun. Many comets have shown evidence for significant contributions of the super-volatile CO molecule to their outgassing activity (e.g., Bockelée-Morvan et al. 2004). This may be taken as an indication of a chemically pristine nature of the material that comet nuclei are made of, which supports the idea of a very gentle accretion process.

However, the problem of collisional evolution in the population of icy planetesimals has also been discussed in several papers. Davis & Farinella (1997) modeled the collisional evolution of the Edgeworth-Kuiper belt (EKB), which at the time was thought to be the source region for the Jupiter family comets (JFCs). They found that while large EKB members are most likely primordial objects, those with radii of a few km, like typical JFCs, are multigenerational fragments formed by the splitting of larger objects (see also Schlichting et al. 2013). Simultaneous with this first investigation, the scattered disk was discovered (Luu et al. 1997), and it was rapidly recognized to be a more efficient source of JFCs than the EKB. As a result of the longer orbital periods of its typical orbits, the scattered disk is believed to be less collisionally evolved (Rickman 2004), so that it is expected that the observed JFCs have a higher probability of being primordial planetesimals than in the analysis of Farinella & Davis (1997).

Another scenario was considered by Stern & Weissman (2001): the formation of the Oort cloud. This was modeled in the classical picture of gravitational ejection of icy planetesimals from the growth region of the giant planets (Safronov 1977). Stern and Weissman showed that during the course of this process, comet nuclei would be destroyed by collisions with small-scale debris. The authors concluded that the storage into the Oort cloud would be delayed until the comet source region had been cleared of material so that the collisional lifetime becomes longer than the ejection lifetime. Naturally, most of the Oort cloud comets would still bear the marks of collisional erosion in this scenario. A similar conclusion would also apply to the origin of the scattered disk.

Charnoz & Morbidelli (2003, 2007; CM03/07 hereafter) introduced a new algorithm for evaluating the effects of collisions in a population of small bodies that are subject to a complex and rapid dynamical evolution through gravitational perturbations, as is the case for planetesimals ejected from the region of the giant planets. This replaced the particle-in-a-box models used before. CM03/07 showed that for some appropriate initial size distributions (similar to that currently observed in the trans-Neptunian region), a sufficiently large number of comet-sized bodies would have reached the Oort cloud and the scattered disk. Although this was not explicitly discussed in these papers, they found that the vast majority of the Oort cloud objects that are larger than 1 km in radius would be pristine planetesimals. However, in the scattered disk, only 2% of the final population of objects of this size would be primordial, the rest would be collisional fragments.

Since these papers were produced, the Nice model for the early evolution of the solar system has been introduced (Tsiganis et al. 2005; for the latest version, see Levison et al. 2011). This changes the picture of the origin and evolution of comets in important ways. One central concept is the transplanetary disk of icy planetesimals. This was the disk of objects formed during the infant stages of the solar system beyond the original orbits of all giant planets, which were originally closer to the Sun. This disk extended out to about 30 AU and had a total mass of 20–50 Earth masses. It remained relatively quiescent until it was eventually dispersed as a consequence of a dynamical instability among the giant planets and of the planets’ subsequent migration toward their current orbits. The trans-planetary disk, upon its dispersal, is thought to have given rise to both the scattered disk and the Oort cloud (Brasser & Morbidelli 2013). Thus, this disk may once have been the repository for all the comets observed today. This would be compatible with the lack of evidence for any clear-cut differences in molecular composition (A’Hearn et al. 2012) or D/H ratios (Altwegg et al. 2015) between JFCs and comets of Oort cloud provenance, that is, long-period comets (LPCs) and Halley-type comets (HTCs). There currently is no alternative model capable of fully explaining the structure of the outer solar system that does not invoke an instability of the giant planets associated with the dispersal of the transplanetary disk that is similar to the instability in the concept of the Nice model.

In this paper we investigate the collision rates involving the members of the transplanetary disk during its whole evolution. First, we make the standard assumption that the dispersal of the disk coincided with the beginning of the so-called Late Heavy Bombardment (Gomes et al. 2005; Morbidelli et al. 2012), which means a lifetime of about 450 My for the disk before its dynamical dispersal. The dynamical state of the disk was most likely significantly excited. In fact, the probability of an object surviving the dynamical dispersal of the disk, remaining permanently trapped into the hot EKB population (including mean-motion resonances with Neptune), is lower than 10^-3 (see Nesvorný 2015, for the most updated estimate); this means that about 1000 Pluto-sized objects probably existed in the primordial disk (Stern 1991). These bodies would have induced a significant excitation in the disk, causing a velocity dispersion on the order of 0.5−1 km s^-1 (Levison et al. 2011). If the disk remained in this state for hundreds of millions of years, the collisional evolution of comet-sized objects has probably been severe. We quantify this in Sect. 3.1. The conclusions apply to comets in the scattered disk and in the Oort cloud because the transplanetary disk was the progenitor of both these reservoirs (Brasser & Morbidelli 2013).

However, it is not yet certain that the dynamical dispersal of the transplanetary disk occurred late. The formation of the latest basins on the Moon (Imbrium and Orientale, and possibly all Nectarian basins) requires that new projectiles appeared in the terrestrial planet-crossing region several 100 My after terrestrial planet formation (Bottke et al. 2007). The late instability of the giant planets would do this in a natural way (Bottke et al. 2012; Morbidelli et al. 2012). However, alternative models have been proposed. Some (Ćuk 2012; Minton et al. 2015) invoked unlikely spectacular collision events in the inner solar system, generating a large amount of debris that would subsequently have bombarded the Moon and the terrestrial planets. These models have not yet been thoroughly tested against all available solar system constraints. Other models, such as the destabilization of a population of lunar coorbital objects (Ćuk & Gladman 2006) or of a fifth terrestrial planet that would then have partially destabilized the asteroid belt (Chambers 2007) did not pass such tests (Ćuk & Gladman 2009; Brasser & Morbidelli 2011). The main argument independent of the lunar crater record in favor of a late dispersal of the trans-planetary disk is that the impact basins on Iapetus (a satellite of Saturn) have topographies that have relaxed by 25% or less, which argues that the surface layer of Iapetus was already very viscous at the time of basin formation. According to models of the thermal evolution of the satellite, this high viscosity could not be possible earlier than 200 My after the beginning of the solar system (Robuchon et al. 2011), which implies that the basins of Iapetus formed late. Nevertheless, this constraint remains model dependent.

Thus, in the second part of the paper we consider the case of an early dispersal of the transplanetary disk, occurring just after gas removal. In this case, the collisional evolution before the instability would be negligible (the disk would have survived just a few My and its planetesimals would presumably have had a very low velocity dispersion because of gas drag); however, during the dispersal of the disk, the collisional evolution might have been severe (similar to the case studied by Stern & Weissman 2001). In Sect. 3.2 we quantify the collisional evolution of comet-sized bodies that might eventually have been stored in the scattered disk during such a dispersal. Obviously, if the transplanetary disk dispersed late, the real collisional evolution of comets now in the scattered disk would be the sum of those studied in Sects. 3.1 and 3.2.

This paper is structured as follows. We start by presenting the logic of our reasoning, our methods, assumptions, and choice of parameters in Sect. 2; the results are presented in Sects. 3.1 and 3.2; a discussion and the summary of our conclusions are given in Sect 4.

2. Methods and principles

We here follow the approach of CM03/07, but with some important variants detailed in this section. The approach of CM03/07 is very suitable for studying the collisional evolution of a population of bodies undergoing a significant dynamical evolution, with orbital histories that can be quite different from one particle to the other. It consists of three steps. First, the dynamical history of each particle is followed through a numerical simulation. Second, the intrinsic collision probability and impact velocity of each particle with all others at each output time-step is computed from the output of the numerical simulation. Third, each particle in the simulation is assumed to be a tracer of a swarm of particles with a given initial size distribution. Then, using the information computed in the second step, the size distributions associated with each tracer are derived from one output step to another. This involves computing the minimal projectile size for a catastrophic impact on targets of any given size, and the size distribution of the generated fragments.

Below, we detail how we performed each of these three steps, and in particular, how we simplified the third step to reduce the parameter space we need to explore while still satisfying our needs and achieving our goals.

2.1. Numerical simulations

We used two pre-existing simulations that represent well the two phases of the evolution of the transplanetary disk described in the Introduction: the pre-instability phase and the dynamical dispersal phase.

2.1.1. Pre-instability disk

For the pre-instability phase we used one of the simulations of Levison et al. (2011). In these simulations, the planet instability occurs late and the disk is modeled in such a way as to mimic the self-excitation it would suffer if it contained 1000 Pluto-mass bodies. As mentioned in the Introduction, this is a realistic number for the bodies of this reference mass in the original transplanetary disk. We refer to Sect. 3 of Levison et al. for a technical description of the simulation and to Fig. 2 of that paper for a test showing that the self-excitation process is properly reproduced.

We took the state of the disk (i.e., the orbital distribution of the particles) 300 My after the beginning of the simulation. The self-stirring process increases the orbital excitation as $\sqrt{\mathit{t}}$ so that the disk is excited very rapidly, during the first few My, and then the evolution of the excitation slows down. Thus, we take the orbital distribution of the disk at 300 My in the simulation as representative of the real dynamical state of the disk during most of its pre-instability history (we test how the results change if the disk state is taken at an earlier time in Sect. 3.1). The top panel of Fig. 1 shows the (a,e) distribution we considered. As seen, there is a clear gradient of excitation with semi-major axis. This is because (i) the shorter orbital periods in the inner part of the disk produce more frequent encounters with the massive bodies and (ii) the orbital density of the massive bodies is higher (in this simulation the initial surface density of the population of bodies in the disk is assumed to be proportional to 1 /r)¹.

Fig. 1 Top: semi-major axis vs. eccentricity distribution of the transplanetary disk under the stirring effect of an embedded population of 1000 Pluto-mass bodies according to Levison et al. (2011). This snapshot of the distribution is taken after 300 My of evolution. Bottom: semi-major axis vs. eccentricity distribution of the scattered disk produced by the dispersal of the transplanetary disk by the giant planet instability according to Gomes et al. (2005). Here the snapshot of the scattered disk orbital distribution is taken 350 My after the beginning of the planet instability, when the scattered disk contains 5% of the original disk’s particles.

The most advanced estimate for the time of the instability, achieved by calibrating the Nice model on various constraints (Bottke et al. 2012; Morbidelli et al. 2012; Marchi et al. 2013) is ~4.1 Gy ago, namely 450 My after the disappearance of the gas from the protoplanetary disk (4.56 Gy ago). Given the uncertainty on this estimate and to remain conservative, we assume in the following that the pre-instability phase of the disk lasts 400 My.

2.1.2. Dispersal of the disk and origin of the scattered disk

To study the dispersal of the transplanetary disk and the formation of the scattered disk, we used one of the simulations presented in Gomes et al. (2005). In that simulation the planets become unstable at 887 My instead of at the preferred date of ~450 My. We considered only the dynamical histories of the particles after the beginning of the instability, ignoring the previous evolution, so that the actual instability date in the simulation has no influence on our results. In fact, the dispersal of the disk in the Nice model is a very violent event, and thus the memory of what occurred in the pre-instability phase is quickly erased. In addition, the overall evolution of the planetesimal population is quite insensitive to the exact evolution of the giant planets’ orbits during the instability. This is shown by the fact that radically different simulations provide scattered disk populations that decay in time in very similar ways down to ~1% of the original disk population, as reviewed in Fig. 5 of Brasser & Morbidelli (2013). Thus we think that the simulation that we consider is sufficient for our purposes.

The simulation covers a time-span of 350 My after the instability, identifies each particle individually, and produces a scattered disk made of 5% of the original particles at this date, whose (a,e) distribution is depicted in the bottom panel of Fig. 1.

2.2. Collision probabilities and velocities

The simulations we considered record the orbital elements of all particles at regular output intervals dt_out. Given any pair of particles, we computed their intrinsic collision probability and impact velocity using the Öpik-like algorithm described in Wetherill (1967), implemented in a fortran code by Farinella and Davis and kindly provided to us. The algorithm considers the orbital elements a,e,i (semi-major axis, eccentricity, and inclination) of each of the two particles and assumes that the orbital angles ω,Ω,M (argument of perihelion, longitude of node, mean anomaly) precess linearly with time (so that their values are random on a sufficiently long time interval), without inducing changes on (a,e,i). Clearly, these are approximations, but they are basically correct until the particles reach very high inclinations and undergo large-amplitude Kozai cycles (Vokrouhlický et al. 2012; Pokorný & Vokrouhlický 2013), which is not the case for the pre-instability disk or for most of the particles in the scattered disk (Kozai cycles for transplanetary orbits are pronounced only in mean-motion resonances or for planet-crossing orbits; Thomas & Morbidelli 1996). The use of a collisional probability algorithm like Wetherill’s on the output of a numerical integration is standard practice and leads to quite accurate results (e.g., Levison et al. 2000; Rickman et al. 2014).

The code returns the intrinsic collision probability P_i, which is the probability of a point-like projectile hitting a target with in an year. Thus, the probability that two objects of radii R₁ and R₂ (expressed in km) collide over a time interval δt is P_coll = P_i(R₁ + R₂)²δt. The impact velocity v_coll is the mean of the relative velocities between the two orbits over all collision configurations. It corresponds to the velocity of approach before any acceleration due to the mutual attraction between the two bodies. The latter is negligible for planetesimals.

For the pre-instability disk, it is not necessary to consider the collision probability of each of the simulated particles. Averaged values are enough. However, given the radial excitation gradient shown in the top panel of Fig. 1, we divided the disk into three zones: zone I with 15 <a< 20 AU, zone II with 20 <a< 25 AU, and zone III with 25 <a< 30 AU. Then, denoting by k the particles in one zone and m those in the other zone (possibly the same zones), we computed the mean intrinsic collision probability as $\mathit{P̅}\mathit{i}\mathrm{=}\frac{\mathrm{1}}{\mathit{KM}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}\sum _{\mathit{m}\mathrm{=}\mathrm{1}}^{\mathit{M}}{\mathit{P}}_{\mathit{i}}\mathrm{\left(}\mathit{k,m}\mathrm{\right)}\mathit{,}$(1)where K and M are the total numbers of particles in the considered disk zones and P_i(k,m) is the intrinsic probability between particles k and m. Similarly, the mean impact velocity (weighted by collision probability) is $\mathit{v̅}\mathrm{coll}\mathrm{=}\frac{\mathrm{1}}{\mathit{KM}\mathit{P̅}\mathit{i}}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}\sum _{\mathit{m}\mathrm{=}\mathrm{1}}^{\mathit{M}}{\mathit{P}}_{\mathit{i}}\mathrm{\left(}\mathit{k,m}\mathrm{\right)}{\mathit{v}}_{\mathrm{coll}}\mathrm{\left(}\mathit{k,m}\mathrm{\right)}\mathit{,}$(2)where v_coll(k,m) is the collision velocity between particles k and m.

For the simulation of the disk dispersal, we instead individually computed the collision probability for each particle that will be a dynamical survivor in the scattered disk at the end of the simulation, against all other particles. Denoting by j a scattered disk particle and by l any other particle, the mean collision probability of particle j at time t is $\mathit{P̅}\mathit{i}\mathrm{\left(}\mathit{j,t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{L}\mathrm{\left(}\mathit{t}\mathrm{\right)}}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{\mathit{L}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathit{P}}_{\mathit{i}}\mathrm{\left(}\mathit{j,l}\mathrm{\right)}\mathit{,}$(3)where L(t) is the total number of particles surviving in the integration at time t. For the velocity, we have $\mathit{v̅}\mathrm{coll}\mathrm{\left(}\mathit{j,t}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{L}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{P̅}\mathit{i}\mathrm{\left(}\mathit{j,t}\mathrm{\right)}}\sum _{\mathit{l}\mathrm{=}\mathrm{1}}^{\mathit{L}\mathrm{\left(}\mathit{t}\mathrm{\right)}}{\mathit{P}}_{\mathit{i}}\mathrm{\left(}\mathit{j,l}\mathrm{\right)}{\mathit{v}}_{\mathrm{coll}}\mathrm{\left(}\mathit{j,l}\mathrm{\right)}\mathit{.}$(4)

2.3. Size distributions and disruption probabilities

In the method introduced in CM03/07, an initial size distribution for the swarm of planetesimals represented by each simulation particle is defined. Using the pre-computed collision probabilities among the simulation particles, the number of collisions occurring between planetesimals of any given size is computed at each step. From the impact velocities and masses of projectile and target, the consequences of the collisions (cratering event, catastrophic break-up) and the size distribution of the generated fragments are computed. In this way, the evolution of the size distributions associated with each simulation particle is computed. In the end, the fraction of the planetesimals of a given size is evaluated that are survivors of the original population or collisional fragments of larger planetesimals.

This approach is correct, but it is too computationally expensive for our goal in this paper, which is just to assess whether a comet-sized object might have avoided catastrophic collisions. Moreover, it requires exploring a variety of initial size distributions, demanding a quite tedious exploration of the parameter space defining them.

We therefore modified and simplified the approach as described below. We used the approach typical of a mathematical demonstration ad absurdum (by reduction to the absurd). That is, we start by assuming that the planetesimals down to comet-sized objects are not significantly affected by collisions. This means that the planetesimal size distribution does not evolve with time and that the initial distribution has to be the same as the current distribution in the scattered disk, but scaled up by the inverse of the implantation efficiency (the fraction of the disk population surviving in the end in the scattered disk). We detail this size distribution below.

Then, using this distribution and the pre-computed collision probabilities and velocities, we evaluate the minimum size of a projectile that is capable of disrupting a comet-sized body and thus the number of catastrophic collisions n_coll that each comet-sized body might suffer (we detail this evaluation below). The probability that a comet-sized body has escaped all catastrophic collisions is then P_intact = exp(−n_coll). If P_intact is close to unity, then our assumption of a negligible collisional role is verified. But if P_intact is low, then we reach the absurd situation that by assuming that the planetesimal population was not affected by collisions, we conclude that most planetesimals should have been destroyed! This means that the assumption was incorrect, and hence, that the planetesimal size distribution was significantly affected by collisions.

With this approach we cannot compute the actual probability that a comet-like body has escaped all catastrophic collisions, but we know that it has to be lower than P_intact. In fact, any initial size distribution evolving by collisions toward the current scattered disk distribution must originally have had more bodies than we assumed (because of collisional comminution), and therefore the probability that a given body was catastrophically disrupted must be higher than we computed (because of a larger initial number of projectiles). If the value of P_intact that we computed is already low, this is enough for our purposes.

We here consider as comet-sized bodies objects with radius R = 2 km, which is appropriate for comet 67P/Churyumov-Gerasimenko, the target of the Rosetta mission (see Rickman et al. 2015).

2.3.1. Scattered disk size distribution and the minimal number of comet-sized objects in the original transplanetary disk

The most recent estimate of the scattered disk population has been presented in Brasser & Morbidelli (2013). In that work, as in Duncan & Levison (1997), the number of comet-sized bodies in the scattered disk was evaluated from (a) the number of known Jupiter-family comets in some given range of orbits and magnitudes for which the JFC sample is assumed complete, and (b) the numerical relationship between the scattered disk population and the Jupiter-family population that the former sustains, obtained from numerical simulations. With respect to previous estimates (e.g., Duncan & Levison 1997), the estimate in Brasser and Morbidelli is improved in two respects: it uses the most recent conversion from total magnitude to nuclear size from Fernández & Sosa (2012), and it is based on new simulations deriving Jupiter-family comets from a scattered disk that is excited in inclination (the original work by Duncan and Levison assumed that inclinations in the scattered disk are only of a few degrees, which has then been refuted by observations). Brasser and Morbidelli concluded that there are 2 × 10⁹ bodies in the scattered disk today that are larger than 2.3 km in diameter. Given a scattered disk implantation efficiency of 1%, this means that the original transplanetary disk, if not affected by collisional comminution, probably contained 2 × 10¹¹ of these bodies.

This is most likely a lower bound for the original disk population for three reasons. First, in a subsequent work accounting for an improved fading law (probability of a comet not surviving more than n perihelion passages), Brasser & Wang (2015) raised the estimate of the scattered disk population to 6 × 10⁹ bodies larger than 2.3 km in diameter. Second, serendipitous stellar occultation observations by the HST guiding sensors (Schlichting et al. 2009) and by the Corot survey (Liu et al. 2015) suggested that the average sky density of bodies larger than 250 m in radius over a ± 5° ecliptic band is 2 × 10⁷/deg². This means that there are at least 7 × 10¹⁰ bodies of this size in the trans-Neptunian region; with a cumulative size distribution proportional to R^-2, this implies 3.5 × 10⁹ objects with D> 2.3 km. It is unclear which population (cold EKB, hot EKB, scattered disk) the detected objects belong to. But, given that the scattered disk outnumbers the others (compare Trujillo et al. 2000 with Fraser et al. 2014 for the observational point of view and Brasser & Morbidelli 2013 with Nesvorný 2015 for the modeling point of view), the number above can be considered to be an estimate – if not a lower bound – of the scattered disk population. Third, repeating the same exercise for the Oort cloud population, Brasser & Morbidelli (2013) estimated that the primordial transplanetary disk probably contained 10¹² objects with D> 2.3 km; the reality therefore most likely lies in between 2 × 10¹¹ and 10¹².

Thus, to remain conservative (i.e., underestimate the total number of collisions), we assumed that the transplanetary disk contained 2 × 10¹¹ objects with D> 2.3 km. As for the size distribution, we again turned to comet observations. Estimates of the JFC size distribution vary significantly among authors, from quite steep (exponent of the differential distribution close to −3.5 – Fernández et al. 1999; Tancredi et al. 2006 – or even steeper – Belton 2015) to shallow (differential slope of −2.6; Lowry & Weissman 2003). Consequently, we assumed for the nominal differential slope the value −3 (Meech et al. 2004; Lamy et al. 2004; Snodgrass et al. 2011), but we also studied the dependence of the results on exponents for the differential distribution ranging from −2.5 to −3.5.

A shallow size distribution is preferred according the the most recent planetesimal formation models (Johansen et al. 2015: q = −2.8). TNO surveys (e.g. Fraser et al. 2014) also suggest that the size distribution of objects smaller than 50 km in radius is shallow (q between −3.1 and −2.5, although it may steepen up for not yet detectable comet-size bodies).

For reference, a disk with a size distribution similar to that of the hot EKB (Fraser et al. 2014), that is, with a differential slope of −3 for R< 50 km and −5 for R> 50 km, a total number of 2 × 10¹¹ objects with R> 1.15 km and a density of 1 g/cm³ would have a total mass of 35 Earth masses, in good agreement with the mass required by the Nice model (Gomes et al. 2005; Morbidelli et al. 2007; Nesvorný & Morbidelli 2012).

2.3.2. Minimum size of catastrophic projectiles

The kinetic energy of an impact that can catastrophically destroy an object is ${\mathit{E}}_{\mathrm{disrupt}}\mathrm{=}\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi \rho }{\mathit{R}}^{\mathrm{3}}{\mathit{Q}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathit{,}$(5)where ρ is the bulk density of the target of radius R; Q^∗(R) is the specific energy for disruption and is size dependent. There are several Q^∗(R) laws proposed in the literature for various materials. Benz & Asphaug (1999) proposed two such laws for bodies made of “strong ice”, which is hit at 1 km s^-1 and 3 km s^-1, respectively. As we show in Sect. 3, the former velocity is well adapted to the pre-instability disk, while the second is suitable for the disk dispersal phase. For a R = 2 km body the two values of Q^∗(R) are similar. Leinhardt & Stewart (2009) produced a Q^∗ law for bodies made of “weak ice”, which is hit at 1 km s^-1. Their Q^∗(R) function follows the general trend of those in Benz & Asphaug (1999), but the value for R = 2 km is about an order of magnitude lower (see their Fig. 11). The scaling of Q^∗ with velocity given in Eq. (2) of Stewart & Leinhardt (2009) gives a value 2.4 times higher for a velocity of 3 km s^-1, which is still four times lower than that reported in Benz & Asphaug (1999) for the same speed. Even if the Leinhardt and Stewart value may be more appropriate for pristine, low-density planetesimals, we used the values reported by Benz and Asphaug values to be conservative once again. This probably overestimates the minimal size of a projectile that is capable of disrupting the target and thus underestimates the number of catastrophic collisions.

When E_disrupt is known, the minimal size of a catastrophic projectile r_p is given by the equation $\frac{\mathrm{4}}{\mathrm{3}}\mathit{\pi \rho }{\mathit{r}}_{\mathrm{p}}^{\mathrm{3}}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{v}}_{\mathrm{coll}}^{\mathrm{2}}\mathrm{=}{\mathit{E}}_{\mathrm{disrupt}}\mathit{,}$(6)which means that the higher E_disrupt, the larger is r_p. If one assumes that the bulk density of projectile and target is the same, ρ simplifies from the right-hand and left-hand sides of Eq. (6) and the result is independent of ρ.

2.3.3. Total number of catastrophic events

When the minimum size of a catastrophic projectile is known, the total number of catastrophic impacts for a target of radius R_T is computed as ${\mathit{N}}_{\mathrm{coll}}\mathrm{=}\mathrm{\left(}\mathit{P̅}\mathit{i}\mathit{\delta t}\mathrm{\right)}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{p}}}^{{\mathit{R}}_{\max }}\mathrm{(}{\mathit{R}}_{\mathrm{T}}\mathrm{+}{\mathit{R}}_{\mathrm{p}}{\mathrm{)}}^{\mathrm{2}}\mathit{N}\mathrm{(}{\mathit{R}}_{\mathrm{p}}\mathrm{)}\mathrm{d}{\mathit{R}}_{\mathrm{p}}\mathit{,}$(7)where is the considered intrinsic probability (averaged over the ensemble of potential projectiles, as explained in Sect. 2.2), δt is the considered time-span, N(R_p)dR_p is the differential size distribution, r_p is the minimum size for a catastrophic projectile, and R_max is the maximum size for which the considered size distribution is valid. Given that the size distribution of the trans-Neptunian populations changes from steep (at the large-size end) to shallow (at the small-size end) at a size of approximately R ~ 50 km (Bernstein et al. 2004; Fuentes et al. 2009; Fraser et al. 2014), we assumed R_max = 50 km. We neglected the relatively small contribution by projectiles of even larger sizes.

Equation (7) is often approximated by ${\mathit{N}}_{\mathrm{coll}}\mathrm{=}\mathit{P̅}\mathrm{i}\mathit{\delta t}{\mathit{R}}_{\mathrm{T}}^{\mathrm{2}}{\mathrm{\int }}_{{\mathit{r}}_{\mathrm{p}}}^{{\mathit{R}}_{\max }}\mathit{N}\mathrm{\left(}{\mathit{R}}_{\mathrm{p}}\mathrm{\right)}\mathrm{d}{\mathit{R}}_{\mathrm{p}}\mathrm{=}\mathit{P̅}\mathrm{i}\mathit{\delta t}{\mathit{R}}_{\mathrm{T}}^{\mathrm{2}}\mathit{N}\mathrm{(}\mathit{>}{\mathit{r}}_{\mathrm{p}}\mathrm{)}\mathit{,}$(8)where N( >r_p) is the cumulative number of bodies larger than r_p. This approximation is good for steep size distributions or in the limit r_p → 0. However, for shallow size distributions like the one we consider here and r_p not much smaller than R_T (as is the case for low-velocity collisions), the approximation is not very precise. Thus, we solved the integral Eq. (7) exactly. That is, denoting by q the exponent of the differential size distribution, the primitive of the integrand in (7) is $\begin{array}{c}\\ \mathit{q}& \mathrm{=}& \mathrm{-}\mathrm{2.5}\mathrm{:} \mathrm{-}\mathrm{2}\mathrm{(}{\mathit{R}}_{\mathrm{T}}^{\mathrm{2}}\mathrm{+}\mathrm{6}{\mathit{R}}_{\mathrm{T}}{\mathit{R}}_{\mathrm{p}}\mathrm{-}\mathrm{3}{{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{/}\mathrm{(}\mathrm{3}{{\mathit{R}}_{\mathrm{p}}^{\mathrm{3}\mathit{/}\mathrm{2}}}^{\mathrm{)}}\\ \mathit{q}& \mathrm{=}& \mathrm{-}\mathrm{3}\mathrm{:} \mathrm{-}{\mathit{R}}_{\mathrm{T}}^{\mathrm{2}}\mathit{/}\mathrm{(}\mathrm{2}{{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{-}\mathrm{2}{\mathit{R}}_{\mathrm{T}}\mathit{/}{\mathit{R}}_{\mathrm{p}}\mathrm{+}\log \mathrm{(}{\mathit{R}}_{\mathrm{p}}\mathrm{)}\\ \mathit{q}& \mathrm{=}& \mathrm{-}\mathrm{3.5}\mathrm{:} \mathrm{-}\mathrm{2}\mathrm{(}\mathrm{3}{\mathit{R}}_{\mathrm{T}}^{\mathrm{2}}\mathrm{+}\mathrm{10}{\mathit{R}}_{\mathrm{T}}{\mathit{R}}_{\mathrm{p}}\mathrm{+}\mathrm{15}{{\mathit{R}}_{\mathrm{p}}^{\mathrm{2}}}^{\mathrm{)}}\mathit{/}\mathrm{(}\mathrm{15}{{\mathit{R}}_{\mathrm{p}}^{\mathrm{5}\mathit{/}\mathrm{2}}}^{\mathrm{)}}\mathit{.}\end{array}$(9)

3. Results

We report here the results obtained for the pre-instability disk and the disk dispersal phase, obtained by applying the methods described in the previous section.

3.1. Disruptive collisions in the pre-instability disk

We show in Table 1 the results for the intrinsic collision probability , the collision velocity v_coll, and the minimum size of a catastrophic projectile r_p for a target with R = 2 km, considering all possible combinations of disk zones for projectile and target.

As seen, the intrinsic collision probability is higher if both target and projectile are in the inner part of the disk than it is in the outer part. The collision velocity is also higher. The mean intrinsic collision probability is lower if target and projectile belong to different disk zones, because not all particle orbits from the two zones intersect.

Table 2 reports the number of catastrophic collisions expected for a 2 km target for the three considered values of the exponent q of the differential size distribution. The calculation was made assuming δt = 400 My (the expected lifetime of the pre-instability phase), and applying Eqs. (7) and (9) to the numbers reported in Table 1 for projectiles in each disk zone. Because the surface density of the disk is assumed to be proportional to 1 /r in the simulation by Levison et al. (2011) that we used, an equal number of projectiles of a given size was assumed to initially exist in each of the disk zones.

We note that the total number of collisions for comets in zone III of the disk has a very weak dependence on q because the size of the minimum catastrophic projectile is quite large (0.5 km in radius). The opposite is true for comet-sized targets in zone I of the disk. Clearly, in all cases the total number of collisions is larger than 1. This means that the probability of a 2 km body escaping from all catastrophic collisions is low. According to the numbers in Table 2, this probability is always lower than exp(−12) = 6 × 10^-6; using the numbers in parentheses, we derive exp(−7.9) = 4 × 10^-4.

Table 3 reports the radius of the comets for which the probability of a catastrophic impact over the lifetime of the disk is 50%, again for each disk zone and assumed value of q. This size is extremely large, in most cases exceeding 50 km. This is because the number of catastrophic impacts only weakly depends on the size of the target. If one assumes for simplicity that Q^∗ is independent of size, the size of the minimal catastrophic projectile scales linearly with the target size R_T; then, the number of catastrophic impacts, using Eq. (8), scales as ${\mathit{R}}_{\mathrm{T}}^{\mathrm{2}\mathrm{+}\mathit{q}\mathrm{+}\mathrm{1}}$, which for q = −3 eliminates the dependence on R_T. For the compilation of Table 3 we nevertheless used the dependence of Q^∗ on radius given in Benz & Asphaug (1999) and the non-approximated formulae (9).

This result is valid for both scattered disk comets (JFCs) and Oort cloud comets (LPCs/HTCs) because both reservoirs form from the same disk (Brasser & Morbidelli 2013). Thus we conclude that if the giant planet instability occurred as late as 4.1 Gy ago, the possibility of a 2 km comet to be a pristine planetesimal and not a collisional fragment is very slim.

3.2. Disruptive collisions during disk dispersal

As described in Sect. 2.2, for the disk dispersal phase we considered each individual particle ending in the scattered disk because their orbital histories can be very diverse. However, because the assumption in our ad absurdum approach is that the disk is not collisionally active and its size distribution does not evolve, for each target particle j we can average over time the value of given in Eq. (3) as$\mathit{P̅}\mathrm{i}\mathrm{\left(}\mathit{j}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathit{L}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\sum _{\mathit{t}}\mathit{L}\mathrm{\left(}\mathit{t}\mathrm{\right)}\mathit{P̅}\mathrm{i}\mathrm{\left(}\mathit{j,t}\mathrm{\right)}$(10)and apply the result over a time interval δt = 350 My, which is the integration time-span.

However, this 350 My simulation time-span covers only a small fraction of the lifetime of the solar system, and in principle, there may still be a significant collisional evolution in the scattered disk over the remaining ~4 Gy of solar system history. To estimate the collision probability over this remaining time, we proceeded as follows. We assumed that the orbital distribution in the scattered disk does not evolve with time and remains equal to that shown in the bottom panel of Fig. 1, but we assumed that the scattered disk population decays with time. The scattered disk at the end of the simulation accounts for 5% of the initial disk particles, and for the rest of the solar system lifetime it would decay to about 1% (Brasser & Morbidelli 2013). We therefore assumed that number of scattered disk particles decays as exp(−t/τ) with τ such that after 4 Gy the population is reduced by a factor of 5. The integral of the exponential function over 4 Gy with such a value of τ is 0.5. Thus, we took the last computed value of for each scattered disk particle (t = 350 My) and multiplied it by δt = 4 Gy and then divided by 2. The integrated collision probability for the last 4 Gy is a fraction (typically from 10% to 80%) of that integrated for the first 350 My. This is because the disk is much more populated at the beginning and, moreover, the most collisionally active phase is when the disk is just stirred up by the planets’ action (Stern & Weissman 2001). Instead, once the scattered disk is formed, the collision probability per unit time per particle is strongly reduced as a result of the large orbital space that the scattered disk fills and the long orbital periods.

Given the typical collision velocities of 2−4 km s^-1, the typical size of the smallest catastrophic projectile for a target of R = 2 km is ~100 m. One may wonder whether bodies this small existed in the disk. The occultation observations mentioned above (Schlichting et al. 2009; Liu et al. 2015) show that bodies of this size exist today. Following our ad absurdum approach, we then need to assume that they existed in the original disk, because if the disk did not evolve collisionally, no objects could be generated.

Fig. 2 Number of expected catastrophic collisions for each particle surviving in the scattered disk at the end of the disk dispersal simulation. The symbols depict different values for the exponent of the differential size distribution q, as labeled in the plot.

Figure 2 shows for each scattered disk particle the expected number of catastrophic collisions that a 2 km target should have suffered; the colors refer to different values of q. The number of catastrophic collisions changes considerably from particle to particle because the dynamical histories of scattered disk objects can be very diverse. We see that if q = −3.5 or steeper, clearly each comet-sized object should have suffered at least two catastrophic collisions, with an average of 4.7 collisions. The radius of comets for which the average number of catastrophic impacts would be 0.5 is approximately 10 km. This excludes the possibility, proposed by Belton (2015), of the size distribution of the scattered disk being steeper than q = −3.5 below this radius. In fact, the disk would be collisionally evolved and therefore would have acquired a collisional equilibrium size distribution, which implies q = −3.5 (Dohnanyi 1969) or shallower (| q | < 3.5; O’Brien & Greenberg 2003).

If the distribution is very shallow (q = −2.5), Fig. 2 shows that about half of the comets should have had no catastrophic collisions, the average number of catastrophic collisions per object being 0.5. The nominal case q = −3.0 is borderline. Most comets should have had at least one catastrophic collision (the average being 1.5 collisions per comet) but some, with favorable orbital histories, would have had no collision at all. Had we used the four times lower value of Q^∗ from Leinhardt and Stewart (see Sect. 2.3.2), the number of catastrophic impacts would have increased by a factor of ~3 for q = −3.5, a factor ~2.5 for q = −3 and a factor of 2 for q = −2.5.

Figure 3 shows the same results, but using a different representation. The number of collisions N_coll is converted into a probability of avoiding all collisions P(0) = exp(−N_coll). Then, the normalized cumulative distribution of the P(0) values is plotted. The thick curves represent the nominal Q^∗ value from Benz & Asphaug (1999) and the thin curves a Q^∗ value four times lower.

Fig. 3 Fraction of particles ending in the scattered disk with a probability of escaping all catastrophic collisions P(0) lower than that indicated on the horizontal axis. This is an alternative representation of the results shown in Fig. 2. The different line styles refer to different exponents for the differential size distribution q, as labeled in the plot. The thick curves correspond to the value of Q∗ given in Benz & Asphaug (1999) and the thin curves to a Q∗ value four times lower, as in Leinhardt & Stewart (2009) with the velocity scaling provided in Stewart & Leinhardt (2009).

If the disk dispersal occurred late, the number of catastrophic collisions found in this section should be added to those reported in Table 2. This would mean that all original comets should have been destroyed. If instead the dispersal of the disk occurred soon after the disappearance of gas from the disk, only the results illustrated in Figs. 2 and 3 apply. In this case, if the original and current size distributions in the comet-sized range are shallow and a quite high value of Q^∗ applies, there is the possibility of a fraction of the comets being pristine objects that escaped catastrophic collisions.

4. Conclusions

We have estimated the total number of catastrophic collisions that a typical Jupiter-family comet (here assumed to have a radius of R = 2 km) should have had during its dynamical lifetime, first in the transplanetary disk, and then in the scattered disk.

We have shown that if the transplanetary disk beyond the original orbit of Neptune has been dispersed by a late dynamical instability of the giant planets occurring ~4.1 Gy ago, comet-sized objects should have suffered numerous catastrophic collisions in the pre-instability phase. Thus, not only JFCs, but also Oort cloud comets should be fragments of originally larger bodies. Because the late instability of the giant planet system is, at the current level of understanding, the best explanation for the trigger of the Late Heavy Bombardment of the solar system, the formation of late lunar basins on the Moon (Bottke et al. 2012; Morbidelli et al. 2012) and the impact age record on meteorites (Marchi et al. 2013), this is probably the conclusion of our work.

However, in the hypothetical case that the dispersal of the disk occurred early, the collisional evolution of comet-sized bodies ending in the scattered disk would have been less severe. If the size distribution of comet-sized objects in today’s scattered disk and in the primordial transplanetary disk was shallow (differential index | q | ≲ 3), it is possible in principle that a significant fraction of comet-sized objects escaped all catastrophic collisions.

We recall, however, that throughout our study we have taken the most conservative assumptions, so that the number of catastrophic collisions that we computed should be considered as a lower estimate. In fact, we have considered an initial size distribution in the disk that contains the minimum possible number of comet-sized objects (Sect. 2.3.1). In addition, we assumed the specific energy for catastrophic disruption given in Benz & Asphaug (1999), which probably overestimates the energy appropriate for weak icy aggregates (Leinhardt & Stewart 2009; Sect. 2.3.2): the adoption of the four times lower specific energy for disruption of Leinhardt & Stewart (2009) would have multiplied the number of catastrophic impacts shown in Fig. 2 by a factor 2.5 for q = −3 and a factor of 2 if q = −2.5, while producing the thin curves in Fig. 3. Moreover, the ad absurdum approach that we followed by definiton provides just a minimal estimate of the number of collisions (Sect. 2.3). Finally, for each catastrophic collision, the number of quasi-catastrophic collisions would be much higher (being caused by smaller projectiles, which are more numerous). Thus, even if a comet had not suffered, by chance, any catastrophic collision, its morphology would have been sculpted by numerous large subcatastrophic impacts.

Therefore, we conclude that typical JFCs of the size of 67P/ Churyumov-Gerasimenko most likely are not intact planetesimals, but are either fragments of originally larger bodies (the most likely case), or are planetesimals strongly sculpted by barely subcatastrophic impacts². In the latter case, this sculpting may disappear after the comets have been severely eroded by sublimation or splitting. However, some JFCs are statistically likely to bear these scars.

There are some properties of comets that appear to contradict a collisional origin and instead suggest that comets are primordial survivors. In the Introduction we mentioned the low bulk densities and negligible tensile strengths along with the large abundance of supervolatiles like CO. These properties are certainly compatible with an origin of comets as primordial survivors of the icy planetesimal population. Hence, our results lead to the obvious question of whether collisions are able to conserve these primordial properties in the fragments they produce. This question is currently unanswered and merits careful consideration.

If a detailed modeling of the collisions between icy planetesimals would show that the primordial-like features of comets are not preserved in the fragments, one may suspect that our current vision of outer solar system evolution is not appropriate. For instance, there might have been no delay in the dynamical instability, or the disk remained less self-excited due to a smaller number of large bodies than we envision, or there was a drastic cut-off in the size distribution affecting sub-km objects, thus limiting the number of catastrophic projectiles.

Acknowledgments

A.M. was supported by ANR, pro ject number ANR-13–13-BS05-0003-01 projet MOJO(Modeling the Origin of JOvian planets). H.R. was supported by Grant No. 2011/01/B/ST9/05442 of the Polish National Science Center. The authors thank the anonymous referee for constructive comments.

References

All Tables

All Figures

Fig. 1 Top: semi-major axis vs. eccentricity distribution of the transplanetary disk under the stirring effect of an embedded population of 1000 Pluto-mass bodies according to Levison et al. (2011). This snapshot of the distribution is taken after 300 My of evolution. Bottom: semi-major axis vs. eccentricity distribution of the scattered disk produced by the dispersal of the transplanetary disk by the giant planet instability according to Gomes et al. (2005). Here the snapshot of the scattered disk orbital distribution is taken 350 My after the beginning of the planet instability, when the scattered disk contains 5% of the original disk’s particles.

In the text

	Fig. 2 Number of expected catastrophic collisions for each particle surviving in the scattered disk at the end of the disk dispersal simulation. The symbols depict different values for the exponent of the differential size distribution q, as labeled in the plot.
In the text

Fig. 3 Fraction of particles ending in the scattered disk with a probability of escaping all catastrophic collisions P(0) lower than that indicated on the horizontal axis. This is an alternative representation of the results shown in Fig. 2. The different line styles refer to different exponents for the differential size distribution q, as labeled in the plot. The thick curves correspond to the value of Q∗ given in Benz & Asphaug (1999) and the thin curves to a Q∗ value four times lower, as in Leinhardt & Stewart (2009) with the velocity scaling provided in Stewart & Leinhardt (2009).

In the text