© 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 bi-lobe structure of comet 67P/Churyumov-Gerasimenko (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 67P-like object is presented in a companion paper (Jutzi et al. 2017, hereafter Paper I). It is found that even in the scenario without a long-lasting primordial disc, 67P could not have conserved its primordial shape and a large number of shape-changing collisions should have occurred during its lifetime. This conclusion also holds true for generic bi-lobe structures, which might evolve further through a fission-merging 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, 67P-like bi-lobe 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 kilometer-sized bodies and therefore direct bi-lobe formation by collisional mergers of similar-sized bodies (Jutzi & Asphaug 2015) is no longer a viable mechanism. However, suitable low relative velocities could be found again between fragments/aggregates following a higher-velocity collisional disruption of a larger parent body. In this secondary bi-lobe formation scenario, the dispersed material following a catastrophic collision re-accumulates in two gravitationally bound bodies which subsequently collide at low speed. Re-accumulation 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 large-scale disruption and sub-sequent re-accumulation 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 shape-changing 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 ${\mathit{Q}}_{\mathrm{D}}^{\mathrm{\ast }}$, which leads to the escape of half of the initial mass involved in the collision. As it is shown in Paper I, the ratio ${\mathit{Q}}_{\mathrm{D}}^{\mathrm{\ast }}$/Q_reshape is very large, of the order of 10³. 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, 67P-like bi-lobe structures could also emerge from low energy, sub-catastrophic 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 re-shape 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 two-component 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 bi-lobe 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 two-component object of cometary size must exist to begin with. Hence, even in this case the question of the long-enough survival of bi-lobe shapes is a central question (Paper I).

Observational results, such as the abundance of supervolatiles (CO, CO₂, N₂) (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 bi-lobe 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 low-energy, sub-catastrophic 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 bi-lobe structures reminiscent of those objected for cometary nuclei. In Sect. 2 we present our model approach and describe the setup and initial conditions. Representative bi-lobe 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 67P-like structure being formed in such a manner to avoid destruction by subsequent shape-changing collisions.

Finally we discuss our bi-lobe formation model in the context of the question if comets are primordial (Sect. 5).

2. Bi-lobe formation by sub-catastrophic impacts

2.1. Motivation and assumptions

Remote observations of cometary nuclei suggest that a large fraction of these objects have elongated rather than sphere-like 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 bi-lobed body.

To study under which conditions such splitting can indeed produce such bi-lobed 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 sub-catastrophic 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 break-up of a large parent body (Schwartz et al., in prep.) where the largest re-accumulated 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 self-gravity 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 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 Drucker-Prager 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 P-alpha model (Jutzi et al. 2008) with a simple quadratic crush curve defined by the parameters P_e, P_s, ρ₀, ρ_s0 and α₀. Following the approach in Paper I, we also introduce the density of the compacted material as ρ_compact = 1980 kg/m³ to define the initial distention α₀ = ρ_compact/ρ₀ = 4.5 corresponding to an initial porosity of 1−1 /α₀ ~ 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.

2.3. Setup and initial conditions

There is an infinite number of combinations of impact parameters in the case of non-spherical, 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 ρ₀ ~ 440 kg/m³.

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¹² kg and 1.45 × 10¹³ 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⁹ kg and M_p = 1.4 × 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 time-steps to avoid unphysical oscillations. Even at a relatively modest resolution (~10⁵ 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 time-steps) 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 super-sonic (sound speed ~100 m/s) and the projectile small enough for the impact to be taking place in the so-called 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.

Fig. 2 Comet 67P shape formation by sub-catastrophic collisions. Shown is an example of an SPH calculation of an impact on a rotating ellipsoid. After the initial disruption, subsequent re-accumulation 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 bi-lobed structure. The final shape is also shown in Fig. 4. Initial conditions: off axis (impact point + dy; see Fig. 1) impact of a Rp = 200 m impactor with a velocity of V = 300 m/s on a target with axis ratio 0.7 (mass Mt = 1.4 × 1013 kg) and rotation period of T = 6 h. The rotation axis is along the y axis (see Fig. 1).

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 Rp = 100 m impactor with a velocity of V = 250 m/s on a target with axis ratio 0.4 (mass Mt = 4.7 × 1012 kg) and rotation period of T = 9 h. The rotation axis is along the −y axis (see Fig. 1).

Fig. 4 Shapes resulting from sub-catastrophic 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 Rp = 200 m. The initial target mass is 1.45 × 1013 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].

Fig. 5 Shapes resulting from sub-catastrophic 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 Rp = 100 m. The initial target mass is 4.7 × 1012 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].

3. Results

3.1. Two examples of bi-lobe-forming collisions

We present snapshots from the simulation of two representative cases of bi-lobe-forming collisions in Figs. 2 and 3. Due to the elongated nature of the targets, the immediate post-impact mass distribution is concentrated at two locations. As a result, subsequent re-accumulation 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 re-accumulating 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 67P-like bi-lobe shape.

We note that the formation of a bi-lobe structure in such low-velocity collisions (V ~ V_esc) of two gravitationally bound objects is consistent with previous results (Jutzi & Asphaug 2015). While these studies considered these low-velocity 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 low-velocity 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 left-over planetesimal disc and are significantly lower than average velocity after disc dispersal (Morbidelli & Rickman 2015). However, the critical energies Q_sub for sub-catastrophic bi-lobe 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 post-impact 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 67P-like bi-lobe 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 $\mathit{Q}\mathrm{=}\mathrm{0.5}{\mathit{\mu }}_{\mathrm{r}}{\mathit{V}}^{\mathrm{2}}\mathit{/}\mathrm{(}{\mathit{M}}_{\mathrm{t}}\mathrm{+}{\mathit{M}}_{\mathrm{p}}\mathrm{)}\mathit{,}$(1)where μ_r = M_pM_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 sub-catastrophic impacts considered in this study to the specific impact energies for catastrophic collisions (${\mathit{Q}}_{\mathrm{D}}^{\mathrm{\ast }}$) as well as for shape-changing impacts on 67P (Q_reshape) (see Paper I). Following the scaling law already applied for ${\mathit{Q}}_{\mathrm{D}}^{\mathrm{\ast }}$ and Q_reshape, we can write ${\mathit{Q}}_{\mathrm{sub}}\mathrm{=}\mathit{a}{\mathit{R}}^{\mathrm{3}\mathit{\mu }}{\mathit{V}}^{\mathrm{2}\mathrm{-}\mathrm{3}\mathit{\mu }}\mathit{,}$(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.

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 bi-lobed 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 super-volatiles, 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 (${\mathit{Q}}_{\mathrm{D}}^{\mathrm{\ast }}$), shape changing collisions (Qreshape) and sub-catastrophic impacts (Qsub). 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.

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 re-accumulation 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 re-accumulation, thereby introducing macroporosity. Consequently, the final average porosities are slightly higher than the initial porosity (Fig. 8).

Figure 9 displays cross-sections (in the x-y 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 re-accumulation 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 re-accumulated 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.

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 re-accumulation) is 85.1% (axis ratio 0.7) and 82.3% (axis ratio 0.4), while the initial porosity was 77.8%.

Fig. 9 Cross-sections 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).

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 shape-changing collisions since the time of its formation.

Assuming that the formation of a bi-lobed structure of the size of 67P requires at least a specific impact energy of Q_form, we can estimate the average number of subsequent shape-changing 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 ${\mathit{R}}_{\min }\mathrm{=}\mathrm{(}\mathrm{2}{\mathit{Q}}_{\min }\mathit{V}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{3}}{\mathit{R}}_{\mathrm{t}}\mathit{/}\mathit{V,}$(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 ${\mathit{Q}}_{\min }\mathrm{=}{\mathit{Q}}_{\mathrm{D}}^{\mathrm{\ast }}$). 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): $\mathit{N}\mathrm{=}{\mathit{P}}_{\mathrm{i}}\mathit{\delta t}{\mathrm{\int }}_{{\mathit{R}}_{\min }}^{{\mathit{R}}_{\max }}\mathit{\pi }\mathrm{(}{\mathit{R}}_{\mathrm{t}}\mathrm{+}{\mathit{R}}_{\mathrm{p}}{\mathrm{)}}^{\mathrm{2}}{\mathit{N}}_{\mathrm{p}}\mathrm{(}{\mathit{R}}_{\mathrm{p}}\mathrm{)}\mathrm{d}{\mathit{R}}_{\mathrm{p}}\mathit{,}$(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 67P-like 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: ${\mathit{N}}_{\mathrm{rs}\mathit{,}\mathrm{norm}}\mathrm{=}\frac{\mathit{N}\mathrm{(}{\mathit{Q}}_{\min }\mathrm{=}{\mathit{Q}}_{\mathrm{reshape}}\mathrm{)}\mathrm{-}\mathit{N}\mathrm{(}{\mathit{Q}}_{\min }\mathrm{=}{\mathit{Q}}_{\mathrm{form}}\mathrm{)}}{\mathit{N}\mathrm{(}{\mathit{Q}}_{\min }\mathrm{=}{\mathit{Q}}_{\mathrm{form}}\mathrm{)}}\mathit{,}$(5)or rewritten: ${\mathit{N}}_{\mathrm{rs}\mathit{,}\mathrm{norm}}\mathrm{=}\frac{\mathit{N}\mathrm{(}{\mathit{Q}}_{\min }\mathrm{=}{\mathit{Q}}_{\mathrm{reshape}}\mathrm{)}}{\mathit{N}\mathrm{(}{\mathit{Q}}_{\min }\mathrm{=}{\mathit{Q}}_{\mathrm{form}}\mathrm{)}}\mathrm{-}\mathrm{1.}$(6)We compute N_rs,norm for the scenario of the formation of the 67P shape as a result of either a catastrophic break-up event (where we assume that a specific energy of at least ${\mathit{Q}}_{\mathrm{form}}\mathrm{=}{\mathit{Q}}_{\mathrm{D}}^{\mathrm{\ast }}$ is required) or a sub-catastrophic 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).

4.2. Probabilities of avoiding subsequent collisions

To compute the survival probability of a 67P-like 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. ${\mathit{Q}}_{\mathrm{form}}\mathrm{=}{\mathit{Q}}_{\mathrm{D}}^{\mathrm{\ast }}$ or Q_form = Q_sub)). We then take into account that on average, only a fraction f of the shape-changing collisions N_rs,norm take place after the structure formed: ${\mathit{n}}_{\mathit{f}}\mathrm{=}\mathit{f}{\mathit{N}}_{\mathrm{rs}\mathit{,}\mathrm{norm}}\mathit{.}$(7)We further assume that these n_f events follow a Poisson distribution. In this case, the probability that all subsequent shape-changing collisions are avoided is given by $\mathit{P}\mathrm{\left(}\mathit{f}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}{\mathit{n}}_{\mathit{f}}}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\mathit{f}{\mathit{N}}_{\mathrm{rs}\mathit{,}\mathrm{norm}}}\mathit{.}$(8)To obtain the total survival probability P_survival we integrate over all possible orders in which the collisions may take place ${\mathit{P}}_{\mathrm{survival}}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{1}}\frac{\mathit{P}\mathrm{\left(}\mathit{f}\mathrm{\right)}}{{\mathit{N}}_{\mathrm{rs}\mathit{,}\mathrm{norm}}\mathrm{+}\mathrm{1}}\mathrm{d}\mathit{f,}$(9)where we use the weight factor 1/(N_rs,norm + 1). By construction, P_survival = 1 for N_rs,norm = 0.

Fig. 10 Scenarios for the formation of a 67P-like bi-lobe structure. For each scenario, the probability of avoiding all subsequent shape-changing 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.

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 bi-lobe 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 67P-like bi-lobed 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 shape-changing 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 shape-changing 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 sub-catastrophic 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 fission-merging cycle (Scheeres et al. 2016). In this scenario, the bi-lobe forming collisions presented here would not form 67P directly, but rather a two-component structure with the right size ratio. This two-component body would then later evolve into the final shape via (one or several) fission-merging 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 two-component structures require a slightly higher impact energy to be “destroyed” (Paper I).

It is likely that the parent body from which the 67P-like 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 cometary-sized bodies are such (V ≫ V_esc) that the formation mechanism invoked previously, namely the collisional mergers at low velocity (V ≈ V_esc) of similar-sized 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 low-energy, sub-catastrophic 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 bi-lobe 67P-like shapes. This mechanism might not only explain the bi-lobe 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 shape-forming impact. The sub-catastrophic collisions investigated here provide the possibility of bi-lobe formation with small impact energies. Such small-scale impacts are much more frequent than catastrophic disruptions and the probability for such a shape-forming event to occur without a subsequent shape-destroying event occurring until today is estimated to be reasonably high. We note that the two-component structure resulting from the type of collisions investigated here might further evolve by fission-merging 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 global-scale 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 67P-like 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 similar-sized 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 short-lived 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 (shape-changing 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 large-scale heating of the material bound in the largest remnant. A more detailed study of the outcome of large-scale 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 shape-forming 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

All Tables

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.
In the text

Fig. 2 Comet 67P shape formation by sub-catastrophic collisions. Shown is an example of an SPH calculation of an impact on a rotating ellipsoid. After the initial disruption, subsequent re-accumulation 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 bi-lobed structure. The final shape is also shown in Fig. 4. Initial conditions: off axis (impact point + dy; see Fig. 1) impact of a Rp = 200 m impactor with a velocity of V = 300 m/s on a target with axis ratio 0.7 (mass Mt = 1.4 × 1013 kg) and rotation period of T = 6 h. The rotation axis is along the y axis (see Fig. 1).

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 Rp = 100 m impactor with a velocity of V = 250 m/s on a target with axis ratio 0.4 (mass Mt = 4.7 × 1012 kg) and rotation period of T = 9 h. The rotation axis is along the −y axis (see Fig. 1).

In the text

Fig. 4 Shapes resulting from sub-catastrophic 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 Rp = 200 m. The initial target mass is 1.45 × 1013 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].

In the text

Fig. 5 Shapes resulting from sub-catastrophic 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 Rp = 100 m. The initial target mass is 4.7 × 1012 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].

In the text

	Fig. 6 Critical specific impact energies (adapted from Paper I): shown are specific energies for catastrophic disruptions (${\mathit{Q}}_{\mathrm{D}}^{\mathrm{\ast }}$), shape changing collisions (Qreshape) and sub-catastrophic impacts (Qsub). 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.
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.
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 re-accumulation) is 85.1% (axis ratio 0.7) and 82.3% (axis ratio 0.4), while the initial porosity was 77.8%.

In the text

	Fig. 9 Cross-sections 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).
In the text

	Fig. 10 Scenarios for the formation of a 67P-like bi-lobe structure. For each scenario, the probability of avoiding all subsequent shape-changing 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.
In the text