Formation of bilobed shapes by subcatastrophic collisions
A late origin of comet 67P’s structure
Physics Institute, University of Bern, NCCR PlanetS, Sidlerstrasse 5, 3012 Bern, Switzerland
email: martin.jutzi@space.unibe.ch; willy.benz@space.unibe.ch
Received: 19 May 2016
Accepted: 9 November 2016
Context. The origin of the particular shape of comet 67P/ChuryumovGerasimenko (67P) is a topic of active research. How and when it acquired its peculiar characteristics has distinct implications on the origin of the solar system and its dynamics.
Aims. We investigate how shapes such as that of comet 67P can result from a new type of lowenergy, subcatastrophic impact involving elongated, rotating bodies. We focus on parameters potentially leading to bilobed structures. We also estimate the probability of such structures surviving subsequent impacts.
Methods. We used a smooth particle hydrodynamics (SPH) shock physics code to model the impacts, the subsequent reaccumulation of material and the reconfiguration into a stable final shape. The energy increase as well as the degree of compaction of the resulting bodies were tracked in the simulations.
Results. Our modelling results suggest that the formation of bilobed structures like 67P is a natural outcome of the lowenergy, subcatastrophic collisions considered here.
Conclusions. Subcatastrophic impacts have the potential to alter the shape of a small body significantly, without leading to major heating or compaction. The currently observed shapes of cometary nuclei, such as 67P, may be a result of such a major shape forming impact.
Key words: comets: general / comets: individual: 67P/ChuryumovGerasimenko / Kuiper belt: general / planets and satellites: formation
© ESO, 2016
1. Introduction
Whether cometary nuclei structures as observed today are pristine and preserve a record of their original accumulation, or are a result of later collisional or other evolutionary processes is still much debated (e.g. Weissman et al. 2004; Mumma et al. 1993; Sierks et al. 2015; Rickman et al. 2015; Morbidelli & Rickman 2015; Davidsson et al. 2016).
Based on data from the European Space Agency’s Rosetta rendezvous mission (Sierks et al. 2015), it was suggested that the particular bilobe structure of comet 67P/ChuryumovGerasimenko (67P) was formed during the early stages of the solar system (Massironi et al. 2015), possibly by low velocity accretionary collisions (Jutzi & Asphaug 2015) and therefore should be considered as a primordial body. On the other hand, Morbidelli & Rickman (2015) show that in the “standard scenario” of the early dynamical evolution of the solar system, an object of the size of 67P would have experienced a high number of catastrophic collisions and thus could not have survived. This study has been improved and a detailed analysis of the survival probability of a 67Plike object is presented in a companion paper (Jutzi et al. 2017, hereafter Paper I). It is found that even in the scenario without a longlasting primordial disc, 67P could not have conserved its primordial shape and a large number of shapechanging collisions should have occurred during its lifetime. This conclusion also holds true for generic bilobe structures, which might evolve further through a fissionmerging cycle as recently suggested by Scheeres et al. (2016). These results strongly indicate that neither the current shape of 67P nor its two individual lobes can actually be primordial.
Alternatively, 67Plike bilobe structures could be the result of collisional disruptions of somewhat larger bodies taking place at a later stage in the history of the solar system (e.g. Rickman et al. 2015; Morbidelli & Rickman 2015; Marchi et al. 2015, Paper I). During these later stages, typical relative velocities between small bodies have grown much larger than their mutual escape velocity (V ≫ V_{esc} ≈ 1 m/s) of kilometersized bodies and therefore direct bilobe formation by collisional mergers of similarsized bodies (Jutzi & Asphaug 2015) is no longer a viable mechanism. However, suitable low relative velocities could be found again between fragments/aggregates following a highervelocity collisional disruption of a larger parent body. In this secondary bilobe formation scenario, the dispersed material following a catastrophic collision reaccumulates in two gravitationally bound bodies which subsequently collide at low speed. Reaccumulation has been extensively studied for asteroids and results show that such behaviour is indeed occurring (e.g. Michel et al. 2001, 2003; Michel & Richardson 2013). A study of largescale disruption and subsequent reaccumulation in the case of cometary parent bodies is currently in progress (Schwartz et al., in prep.).
However, if the currently observed 67P structure did form as a result of a late catastrophic disruption of a larger parent body, the resulting shape must then survive until today. In other words, from the time of formation until present day it cannot have experienced a single subsequent shapechanging collision. In this context, the difference between the specific impact energy Q required to catastrophically disrupt the parent body of 67P and the energy sufficient to change its shape Q_{reshape} is a crucial quantity, as it determines the relative number of such events. On average, the larger the specific energy differences, the larger the number of shape changing impacts compared to the number of disruption events. Catastrophic disruptions are usually characterised by the specific impact energy , which leads to the escape of half of the initial mass involved in the collision. As it is shown in Paper I, the ratio /Q_{reshape} is very large, of the order of 10^{3}. Hence, there are on average many more shape changing collisions than disruptive collisions and it is unlikely that 67P could have conserved its shape, unless the disruption of the parent body has occurred very recently.
As we show in this paper, 67Plike bilobe structures could also emerge from low energy, subcatastrophic impacts on parent bodies of cometary size. It turns out that the specific energy Q_{sub} for these types of impacts is much closer to the specific reshape energy Q_{reshape} and hence the collision frequency difference between the two events is much smaller. This translates directly to a larger survival rate and therefore a higher probability of observing a shape such as that of 67P.
As suggested recently (Scheeres et al. 2016), comets with twocomponent shapes might enter a fission/merging cycle, once they enter the inner solar system and experience changes in their spin rate. Hence, its is possible that bilobe structures undergo episodic shape changes and the present day observed shape is simply the result of the last of such episodes. Even though this scenario adds the possibility for a time dependent macroscopic shape change, it is important to realise that a twocomponent object of cometary size must exist to begin with. Hence, even in this case the question of the longenough survival of bilobe shapes is a central question (Paper I).
Observational results, such as the abundance of supervolatiles (CO, CO_{2}, N_{2}) (e.g. Hässig et al. 2015; Le Roy et al. 2015), the detection of primordial molecules (Bieler et al. 2015) and the evidence for a low formation temperature (Rubin et al. 2015) suggest that 67P cannot have experienced any substantial, global scale heating after its formation. Further constraints include the high porosity and the observed homogeneity of the nucleus, which appears to be constant in density on a global scale without large voids (Pätzold et al. 2016). As a consequence, any proposed formation mechanism of bilobe shapes must be able to operate within very tight constraints on energy input and compaction of porous material.
In this paper we investigate the final shapes resulting from a new type of lowenergy, subcatastrophic impact on elongated, rotating bodies that meets the constraints mentioned above. We carry out a set of 3D smooth particle hydrodynamics simulations of impacts to investigate the possibility of forming bilobe structures reminiscent of those objected for cometary nuclei. In Sect. 2 we present our model approach and describe the setup and initial conditions. Representative bilobe forming collisions are presented in Sect. 3.1. The results (final shapes) of our parameter space exploration are shown in Sect. 3.2; heating and compaction effects are discussed in Sect. 3.3. In Sect. 4 we investigate the probabilities of a 67Plike structure being formed in such a manner to avoid destruction by subsequent shapechanging collisions.
Finally we discuss our bilobe formation model in the context of the question if comets are primordial (Sect. 5).
2. Bilobe formation by subcatastrophic impacts
2.1. Motivation and assumptions
Remote observations of cometary nuclei suggest that a large fraction of these objects have elongated rather than spherelike shapes (Lamy et al. 2004). This is interesting as it turns out that elongated bodies are naturally easier to “split” into two components than spherically symmetric bodies. This is even more important when they are rotating around their short axis, and centrifugal forces act in opposite directions at each end of the body. Impacts on such elongated rotation bodies might therefore act as a splitting mechanism leading first to two distinct bodies, which can potentially form a binary system or eventually merge together forming a bilobed body.
To study under which conditions such splitting can indeed produce such bilobed structures, we investigate the effects of impacts on rotating ellipsoids. We use axis ratios and rotation periods which are consistent with the observed range of values typical for cometary nuclei (Lamy et al. 2004). We consider impacts in the subcatastrophic regime for which most of the mass remains bound and accumulates on the main body (possibly made of two components). The impact scenario investigated here is very different from the case of a catastrophic breakup of a large parent body (Schwartz et al., in prep.) where the largest reaccumulated remnants contain only a small fraction of the initial mass. However, such remnants of catastrophic disruptions might have properties (elongated and rotating) similar to the targets considered here.
2.2. Modelling approach
The modelling approach used here is the same as in Paper I. We use a parallel smooth particle hydrodynamics (SPH) impact code (Benz & Asphaug 1995; Nyffeler 2004; Jutzi et al. 2008; Jutzi 2015) which includes selfgravity as well as material strength models. To avoid numerical rotational instabilities, the scheme suggested by Speith (2006) has been implemented.
In our modelling, we include an initial cohesion Y_{0}> 0 and use a tensile fracture model (Benz & Asphaug 1995), using parameters that lead to an average tensile strength Y_{T} ~ 100 Pa. To model fractured, granular material, a pressure dependent shear strength (friction) is included by using a standard DruckerPrager yield criterion (Jutzi 2015). As shown in Jutzi (2015), Jutzi et al. (2015), granular flow problems are well reproduced using this method. Porosity is modelled using a Palpha model (Jutzi et al. 2008) with a simple quadratic crush curve defined by the parameters P_{e}, P_{s}, ρ_{0}, ρ_{s0} and α_{0}. Following the approach in Paper I, we also introduce the density of the compacted material as ρ_{compact} = 1980 kg/m^{3} to define the initial distention α_{0} = ρ_{compact}/ρ_{0} = 4.5 corresponding to an initial porosity of 1−1 /α_{0} ~ 78%. We apply a modified Tillotson equation of state (EOS; e.g. Melosh 1989) with parameters appropriate for water ice. As in Paper I, we use a reduced bulk modulus, corresponding to an elastic wave speed of c_{e} ~ 0.1 km s^{1}.
The relevant material parameters used in the simulations are provided in Table 1.
Material parameters used in the simulations.
2.3. Setup and initial conditions
There is an infinite number of combinations of impact parameters in the case of nonspherical, rotating targets. Therefore, for practical reasons, we have to limit the size of the initial parameter space that can be investigated. The aim of this paper being more a demonstration of principle than a complete investigation of all possibilities, we limit ourselves to a few promising cases of relatively central collisions. However, a larger sample of various rotation rates and orientations, as well as a range of impact locations, have been investigated. The parameters defining the impact geometries used in this study are illustrated in Fig. 1.
Target and impactor have both the same initial material properties including their initial bulk density which is set to ρ_{0} ~ 440 kg/m^{3}.
We use two different ellipsoidal targets both of length L = 5.04 km but differing from one another by their axis ratios (0.4 and 0.7). With the assumed initial density, the two ellipsoids have initial masses of 4.7 × 10^{12} kg and 1.45 × 10^{13} kg, respectively. The impactor size is R_{p} = 100 m (in collisions involving the target with an axis ratio of 0.4) and R_{p} = 200 m (in collisions involving the target with an axis ratio of 0.7). The corresponding projectile masses are M_{p} = 1.8 × 10^{9} kg and M_{p} = 1.4 × 10^{10} kg, respectively. The impact velocities were chosen in the range of 200–300 m/s. We note that the sizes of the impactors as well as the collision velocities were motivated mainly by numerical rather than physical reasoning. Smaller impactors at higher speeds would not be well resolved (spatially) in our simulations and would also require smaller timesteps to avoid unphysical oscillations. Even at a relatively modest resolution (~10^{5} SPH particles), these simulations are very challenging even for a parallelised code. This is because the simulations have to extend over a very large number of dynamical timescales (and hence involve a very large number of timesteps) before the final structure of the resulting object can be determined (typically one day real time). We note, however, that the impact velocities considered here are supersonic (sound speed ~100 m/s) and the projectile small enough for the impact to be taking place in the socalled point source regime. Hence, our results can be scaled to larger impact speeds using appropriate scaling laws.
Fig. 1 Impact geometries. We use an ellipsoid with a constant length of L = 5.04 km and two different axis ratios (0.4 and 0.7). The ellipsoids are rotating with a rotation axis either along the y or z direction. Various impact points are investigated, as illustrated in the plot. The distance  dx  =  dy  = 670 m. 

Open with DEXTER 
Fig. 2 Comet 67P shape formation by subcatastrophic collisions. Shown is an example of an SPH calculation of an impact on a rotating ellipsoid. After the initial disruption, subsequent reaccumulation leads to the formation of two lobes. This processes may include the possible formation of layers. The two lobes are gravitationally bound and collide with each other within ~ one day forming a bilobed structure. The final shape is also shown in Fig. 4. Initial conditions: off axis (impact point + dy; see Fig. 1) impact of a R_{p} = 200 m impactor with a velocity of V = 300 m/s on a target with axis ratio 0.7 (mass M_{t} = 1.4 × 10^{13} kg) and rotation period of T = 6 h. The rotation axis is along the y axis (see Fig. 1). 

Open with DEXTER 
Fig. 3 As in Fig. 2 but showing a case with a more elongated target and different rotation axis and impactor properties. The final shape is shown in Fig. 5. Initial conditions: off axis (impact point + dx; see Fig. 1) impact of a R_{p} = 100 m impactor with a velocity of V = 250 m/s on a target with axis ratio 0.4 (mass M_{t} = 4.7 × 10^{12} kg) and rotation period of T = 9 h. The rotation axis is along the −y axis (see Fig. 1). 

Open with DEXTER 
Fig. 4 Shapes resulting from subcatastrophic disruptions of rotation ellipsoids. Shown are the results for different impact positions, rotation axes and periods (see Fig. 1 for the impact geometries). The impact velocity is V = 300 m/s and the impactor size R_{p} = 200 m. The initial target mass is 1.45 × 10^{13} kg and the target axis ratio is 0.7. The mass of the final bodies is of the order of ~70–80% of the initial mass. The case shown in Fig. 2 has the coordinates [+y(6h), +dy]. 

Open with DEXTER 
Fig. 5 Shapes resulting from subcatastrophic disruptions of rotation ellipsoids. Shown are the results for different impact positions, velocities and rotation periods and axes (see Fig. 1). The impact velocity is labelled in the plot; the impactor size is R_{p} = 100 m. The initial target mass is 4.7 × 10^{12} kg and the target axis ratio is 0.4. The mass of the final bodies is of the order of ~80–90 % of the initial mass. The case shown in Fig. 3 has the coordinates [–y(12h), +dx with v = 250/ms]. 

Open with DEXTER 
3. Results
3.1. Two examples of bilobeforming collisions
We present snapshots from the simulation of two representative cases of bilobeforming collisions in Figs. 2 and 3. Due to the elongated nature of the targets, the immediate postimpact mass distribution is concentrated at two locations. As a result, subsequent reaccumulation due to gravity leads to the formation of two distinct bodies. Centrifugal forces due to target rotation enhance the initial separation of these two masses. Material reaccumulating at low relative velocity during this process might also lead to the formation of layered structures on each individual body, such as observed on 67P (Massironi et al. 2015) and computed by Jutzi & Asphaug (2015). Finally, the two main gravitationally bound bodies eventually collide with each other at low velocity (determined by their mutual gravity and hence of order ~1 m/s) within approximately one day forming a 67Plike bilobe shape.
We note that the formation of a bilobe structure in such lowvelocity collisions (V ~ V_{esc}) of two gravitationally bound objects is consistent with previous results (Jutzi & Asphaug 2015). While these studies considered these lowvelocity impacts to occur during the early days of the solar system before small bodies became scattered by growing planets, the impact scenario presented here provides additional possibilities for such lowvelocity collisions to occur much later in the history of the solar system.
Due to numerical limitations (Sect. 2.3), the velocities considered are at the lower limit of the expected average velocities in the initial leftover planetesimal disc and are significantly lower than average velocity after disc dispersal (Morbidelli & Rickman 2015). However, the critical energies Q_{sub} for subcatastrophic bilobe formation for different impact velocities (and the corresponding projectile sizes) can be obtained using appropriate scaling (Sect. 3.3 below). It is important to point out that the postimpact focusing of the mass into two distinct locations, which has not been observed in previous simulations, is not due to the low velocity of the impacts, but rather the result of the specific target properties, namely their elongated shapes, the initial rotation and the relatively high porosity. All these properties are typical for real comets (Lamy et al. 2004), but were not considered in previous impact studies. Moreover, the collision energies considered here are between the cratering and the catastrophic regime, an area that has not yet been well explored.
3.2. Resulting final shapes
A summary of our simulations, using various rotation rates and orientations, as well as a range of impact points, is presented in Figs. 4 and 5, which show the final shapes resulting from the collisions considered. We find that for the conditions investigated here, there is a reasonably large fraction of shapes that have 67Plike bilobe structures. While this is true for the two different targets (different axis ratios) considered here, bodies with two distinct lobes are more probable outcomes in collisions involving the more elongated target.
3.3. Specific impact energies, heating and compaction
The specific impact energy is defined as (1)where μ_{r} = M_{p}M_{t}/ (M_{t} + M_{p}) is the reduced mass. We note that μ_{r} ≃ M_{p} for M_{p} ≪ M_{t}, as is the case here. In Fig. 6, we compare the specific energies Q_{sub} of the subcatastrophic impacts considered in this study to the specific impact energies for catastrophic collisions () as well as for shapechanging impacts on 67P (Q_{reshape}) (see Paper I). Following the scaling law already applied for and Q_{reshape}, we can write (2)with R = 2.52 km and μ = 0.42 (Paper I). The parameter a is determined by the impact conditions used for the two different targets (Sect. 2.3). For the case with the axis ratio of 0.4 we use an impact velocity of 250 m/s. Table 2 lists the values of parameter a for the various cases.
Parameters (SI units) for the scaling law Q_{crit} = aR^{3μ}V^{2−3μ}, where R is the target radius and V the impact velocity.
3.4. Impact heating
Also shown in Fig. 6 is the maximal global temperature increase dT resulting from the impacts. To estimate an upper limit for the global dT, we assumed that all kinetic impact energy is converted into internal energy: dT = Q/c_{p} for which a constant heat capacity c_{p} = 100 J/kg/K has been adopted (see Paper I). This approximate estimate already indicates that for impacts with energies comparable to Q = Q_{sub}, the maximal global temperature increase must remain relatively small. The actual temperature distribution depends on how much of the available kinetic energy is converted into heating and what fraction of the heated material remains on the body. Moreover, due to the highly dissipative characteristics of porous material, the compressed and therefore heated region remains very localised to the vicinity of the impact and therefore such an impact does not significantly affect the bulk content of volatile elements or the bulk porosity.
We demonstrate this in Fig. 7 by plotting the fraction of material that actually experienced a temperature increase larger than a certain dT for the two cases presented in Figs. 2 and 3. As can be seen, only ≲1% of the mass in the final bodies experienced a temperature increase larger than a few K. According to the scaling laws (Fig. 6) higher impact velocities would require slightly higher specific impact energies in order to result in a similar final bilobed configuration. For km s^{1} impact velocities, we would therefore expect the curves in Fig. 7 to be slightly shifted towards higher dT. However, as discussed above, the heating would remain localised and, on a global scale, dT remains limited to relatively small values even at high impact speeds. We note that the bulk of the dT values found here are generally much smaller than the sublimation temperatures of the observed supervolatiles, which are typically T_{s} > 20 K (Yamamoto 1985; Meech & Svoren 2004).
Fig. 6 Critical specific impact energies (adapted from Paper I): shown are specific energies for catastrophic disruptions (), shape changing collisions (Q_{reshape}) and subcatastrophic impacts (Q_{sub}). For the latter, a target radius R = 2.52 km and impact velocities of V = 300 m/s (axis ratio 0.7) and V = 250 m/s (axis ratio 0.4) have been assumed. 

Open with DEXTER 
3.5. Porosity evolution
Using the same procedure as in Paper I, we compute the cumulative distribution of the porosity in the final body, which takes into account compaction as well as the addition of macroporosity by reaccumulation of ejected material (Fig. 8).
Only a small fraction of the target mass (≲5%) significantly compacted (the porosity is reduced by ≳10%). On the other hand, a considerable amount of material experiences ejection and subsequent reaccumulation, thereby introducing macroporosity. Consequently, the final average porosities are slightly higher than the initial porosity (Fig. 8).
Figure 9 displays crosssections (in the xy plane) of the two bodies resulting from the impact simulations, showing the distribution of the final porosity. The porosity variations are generally relatively small, consistent with the observation that the interior of the nucleus is homogeneous on a global scale (Pätzold et al. 2016). On a large scale, porosity slightly decreases with depth. On smaller scales, layers of varying porosities are observed, suggesting that stratification, such as that observed on 67P (Massironi et al. 2015), may be produced by the reaccumulation of material on each single lobe (see also Figs. 2 and 3). In some locations, porosity decreases with increasing depths, as suggested by CONSERT measurements of the first ~100 m of the subsurface (Ciarletti et al. 2015). However, we stress that the spatial resolution of our simulations (of the order of ~100 m) does not allow for direct comparison of our results with these measurements. It also prevents us from making predictions regarding the size distribution of reaccumulated boulders.
Fig. 7 Fraction of material in the final body that experienced a temperature increase larger than a certain value dT. The curves correspond to the cases shown in Fig. 2 (axis ratio 0.7) and Fig. 3 (axis ratio 0.4), respectively. 

Open with DEXTER 
Fig. 8 Cumulative distribution of the porosity in the final body. The curves correspond to the cases shown in Fig. 2 (axis ratio 0.7) and Fig. 3 (axis ratio 0.4), respectively. The porosity calculation takes into account compaction as well as the increase of macroporosity. For comparison, the porosity distributions resulting from compaction only are shown as well. The final average porosity (compaction plus macroporosity by reaccumulation) is 85.1% (axis ratio 0.7) and 82.3% (axis ratio 0.4), while the initial porosity was 77.8%. 

Open with DEXTER 
Fig. 9 Crosssections showing the distribution of the porosity in the final body. The plots correspond to the cases shown in Fig. 2 (top) and Fig. 3 (bottom). 

Open with DEXTER 
4. Survival probabilities
4.1. Number of collisions
If the currently observed 67P structure did indeed form as a result of a collision event as shown above, for it to be observed today implies that it has not suffered any shapechanging collisions since the time of its formation.
Assuming that the formation of a bilobed structure of the size of 67P requires at least a specific impact energy of Q_{form}, we can estimate the average number of subsequent shapechanging collisions which have the minimum energy Q_{reshape} (as determined in Paper I), as well as the probability of avoiding all of these collisions. To allow for this, from Eq. (1) we calculate the minimal projectile radius delivering a specific impact energy Q ≥ Q_{min} as (3)where we assume that target and impactor have the same density and that M_{p} ≪ M_{t} (for impact velocities of a few hundred m/s, this is true even for ). The number of impacts on a target of size R_{t} by projectiles within the size range R_{min} ≤ R_{t} ≤ R_{max} during a time interval δt, can be written as (e.g. Morbidelli & Rickman 2015): (4)where P_{i} is the average intrinsic collision probability. Following Morbidelli & Rickman (2015) we use R_{max} = 50 km. N_{p}(R_{p}) represents the number of projectiles with a radius between R_{p} and R_{p} + dR_{p}. This distribution is not known precisely but in agreement with other studies of small body size distribution, we assume a differential size distribution of bodies given by dN/ dr ~ r^{q}. For the exponent q we use the same values as in calculations performed in Paper I: q = −2.5,−3.0 and −3.5.
Equation (4) can be used to compute the number of collisions N(Q_{min} = Q_{form}) within the time interval δt having a minimal specific impact energy Q_{form} allowing the formation of a 67Plike structure. In addition to these collisions, during the same time interval, a number of reshaping collisions N(Q_{min} = Q_{reshape}) will take place. Given that Q_{form} is larger than Q_{reshape}, N(Q_{min} = Q_{reshape}) will be larger than N(Q_{min} = Q_{form}). In other words, for one formation event (N(Q_{min} = Q_{form}) = 1) there will be on average N_{rs,norm} reshaping collisions. This can be expressed mathematically as: (5)or rewritten: (6)We compute N_{rs,norm} for the scenario of the formation of the 67P shape as a result of either a catastrophic breakup event (where we assume that a specific energy of at least is required) or a subcatastrophic formation event as illustrated above (where Q_{form} = Q_{sub}). The numbers obtained are shown in Table 3 for two different projectile velocities V = 500 m/s and V = 2 km s^{1}, which are representative of the relative velocities of small bodies within the planetesimal disc before and after dispersal, respectively (Morbidelli & Rickman 2015).
Average number of shapechanging collisions N_{rs,norm} for one formation event [N(Q_{min} = Q_{form}) = 1].
4.2. Probabilities of avoiding subsequent collisions
To compute the survival probability of a 67Plike structure formed by a collision event of a certain type, we can use the number of reshaping collisions per unit number of formation collisions N_{rs,norm} derived in the previous section. In the following, we estimate this survival probability P_{survival} for the different formation scenarios by assuming that a 67P structure was formed as a result of the last possible collision (the last collision involving the required specific energy (e.g. or Q_{form} = Q_{sub})). We then take into account that on average, only a fraction f of the shapechanging collisions N_{rs,norm} take place after the structure formed: (7)We further assume that these n_{f} events follow a Poisson distribution. In this case, the probability that all subsequent shapechanging collisions are avoided is given by (8)To obtain the total survival probability P_{survival} we integrate over all possible orders in which the collisions may take place (9)where we use the weight factor 1/(N_{rs,norm} + 1). By construction, P_{survival} = 1 for N_{rs,norm} = 0.
Survival probabilities in the different scenarios.
Fig. 10 Scenarios for the formation of a 67Plike bilobe structure. For each scenario, the probability of avoiding all subsequent shapechanging collisions is shown. For the cases of a late formation by a collision, we also indicate the probability that the parent body did not experience any prior catastrophic disruption. 

Open with DEXTER 
The results of these calculations are given in Table 4. We note that the calculations presented here do not take into account the probability of the formation of a bilobe shape in a specific impact event with a given energy. As discussed in Sect. 3.2, this probability is reasonably high, but is not possible to quantify.
For comparison, we also show in Table 4 the survival probabilities in the case of a primordial formation of a 67Plike bilobed structure (Paper I). For the standard scenario with a giant planet instability taking place at ~400 Myr, we adapt the number of disruptive collisions computed in Morbidelli & Rickman (2015) to the number of shapechanging collisions using Eq. (4) and Q_{reshape} (for a tensile strength of 100 Pa; Paper I). The survival probability then corresponds to the probability of avoiding all shapechanging collisions. For the scenario with no transient disc, the survival probabilities are computed in Paper I.
The results vary significantly between the different scenarios and also strongly depend upon the size distribution of available projectiles (characterised by q). Preliminary results from the New Horizons mission to Pluto and Charon based on the crater size frequency distribution suggest that q ~ −3.3 (Singer et al. 2015). We note that based on the most recent results it has been suggested that there may be a deficit of small objects (Singer et al. 2016); see discussion in Paper I, Sect. 4. For size distributions characterised by q ≥ 3, we find that P_{survival} is of the order of 1–10% in the framework of a formation following a subcatastrophic collision while P_{survival} < 0.1% in case of the formation directly following a catastrophic disruption. In the case of a primordial formation, the survival probability P_{survival} < 10^{8}, even with the conservative assumption that there was no initial massive transient planetesimal disc. We note that alternative models of the dynamical evolution (Davidsson et al. 2016, see discussion in Paper I) predict a much smaller collisional evolution and prefer shallower slopes of the size distribution, which would generally increase the survival probabilities computed here. An overview of the various scenarios considered in our study is presented in Fig. 10.
Note that so far we have not taken into account the possibility of a final shape evolution via a fissionmerging cycle (Scheeres et al. 2016). In this scenario, the bilobe forming collisions presented here would not form 67P directly, but rather a twocomponent structure with the right size ratio. This twocomponent body would then later evolve into the final shape via (one or several) fissionmerging cycle(s) once the comet enters the inner solar system. In this scenario, the survival probabilities might be somewhat larger than those computed above since general twocomponent structures require a slightly higher impact energy to be “destroyed” (Paper I).
It is likely that the parent body from which the 67Plike structure formed is not primordial itself (as indicated by the survival probabilities in Fig. 10) and has already experienced some collisional evolution or may be the result of a collisional cascade.
5. Discussion and conclusions
The analysis of the survival of the global structure of 67P (Paper I) strongly suggests that such a shape cannot be primordial. It must have formed as a result of a collision at a subsequent time (most probably within the last Gyr). At this time, the relative velocities between cometarysized bodies are such (V ≫ V_{esc}) that the formation mechanism invoked previously, namely the collisional mergers at low velocity (V ≈ V_{esc}) of similarsized bodies (Jutzi & Asphaug 2015), can no longer work directly.
In this paper, we present an alternative scenario. We investigate the final shapes resulting from a new type of lowenergy, subcatastrophic impact on elongated, rotating bodies, using a 3D SPH shock physics code. Our modelling results suggest that such collisions result in “splitting” events which frequently lead to formation of bilobe 67Plike shapes. This mechanism might not only explain the bilobe shape of some cometary nuclei but could potentially also provide an explanation for structures observed in the asteroid population, as for instance the particular forms of the asteroids 25143 Itokawa or 4179 Toutatis.
According to our model, comets are not primordial in the sense that their shape and structure formed during the initial stages of the formation of the solar system. Rather, the final structure is the result of the last major shapeforming impact. The subcatastrophic collisions investigated here provide the possibility of bilobe formation with small impact energies. Such smallscale impacts are much more frequent than catastrophic disruptions and the probability for such a shapeforming event to occur without a subsequent shapedestroying event occurring until today is estimated to be reasonably high. We note that the twocomponent structure resulting from the type of collisions investigated here might further evolve by fissionmerging cycles once the comet enters the inner solar system, as suggested recently (Scheeres et al. 2016).
Although the individual collisions considered in this work can alter the global shape, their respective energy is small enough not to lead individually to any substantial globalscale heating or compaction. In this sense, our formation model is consistent with the observed “pristinity” of 67P (e.g. Rubin et al. 2015; Bieler et al. 2015). However, it is likely that the parent body from which the 67Plike structure ultimately formed must have undergone significant collisional evolution (Fig. 10; see also Paper I), or is itself the result of a collisional disruption of a larger parent body. Several paths through the collisional cascade being able to lead to the similarsized bodies, the cumulative effects of impact heating and compaction experienced during the 4.6 Gyr of evolution by the material components that eventually form the comets observed today are difficult to establish with certainty. This formation degeneracy implies that it is not possible to uniquely reconstruct the detailed collision history nor the number and sizes of the parent bodies. Nevertheless, an upper limit for the size of parent bodies is given by the fact that larger bodies are subjected to internal heating by shortlived radionuclides (e.g. Prialnik et al. 2008) that significantly alter the pristine nature of the material and are therefore incompatible with observations of porosity and content in highly volatile elements in comets. Interestingly, this upper limit of the size of the parent body coupled with the requirement that the cumulative effects of impacts in terms of compaction and volatile losses can only be very moderate provide strong constraints for the duration and intensity of the collisional bombardment.
Whether these constraints are compatible with a scenario of a massive planetesimal disc phase existing for 450 Myr, as proposed by the Nice model (Tsiganis et al. 2005; Gomes et al. 2005; Morbidelli et al. 2012), remains to be analysed in detail. We note that Davidsson et al. (2016) suggest that the number of objects in the disc was much smaller, leading to less collisions (however, this model has other issues; see discussion in Paper I). On the other hand, our analysis of heating and porosity evolution in the impacts considered here, as well the regimes investigated in Paper I (shapechanging impacts as well as catastrophic disruptions), indicates that collisionally processed objects may still appear “primitive”. It is found that such bodies can still have a high porosity and could have retained their volatiles since these collisions generally do not lead to largescale heating of the material bound in the largest remnant. A more detailed study of the outcome of largescale catastrophic disruptions is currently in progress (Schwartz et al., in prep.).
Nevertheless, given our current understanding of the dynamics of the small bodies in the outer solar system, it is unlikely that the currently observed shape of 67P is primordial (even in the hypothetical scenario in which no initial massive planetesimal disc existed; Paper I). According to the calculations presented here, it may have formed as a result of the last major shapeforming impact. Even so, should future investigations show that the collisional cascade does not preserve the pristine nature of cometary material, we would be facing the conclusion that the current knowledge of the dynamics of small bodies in the outer regions of the solar system is seriously flawed. In this sense, comets provide invaluable tools for probing the origin and evolution of our solar system.
Acknowledgments
M.J. and W.B. acknowledge support from the Swiss NCCR PlanetS. We thank the referees B. Davidsson and N. Movshovitz for their thorough review which helped to improve the paper substantially.
References
 Benz, W., & Asphaug, E. 1995, Comput. Phys. Common., 87, 253 [NASA ADS] [CrossRef] [Google Scholar]
 Bieler, A., Altwegg, K., Balsiger, H., et al. 2015, Nature, 526, 678 [NASA ADS] [CrossRef] [Google Scholar]
 Ciarletti, V., LevasseurRegourd, A. C., Lasue, J., et al. 2015, A&A, 583, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Davidsson, B. J. R., Sierks, H., Güttler, C., et al. 2016, A&A, 592, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Nature, 435, 466 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Hässig, M., Altwegg, K., Balsiger, H., et al. 2015, Science, 347, aaa0276 [Google Scholar]
 Jutzi, M. 2015, Planet. Space Sci., 107, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Jutzi, M., & Asphaug, E. 2015, Science, 348, 1355 [NASA ADS] [CrossRef] [Google Scholar]
 Jutzi, M., Benz, W., & Michel, P. 2008, Icarus, 198, 242 [NASA ADS] [CrossRef] [Google Scholar]
 Jutzi, M., Holsapple, K. A., Wünnemann, K., & Michel, P. 2015, in Asteroids IV (University of Arizona Press), 679 [Google Scholar]
 Jutzi, M., Benz, W., Toliou, A., Morbidelli, A., & Brasser, R. 2017, A&A, 597, A61 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lamy, P. L., Toth, I., Fernandez, Y. R., & Weaver, H. A. 2004, in Comets II, eds. M. C. Festou, H. U. Keller, & H. A. Weaver (Tucson: University of Arizona Press), 223 [Google Scholar]
 Le Roy, L., Altwegg, K., Balsiger, H., et al. 2015, A&A, 583, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Marchi, S., Rickman, H., Massironi, M., et al. 2015, Lunar and Planetary Science Conference, 46, 1532 [NASA ADS] [Google Scholar]
 Massironi, M., Simioni, E., Marzari, F., et al. 2015, Nature, 526, 402 [NASA ADS] [CrossRef] [Google Scholar]
 Meech, K. J., & Svoren, J. 2004, in Comets II, eds. M. C. Festou, H. U. Keller, & H. A. Weaver (Tucson: University of Arizona Press), 317 [Google Scholar]
 Melosh, H. J. 1989, Impact cratering: A geologic process (New York: Oxford University Press) [Google Scholar]
 Michel, P., & Richardson, D. C. 2013, A&A, 554, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Michel, P., Benz, W., Tanga, P., & Richardson, D. C. 2001, Science, 294, 1696 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Michel, P., Benz, W., & Richardson, D. C. 2003, Nature, 421, 608 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Morbidelli, A., & Rickman, H. 2015, A&A, 583, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morbidelli, A., Marchi, S., Bottke, W. F., & Kring, D. A. 2012, Earth Plan. Sci. Lett., 355, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Mumma, M. J., Weissman, P. R., & Stern, S. A. 1993, in Protostars and planets III, 1177 [Google Scholar]
 Nyffeler, B. 2004, Ph.D. Thesis, University of Bern, Switzerland [Google Scholar]
 Pätzold, M., Andert, T., Hahn, M., et al. 2016, Nature, 530, 63 [NASA ADS] [CrossRef] [Google Scholar]
 Prialnik, D., Sarid, G., Rosenberg, E. D., & Merk, R. 2008, Space Sci. Rev., 138, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Rickman, H., Marchi, S., A’Hearn, M. F., et al. 2015, A&A, 583, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rubin, M., Altwegg, K., Balsiger, H., et al. 2015, Science, 348, 232 [NASA ADS] [CrossRef] [Google Scholar]
 Scheeres, D. J., Hirabayashi, T., Chesley, S., et al. 2016, Lunar and Planetary Science Conference, 47, 1615 [NASA ADS] [Google Scholar]
 Sierks, H., Barbieri, C., Lamy, P. L., et al. 2015, Science, 347, aaa1044 [Google Scholar]
 Singer, K. N., Schenk, P. M., Robbins, S. J., et al. 2015, in AAS/Division for Planetary Sciences Meeting Abstracts, 47, 102.02 [Google Scholar]
 Singer, K. N., McKinnon, W. B., Greenstreet, S., et al. 2016, AAS/Division for Planetary Sciences Meeting Abstracts, 48, 213.12 [Google Scholar]
 Speith, R. 2006, Habilitation, University of Tübingen, Germany [Google Scholar]
 Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Nature, 435, 459 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Weissman, P. R., Asphaug, E., & Lowry, S. C. 2004, in Comets II, eds. M. C. Festou, H. U. Keller, & H. A. Weaver (Tucson: University of Arizona Press), 337 [Google Scholar]
 Yamamoto, T. 1985, A&A, 142, 31 [NASA ADS] [Google Scholar]
All Tables
Parameters (SI units) for the scaling law Q_{crit} = aR^{3μ}V^{2−3μ}, where R is the target radius and V the impact velocity.
Average number of shapechanging collisions N_{rs,norm} for one formation event [N(Q_{min} = Q_{form}) = 1].
All Figures
Fig. 1 Impact geometries. We use an ellipsoid with a constant length of L = 5.04 km and two different axis ratios (0.4 and 0.7). The ellipsoids are rotating with a rotation axis either along the y or z direction. Various impact points are investigated, as illustrated in the plot. The distance  dx  =  dy  = 670 m. 

Open with DEXTER  
In the text 
Fig. 2 Comet 67P shape formation by subcatastrophic collisions. Shown is an example of an SPH calculation of an impact on a rotating ellipsoid. After the initial disruption, subsequent reaccumulation leads to the formation of two lobes. This processes may include the possible formation of layers. The two lobes are gravitationally bound and collide with each other within ~ one day forming a bilobed structure. The final shape is also shown in Fig. 4. Initial conditions: off axis (impact point + dy; see Fig. 1) impact of a R_{p} = 200 m impactor with a velocity of V = 300 m/s on a target with axis ratio 0.7 (mass M_{t} = 1.4 × 10^{13} kg) and rotation period of T = 6 h. The rotation axis is along the y axis (see Fig. 1). 

Open with DEXTER  
In the text 
Fig. 3 As in Fig. 2 but showing a case with a more elongated target and different rotation axis and impactor properties. The final shape is shown in Fig. 5. Initial conditions: off axis (impact point + dx; see Fig. 1) impact of a R_{p} = 100 m impactor with a velocity of V = 250 m/s on a target with axis ratio 0.4 (mass M_{t} = 4.7 × 10^{12} kg) and rotation period of T = 9 h. The rotation axis is along the −y axis (see Fig. 1). 

Open with DEXTER  
In the text 
Fig. 4 Shapes resulting from subcatastrophic disruptions of rotation ellipsoids. Shown are the results for different impact positions, rotation axes and periods (see Fig. 1 for the impact geometries). The impact velocity is V = 300 m/s and the impactor size R_{p} = 200 m. The initial target mass is 1.45 × 10^{13} kg and the target axis ratio is 0.7. The mass of the final bodies is of the order of ~70–80% of the initial mass. The case shown in Fig. 2 has the coordinates [+y(6h), +dy]. 

Open with DEXTER  
In the text 
Fig. 5 Shapes resulting from subcatastrophic disruptions of rotation ellipsoids. Shown are the results for different impact positions, velocities and rotation periods and axes (see Fig. 1). The impact velocity is labelled in the plot; the impactor size is R_{p} = 100 m. The initial target mass is 4.7 × 10^{12} kg and the target axis ratio is 0.4. The mass of the final bodies is of the order of ~80–90 % of the initial mass. The case shown in Fig. 3 has the coordinates [–y(12h), +dx with v = 250/ms]. 

Open with DEXTER  
In the text 
Fig. 6 Critical specific impact energies (adapted from Paper I): shown are specific energies for catastrophic disruptions (), shape changing collisions (Q_{reshape}) and subcatastrophic impacts (Q_{sub}). For the latter, a target radius R = 2.52 km and impact velocities of V = 300 m/s (axis ratio 0.7) and V = 250 m/s (axis ratio 0.4) have been assumed. 

Open with DEXTER  
In the text 
Fig. 7 Fraction of material in the final body that experienced a temperature increase larger than a certain value dT. The curves correspond to the cases shown in Fig. 2 (axis ratio 0.7) and Fig. 3 (axis ratio 0.4), respectively. 

Open with DEXTER  
In the text 
Fig. 8 Cumulative distribution of the porosity in the final body. The curves correspond to the cases shown in Fig. 2 (axis ratio 0.7) and Fig. 3 (axis ratio 0.4), respectively. The porosity calculation takes into account compaction as well as the increase of macroporosity. For comparison, the porosity distributions resulting from compaction only are shown as well. The final average porosity (compaction plus macroporosity by reaccumulation) is 85.1% (axis ratio 0.7) and 82.3% (axis ratio 0.4), while the initial porosity was 77.8%. 

Open with DEXTER  
In the text 
Fig. 9 Crosssections showing the distribution of the porosity in the final body. The plots correspond to the cases shown in Fig. 2 (top) and Fig. 3 (bottom). 

Open with DEXTER  
In the text 
Fig. 10 Scenarios for the formation of a 67Plike bilobe structure. For each scenario, the probability of avoiding all subsequent shapechanging collisions is shown. For the cases of a late formation by a collision, we also indicate the probability that the parent body did not experience any prior catastrophic disruption. 

Open with DEXTER  
In the text 