Comets as collisional fragments of a primordial planetesimal disk
^{1} Département Lagrange, University of Nice–Sophia Antipolis, CNRS, Observatoire de la Côte d’Azur, 06304 Nice, France
email: morby@obsnice.fr
^{2} P.A.S. Space Research Center, Bartycka 18A, 00716 Warszawa, Poland
^{3} Dept. of Physics and Astronomy, Uppsala University, Box 516, 75120 Uppsala, Sweden
Received: 17 March 2015
Accepted: 16 April 2015
Context. The Rosetta mission and its exquisite measurements have revived the debate on whether comets are pristine planetesimals or collisionally evolved objects.
Aims. We investigate the collisional evolution experienced by the precursors of current comet nuclei during the early stages of the solar system in the context of the socalled Nice model.
Methods. We considered two environments for the collisional evolution: (1) the transplanetary planetesimal disk, from the time of gas removal until the disk was dispersed by the migration of the ice giants; and (2) the dispersing disk during the time that the scattered disk was formed. We performed simulations using different methods in the two cases to determine the number of destructive collisions typically experienced by a comet nucleus of 2 km radius.
Results. In the widely accepted scenario, where the dispersal of the planetesimal disk occurred at the time of the Late Heavy Bombardment about 4 Gy ago, cometsized planetesimals have a very low probability of surviving destructive collisions in the disk. On the extreme assumption that the disk was dispersed directly upon gas removal, a significant fraction of the planetesimals might have remained intact. However, these survivors would still bear the marks of many nondestructive impacts.
Conclusions. The Nice model of solar system evolution predicts that typical kmsized comet nuclei are predominantly fragments resulting from collisions experienced by larger parent bodies. An important goal for future research is to investigate whether the observed properties of comet nuclei are compatible with such a collisional origin.
Key words: comets: general / Kuiper belt: general / protoplanetary disks / Oort Cloud
© 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/ShoemakerLevy 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 lowvelocity 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 supervolatile CO molecule to their outgassing activity (e.g., BockeléeMorvan 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 EdgeworthKuiper 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 smallscale 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 particleinabox models used before. CM03/07 showed that for some appropriate initial size distributions (similar to that currently observed in the transNeptunian region), a sufficiently large number of cometsized 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 transplanetary 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 clearcut 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, longperiod comets (LPCs) and Halleytype 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 socalled 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 meanmotion resonances with Neptune), is lower than 10^{3} (see Nesvorný 2015, for the most updated estimate); this means that about 1000 Plutosized 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 cometsized 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 planetcrossing 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 transplanetary 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 cometsized 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 timestep 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 preexisting simulations that represent well the two phases of the evolution of the transplanetary disk described in the Introduction: the preinstability phase and the dynamical dispersal phase.
2.1.1. Preinstability disk
For the preinstability 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 selfexcitation it would suffer if it contained 1000 Plutomass 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 selfexcitation 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 selfstirring process increases the orbital excitation as 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 preinstability 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 semimajor 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)^{1}.
Fig. 1 Top: semimajor axis vs. eccentricity distribution of the transplanetary disk under the stirring effect of an embedded population of 1000 Plutomass bodies according to Levison et al. (2011). This snapshot of the distribution is taken after 300 My of evolution. Bottom: semimajor 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. 

Open with DEXTER 
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 preinstability 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 preinstability 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 timespan 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 Öpiklike 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 (semimajor 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 largeamplitude Kozai cycles (Vokrouhlický et al. 2012; Pokorný & Vokrouhlický 2013), which is not the case for the preinstability disk or for most of the particles in the scattered disk (Kozai cycles for transplanetary orbits are pronounced only in meanmotion resonances or for planetcrossing 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 pointlike projectile hitting a target with in an year. Thus, the probability that two objects of radii R_{1} and R_{2} (expressed in km) collide over a time interval δt is P_{coll} = P_{i}(R_{1} + R_{2})^{2}δ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 preinstability 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 (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 (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 (3)where L(t) is the total number of particles surviving in the integration at time t. For the velocity, we have (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 precomputed 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 breakup) 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 cometsized 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 cometsized 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 precomputed collision probabilities and velocities, we evaluate the minimum size of a projectile that is capable of disrupting a cometsized body and thus the number of catastrophic collisions n_{coll} that each cometsized body might suffer (we detail this evaluation below). The probability that a cometsized 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 cometlike 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 cometsized bodies objects with radius R = 2 km, which is appropriate for comet 67P/ChuryumovGerasimenko, the target of the Rosetta mission (see Rickman et al. 2015).
2.3.1. Scattered disk size distribution and the minimal number of cometsized 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 cometsized bodies in the scattered disk was evaluated from (a) the number of known Jupiterfamily 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 Jupiterfamily 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 Jupiterfamily 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^{9} 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^{11} 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^{9} 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^{7}/deg^{2}. This means that there are at least 7 × 10^{10} bodies of this size in the transNeptunian region; with a cumulative size distribution proportional to R^{2}, this implies 3.5 × 10^{9} 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^{12} objects with D> 2.3 km; the reality therefore most likely lies in between 2 × 10^{11} and 10^{12}.
Thus, to remain conservative (i.e., underestimate the total number of collisions), we assumed that the transplanetary disk contained 2 × 10^{11} 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 cometsize 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^{11} objects with R> 1.15 km and a density of 1 g/cm^{3} 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 (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 preinstability 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, lowdensity 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 (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 righthand and lefthand 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 (7)where is the considered intrinsic probability (averaged over the ensemble of potential projectiles, as explained in Sect. 2.2), δt is the considered timespan, 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 transNeptunian populations changes from steep (at the largesize end) to shallow (at the smallsize 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 (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 lowvelocity 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 (9)
Table of results for the preinstability disk.
3. Results
We report here the results obtained for the preinstability disk and the disk dispersal phase, obtained by applying the methods described in the previous section.
3.1. Disruptive collisions in the preinstability 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 preinstability 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.
Number of disruptive collisions expected for a target of R = 2 km located in each disk zone as a function of the exponent q of the differential size distribution.
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 cometsized 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 , 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 nonapproximated 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(10)and apply the result over a time interval δt = 350 My, which is the integration timespan.
However, this 350 My simulation timespan 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. 

Open with DEXTER 
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 cometsized 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). 

Open with DEXTER 
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 cometsized 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 Jupiterfamily 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, cometsized objects should have suffered numerous catastrophic collisions in the preinstability 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 cometsized bodies ending in the scattered disk would have been less severe. If the size distribution of cometsized 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 cometsized 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 cometsized 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 quasicatastrophic 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/ ChuryumovGerasimenko 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^{2}. 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 primordiallike 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 selfexcited due to a smaller number of large bodies than we envision, or there was a drastic cutoff in the size distribution affecting subkm objects, thus limiting the number of catastrophic projectiles.
The special case of comet 67P and constraints on its origin has been treated by Rickman et al. (2015).
Acknowledgments
A.M. was supported by ANR, pro ject number ANR13–13BS05000301 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
 A’Hearn, M. F., Feaga, L. M., Keller, H. U., et al. 2012, ApJ, 758, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Altwegg, K., Balsiger, H., BarNun, A., et al. 2015, Science, 347, 1261952 [CrossRef] [Google Scholar]
 Asphaug, E., & Benz, W. 1996, Icarus, 121, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Belton, M. J. S. 2015, Icarus, 245, 87 [NASA ADS] [CrossRef] [Google Scholar]
 Benz, W., & Asphaug, E. 1999, Icarus, 142, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Bernstein, G. M., Trilling, D. E., Allen, R. L., et al. 2004, AJ, 128, 1364 [NASA ADS] [CrossRef] [Google Scholar]
 Bitsch, B., Johansen, A., Lambrechts, M., & Morbidelli, A. 2015, A&A, 575, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Blum, J., Gundlach, B., Mühle, S., & TrigoRodriguez, J. M. 2014, Icarus, 235, 156 [NASA ADS] [CrossRef] [Google Scholar]
 BockeléeMorvan, D., Crovisier, J., Mumma, M. J., & Weaver, H. F. 2004, Comets II, 391 [Google Scholar]
 Bottke, W. F., Levison, H. F., Nesvorný, D., & Dones, L. 2007, Icarus, 190, 203 [NASA ADS] [CrossRef] [Google Scholar]
 Bottke, W. F., Vokrouhlický, D., Minton, D., et al. 2012, Nature, 485, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Brasser, R., & Morbidelli, A. 2011, A&A, 535, A41 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brasser, R., & Morbidelli, A. 2013, Icarus, 225, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Brasser, R., & Wang, J.H. 2015, A&A, 573, A102 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Chambers, J. E. 2007, Icarus, 189, 386 [NASA ADS] [CrossRef] [Google Scholar]
 Charnoz, S., & Morbidelli, A. 2003, Icarus, 166, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Charnoz, S., & Morbidelli, A. 2007, Icarus, 188, 468 [NASA ADS] [CrossRef] [Google Scholar]
 Ćuk, M. 2012, Icarus, 218, 69 [NASA ADS] [CrossRef] [Google Scholar]
 Ćuk, M., & Gladman, B. J. 2006, Icarus, 183, 362 [NASA ADS] [CrossRef] [Google Scholar]
 Ćuk, M., & Gladman, B. J. 2009, Icarus, 199, 237 [NASA ADS] [CrossRef] [Google Scholar]
 Davidsson, B. J. R., Gutiérrez, P. J., & Rickman, H. 2007, Icarus, 187, 306 [NASA ADS] [CrossRef] [Google Scholar]
 Davis, D. R., & Farinella, P. 1997, Icarus, 125, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531 [NASA ADS] [CrossRef] [Google Scholar]
 Duncan, M. J., & Levison, H. F. 1997, Science, 276, 1670 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Fernández, J. A., & Sosa, A. 2012, MNRAS, 423, 1674 [NASA ADS] [CrossRef] [Google Scholar]
 Fernández, J. A., Tancredi, G., Rickman, H., & Licandro, J. 1999, A&A, 352, 327 [NASA ADS] [Google Scholar]
 Fraser, W. C., Brown, M. E., Morbidelli, A., Parker, A., & Batygin, K. 2014, ApJ, 782, 100 [NASA ADS] [CrossRef] [Google Scholar]
 Fuentes, C. I., George, M. R., & Holman, M. J. 2009, ApJ, 696, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lamy, P. L., Toth, I., Fernandez, Y. R., & Weaver, H. A. 2004, Comets II, 223 [Google Scholar]
 Leinhardt, Z. M., & Stewart, S. T. 2009, Icarus, 199, 542 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., Duncan, M. J., Zahnle, K., Holman, M., & Dones, L. 2000, Icarus, 143, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Levison, H. F., Morbidelli, A., Tsiganis, K., Nesvorný, D., & Gomes, R. 2011, AJ, 142, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, C.Y., Doressoundiram, A., Roques, F., et al. 2015, MNRAS, 446, 932 [NASA ADS] [CrossRef] [Google Scholar]
 Lowry, S. C., & Weissman, P. R. 2003, Icarus, 164, 492 [NASA ADS] [CrossRef] [Google Scholar]
 Luu, J., Marsden, B. G., Jewitt, D., et al. 1997, Nature, 387, 573 [NASA ADS] [CrossRef] [Google Scholar]
 Marchi, S., Bottke, W. F., Cohen, B. A., et al. 2013, Nature Geoscience, 6, 303 [NASA ADS] [CrossRef] [Google Scholar]
 Meech, K. J., Hainaut, O. R., & Marsden, B. G. 2004, Icarus, 170, 463 [NASA ADS] [CrossRef] [Google Scholar]
 Minton, D. A., Jackson, A. P., Asphaug, E., Fassett, C. I., & Richardson, J. E. 2015, Early Solar System Bombardment III, abs. #3033 [Google Scholar]
 Morbidelli, A., Tsiganis, K., Crida, A., Levison, H. F., & Gomes, R. 2007, AJ, 134, 1790 [NASA ADS] [CrossRef] [Google Scholar]
 Morbidelli, A., Marchi, S., Bottke, W. F., & Kring, D. A. 2012, Earth Planet. Sci. Lett., 355, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Nesvorný, D. 2015, AJ, in press [arXiv:1504.06021] [Google Scholar]
 Nesvorný, D., & Morbidelli, A. 2012, AJ, 144, 117 [NASA ADS] [CrossRef] [Google Scholar]
 O’Brien, D. P., & Greenberg, R. 2003, Icarus, 164, 334 [NASA ADS] [CrossRef] [Google Scholar]
 Pokorný, P., & Vokrouhlický, D. 2013, Icarus, 226, 682 [NASA ADS] [CrossRef] [Google Scholar]
 Rickman, H. 1989, Adv. Space Res., 9, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Rickman, H. 2004, Comets II, 205 [Google Scholar]
 Rickman, H., Wiśniowski, T., Wajer, P., Gabryszewski, R., & Valsecchi, G. B. 2014, A&A, 569, A47 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rickman, H., Marchi, S., A’Hearn, M. F., et al. 2015, A&A, 583, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Robuchon, G., Nimmo, F., Roberts, J., & Kirchoff, M. 2011, Icarus, 214, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Safronov, V. S. 1977, IAU Colloq. 39: Comets, Asteroids, Meteorites: Interrelations, Evolution and Origins, 483 [Google Scholar]
 Schlichting, H. E., Ofek, E. O., Wenz, M., et al. 2009, Nature, 462, 895 [NASA ADS] [CrossRef] [Google Scholar]
 Schlichting, H. E., Fuentes, C. I., & Trilling, D. E. 2013, AJ, 146, 36 [NASA ADS] [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 [NASA ADS] [Google Scholar]
 Snodgrass, C., Fitzsimmons, A., Lowry, S. C., & Weissman, P. 2011, MNRAS, 414, 458 [NASA ADS] [CrossRef] [Google Scholar]
 Stern, S. A. 1991, Icarus, 90, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Stern, S. A., & Weissman, P. R. 2001, Nature, 409, 589 [NASA ADS] [CrossRef] [Google Scholar]
 Stewart, S. T., & Leinhardt, Z. M. 2009, ApJ, 691, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Tancredi, G., Fernández, J. A., Rickman, H., & Licandro, J. 2006, Icarus, 182, 527 [NASA ADS] [CrossRef] [Google Scholar]
 Thomas, F., & Morbidelli, A. 1996, Celest. Mech. Dyn. Astron., 64, 209 [NASA ADS] [CrossRef] [Google Scholar]
 Trujillo, C. A., Jewitt, D. C., & Luu, J. X. 2000, ApJ, 529, L103 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Vokrouhlický, D., Pokorný, P., & Nesvorný, D. 2012, Icarus, 219, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Wetherill, G. W. 1967,J. Geophys. Res., 72, 2429 [Google Scholar]
All Tables
Number of disruptive collisions expected for a target of R = 2 km located in each disk zone as a function of the exponent q of the differential size distribution.
All Figures
Fig. 1 Top: semimajor axis vs. eccentricity distribution of the transplanetary disk under the stirring effect of an embedded population of 1000 Plutomass bodies according to Levison et al. (2011). This snapshot of the distribution is taken after 300 My of evolution. Bottom: semimajor 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. 

Open with DEXTER  
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. 

Open with DEXTER  
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). 

Open with DEXTER  
In the text 