© The Authors 2023

1. Introduction

The Cretaceous-Paleogene extinction event wiped out about 75% of all species (Raup & Sepkoski 1982): all non-avian dinosaurs died out, while all ammonites, plesiosaurs, and mosasaurs disappeared from the seas (Fastovsky & Sheehan 2005). This left mammals and birds to dominate the land. Kelly & Dachille (1953) proposed and later Alvarez et al. (1980) and Smit (1980) confirmed the hypothesis that these events were caused by the impact of a large asteroid of 10–15 km diameter on the Earth 66 million years ago (Renne et al. 2013).

Near-Earth objects (NEOs) are asteroids or comets with a perihelion distance, q, less than 1.3 au. Near-Earth asteroids (NEAs), are divided into four groups¹ (Atira, Aten, Apollo, and Amor) based on their semi-major axis and eccentricity. Asteroids with a diameter of more than 140 m in the Apollo and Aten groups having Earth-crossing orbits, are classified as potentially hazardous asteroids (PHAs). We note that comets crossing the Earth’s orbit can also be considered PHAs.

Although the recurrence interval of an impact event with a diameter of 10–15 km is about 0.2 billion years, a smaller impact event with a diameter of 1 km can be expected to occur every 0.6 million years (Collins et al. 2005).

In a recent review by Lubin & Cohen (2023), six methods of mitigating PHA impact are discussed. The most promising technique is the use of an array of penetrators combined with nuclear explosives that disintegrate the asteroid impactor at the so-called interception event into a gravitationally unbound cloud of fragments with sufficient energy to spread, (see, e.g., Kaplinger et al. 2010; Sanchez Cuartielles et al. 2008; King et al. 2021). The shock waves caused by the impact of fragments of an asteroid with a diameter of 1 km can be decorrelated (the impact events are dispersed in time) in the Earth’s atmosphere to a sufficient degree to reduce the threat. Furthermore, assuming an intercept time of more than about 2.5 months, a fragment expansion speed of about 1 m s⁻¹ is sufficient to miss Earth. However, for a 10 km diameter asteroid, the amount of energy released by the shock wave and even the dust production becomes large enough to overwhelm our atmospheric shield. Therefore, a much earlier intercept will be required to make the cloud of fragments grow large enough, so that only a small fraction of them end up hitting the Earth.

To predict the number of fragments that hit Earth, Lubin & Cohen (2023) use an analytic approximation neglecting the effect of orbital dynamics on the trajectory of fragments. In this Letter, we show that the fragment cloud is strongly distorted by the orbital shear, so that classical analytical approximations cannot be used to estimate the size of the impact. Therefore, in agreement with the results of King et al. (2021), using a series of high-precision N-body simulations is necessary.

2. Numerical simulations

A fully gravitationally interacting many-body system is used to model the impact event. The impactor is resolved as 10⁵ spherical particles of equal mass and size. The individual orbits of the impactor fragments, perturbed by the Sun and the major bodies of the Solar System (planets, Moon, and Pluto), are computed with a high-precision direct N-body integrator. We performed two sets of simulations: (1) a gravitationally non-interacting impact cloud to reveal the effect of orbital shear and (2) a fully interacting impact cloud, for which the gravity of the fragments is also considered.

We used our GPU-aided code HIPERION², utilising an eighth-order Hermite scheme described in Nitadori & Makino (2008). Based on the predictor-corrector relative acceleration of bodies, a shared adaptive timestep method is used (see details in Appendix B). Every time step is controlled such that simulations provide better positional accuracy than the diameter of the impactor.

2.1. Ephemerides and impact events

To take into account the perturbation effects of the Solar System bodies on the impactor orbit, we obtained the position and velocity vectors of the planets, Pluto, and the Moon from the NASA JPL Horizons service³ at the time of a hypothetical impact event. Throughout the investigation, 12:00 am 24/06/2030 (JD 2462677.0) was used as the date of the impact event. The Astroquery Python package was used to retrieve the ephemerides (Ginsburg et al. 2019).

We considered an asteroid and a comet-type impactor, which have the same diameter but different internal densities and impact velocities (physical properties are shown in Table 1). We note that the diameter of the impactors is chosen to ensure their survival even in sungrazing orbits (Sekanina 2003). Central collisions between Earth and the impactors are assumed, meaning that the impactor’s velocity vector points towards the Earth’s center of mass.

The properties of the impact are defined by three parameters: v_imp as the impact velocity; δ as the angle between the velocity vector of the Earth and that of the impactor in the plane of the Earth’s orbit; and ϕ as the angle between the velocity vector of the impactor and the orbital plane of the Earth. δ = 0° or δ = 180° correspond to a head-on or a rear-end collision with Earth, respectively. The calculated orbital elements of the impactor models are presented in Appendix >A.

2.2. Backward-forward integration, fragment cloud, and rotation

For the modelling of an impact event, a backward integration in time is first performed by reversing the velocity of the impactor as well as all other Solar System bodies at the impact event. The difference between the interception date (t_int) and the impact event time (t_imp) defines the length of the backward integration. We examine a series of Δt = t_imp − t_int sampled over the interval 1/48 − 3 years by 144 different integration times. At the end of the backward integration, that is, at the interception, the velocity vectors of all celestial bodies are reversed again. Note that the backward-forward integration method is similar to the one used in King et al. (2021).

Using the particle packing algorithm presented in Baranau & Tallarek (2017), the impactor is disintegrated into 1.04739 × 10⁵ individual monodisperse particles assuming a uniform size distribution for simplicity) at the point of impact. The diameter of the individual fragments is about 18 m and the corresponding porosity is measured to be 0.39.

Fragments inherit the impactor velocity components at interception position. To avoid the collapse of the fragment cloud by its own gravity, an additional spherically symmetric velocity distribution proportional to the local escape velocity is superimposed. In the applied expansion model, the velocity vectors of the fragments are given by the equation:

(1)

where R_i and R_i are the position vector and distance of the ith fragment measured from the centre of mass of the impactor and v_esc(R_i) is the local escape velocity at R_i. Four expansion models with different intensities are investigated assuming v₀ = 1,  2,  5, and 10.

In order to model the effect of the axial rotation of the impactor, simulations are performed assuming a minimum rotational period, which is determined by the maximum rotational speed of an impactor made of debris and held in place by gravity alone. This period is calculated based on Pravec & Harris (2000) as , where ρ is the internal density of the impactor. The corresponding rotational velocity at the surface of the impactor with a diameter of D is .

We assume that each fragment receives a tangential kick due to axial rotation at the moment of expansion onset. The kick velocity,

(2)

is added to the velocity of the fragment Eq. (1), where R_i, ⊥ is the tangential vector at a position, R_i, and M(α, β) is the Euclidean rotation matrix. The orientation of the impactor’s rotation axis is defined by two angles: α measures the inclination from the normal to the impactor’s orbital plane, n_o, whereas β measures the angle of rotation along n_o. The effect of the rotation is studied with 50 × 50 different α and β values in the range [0, π].

The shape of the fragment cloud is determined as it approaches Earth. This is done by fitting a triaxial ellipsoid to the particle distribution. The principal axis transformation of the moment of inertia tensor (see details in Appendix C) is used to obtain the best-fit ellipsoid. The parameters to be determined are as follows: first Euler angle θ, the size of the axes of the ellipsoid (a,  b, and c), and the oblateness b/a and c/b of the cloud.

3. Results and discussions

3.1. Importance of orbital dynamics

Our simulations, in agreement with Sanchez Cuartielles et al. (2008), show that orbital shear transforms the initially spherical fragment cloud into a triaxial ellipsoid. The orientation of the largest semi-major axis of the best-fit ellipsoid changes during the cloud’s flight and remains in its orbital plane. We will show that the cloud is shaped by the orbital shear (pericentre or apocentre passages), resulting in a non-trivial dependence of the impact number on the date of interception. Figure 1 shows the fragment clouds immediately before the impact for the models C4 and A4, assuming two different dates of interception. Surprisingly, the size of the late intercept cloud is larger than the early one because it’s more elongated for model A4. As a result, more fragments can hit the Earth for an early interception.

Fig. 1. Fragment clouds immediately before the impact to the impact for models C4 (left) and A4 (right). Two interception dates are shown: Δt = 1.6 yr and 2.2 yr. Note: the Z axis scale is independent of the X and Y axis scales and the fragment clouds are shifted arbitrarily for better visualisation.

First, we consider a comet-type impactor approaching Earth on a hyperbolic orbit (model C4), assuming two interception dates Δt = 1.6 yr and 2.25 yr. The left panels of Fig. 2 show the variations in the parameters of the best-fit ellipsoid as the impactor approaches Earth. As the fragment cloud passes by the apocentre, it is deflected by ∼ − 25°, and later of ∼15° during its final approach to Earth (panel a1). The latter deflection is caused by the Earth’s orbital perturbation (see Fig. A.1). The semi-major axis of the cloud grows to 2.592 × 10⁴ km and 5.281 × 10⁴ km prior to the impact for the two models (panel b1). The fragment cloud becomes an oblate spheroid with c/b ≃ 1 and b/a ≃ 0.5 (panel c1). The semi-major axis of the ellipsoid grows monotonically, and the associated increase in volume is V/V₀ = 1.303 × 10¹³ and 3.579 × 10¹³. As a consequence an earlier interception results in a more diluted cloud (left panel in Fig. 1). The measured number of fragments that impact Earth, that is, the impact numbers are N_imp = 55 529 and 30 240 for Δt = 1.6 yr and 2.25 yr, respectively.

Fig. 2. Variation of the parameters of the ellipsoid that best fits the cloud of fragments as it approaches the Earth. Top panels: first Euler angle θ. Middle panels: size of the semi-major axis, a. Bottom panels: oblateness (ratios of the semi-minor axes b/a and c/b). Two interception dates are compared: Δt = 1.6 yr vs. 2.25 yr for models C4 and A4, Δt = 1.885 yr vs. 2.10 yr for model A1.

Next, we demonstrate the importance of the effect of the orbital shear for an asteroid-type impactor in an elliptical (e = 0.669) orbit in model A4, with Δt = 1.6 yr and 2.25 yr interception dates (middle panels of Fig. 2). We note that the impactor is at pericentre and apocentre during interception. In both models, the oblate cloud completes a full rotation during its orbit around the Sun (panel a2). The semi-major axis of the cloud at impact is a  =  4.023  ×  10⁵ km for the late and a  =  1.3086  ×  10⁵ km for the early interception date (panel b2). Although the fragment cloud expands for a longer period of time, it is almost an order of magnitude smaller before impact in the latter case (right panel in Fig. 1). As a result, the corresponding impact numbers are N_imp = 6951 and 21 210. This can be explained by the fact that the orbital velocity increases (decreases) between the apocentre and the pericentre (and vice versa), causing compression/decompression of the cloud (see previous findings of Sanchez Cuartielles et al. 2008). During the apocentre or pericentre passages, the cloud oblateness (b/c) changes severely and the shortest axis (c) decreases by two orders of magnitude (panels c2).

Finally, we consider a highly eccentric impact orbit (model A1) assuming Δt = 1.885 yr and 2.10 yr (right panels of Fig. 2). The clouds are intercepted at the pericentre and the apocentre, respectively. In these particular cases, the orbital period is short (P = 0.419 yr), thus several pericentre and apocentre passages occur before impact, which are associated with several full rotations of the cloud (panel a3) and a periodic compression and decompression (panel b3). Noticeably, the compression and the following decompression magnitude is more pronounced for the late interception. The final size of fragment cloud is smaller for the earlier interception, that is, a = 3.093 × 10⁵ km and a = 2.422 × 10⁵ km for Δt = 1.885 yr and 2.1 yr, respectively. The corresponding impact numbers are N_imp = 85 415 and N_imp = 103 601.

3.2. Interception dates and fragment’s self-gravity

The number of impacts N_imp as a function of t_int − t_imp assuming different expansion velocities (1 ≤ v₀ ≤ 10) are measured for all models presented in Table A.1 and compared with an analytical prediction of N_imp (Appendix D), as shown in Fig. E.1. The effect of fragment’s self-gravity is also investigated.

The analytical estimation and the numerical measurement of N_imp are very similar for comet-type impactors approaching on a hyperbolic trajectory (models C3 and C4), assuming v₀ ≥ 2 (panels C3 and C4 of Fig. E.1). We note, however, that the analytical predictions underestimate the numerical measurements by about a factor of 2 for self-gravitating models with v₀ = 1, while they slightly overestimate the numerical simulations for non-self-gravitating models. For an asteroid-type impactor on a hyperbolic orbit (model A3) or on an elliptic orbit (models A2 and A4), the numerical measurements show a growing discrepancy with Δt up to about a factor of 2 (panels A2 and A3 of Fig. E.1). We note that comet-type impactors in models C2–C4 do not show this phenomenon.

Models A1 and C1 demonstrate that N_imp cannot be predicted by a simple analytical estimation (e.g., with Eq. (D.1) for NEOs in the Aten group (panels A1 and C1 in Fig. E.1). We find that N_imp fluctuates as a function of Δt with an amplitude proportional to v₀ and a period comparable to the impactor’s orbital period. Variations in the final size of the fragment cloud caused by intercepting at different orbital positions can explain the variation in the impact number. This phenomenon is more robust for fast expanding clouds, namely, for v₀ > 2 models.

Generally, self-gravity of the fragments has a tendency to increase the number of impacts. This effect is strong for models assuming v₀ = 1, while it is very weak for a relatively fast cloud expansion, v₀ > 2, (see the mismatch of the dotted lines and symbols in Fig. E.1). This can be explained by the fact that the volume of the fragment cloud is larger (by a factor of 4.8 for model A4) if the fragment’s self-gravity of the cloud is neglected. The self-gravity of the fragment is much weaker for comet-type impactors because they are less massive than asteroid-type impactors.

3.3. Orientation of axial rotation matters

To investigate the effect of the impactor’s axial rotation, we run additional non-self-gravitating fragment cloud models that account for an additional velocity component due to rotation based on Eq. (2). The expansion velocity is v₀ = 2, in which case the non-self-gravitating models give an accurate solution (see above). Simulations are performed assuming Δt = 0.6 yr and 1 yr to study the combined effect of the interception date and axial rotation.

Figures F.1 and F.2 show the distribution of N_imp in the α − β plane, measured for asteroid-type and comet-type impactors, respectively. In general, a reduction of about 30% in N_imp is achieved by axial rotation with a minimum period. We find that N_imp depends on the orientation of the rotation axis (with a few per cent variation) and there is a well-defined minimum number of impacts for each scenario. A comparison of the left-hand and right-hand columns shows that the effect of the date of interception is weak for models A1, A3, C1, C2, and C3. However, for models A2, A4, and C4, we find that N_imp reaches its minimum value at different points, depending on the interception date.

4. Conclusions

Here, we present a series of numerical simulations predicting the number of fragments that hit Earth from disintegrated asteroid-type (v_imp = 20 km s⁻¹) and comet-type (v_imp = 40 km s⁻¹) impactors with a diameter of 1000 m. We studied the effect of orbital shear, the interception date, and the axial rotation of the impactor. Based on the numerical simulations carried out, we have drawn the following conclusions.

1. The fragment cloud is shaped by orbital shear, making an oblate spheroid (for a hyperbolic orbit) or a triaxial ellipsoid (elliptical orbit). The pericentre and apocentre passages modulate the size and the oblateness of the ellipsoid, as well as, the volume of the cloud, which affects the number of fragments that hit the Earth.

2. To predict the number of impacts on Earth the orbital dynamics must be correctly calculated. Self-gravity of the fragments is important if the cloud expands at a rate slower than twice the local escape velocity. In some cases (180° ≤δ ≤ 270°), the analytical approximation is feasible for comet-type impactors if the cloud expands faster than the local escape velocity.

3. If the impactor is on a highly eccentric orbit, premature disassembly is undesirable. The best time to disintegrate the impactor is when it is at the pericentre of its orbit.

4. The impact number is affected by the orientation of the spin axis of the impactor prior to disassembly. For a well-defined orientation, which can depend on the date of interception, a minimum number of impacts can be found.

Finally, since early interception may not work in all cases, NEOs in the Apollo or Aten groups with orbital periods longer than one year on highly eccentric orbits may pose a major threat. To minimise the number of fragments hitting Earth, a sensitive interception timing is necessary. The impact hazard could be mitigated by the proven existence of an optimal orientation of the impactor’s rotational axis. We propose and explore a possible approach to this issue in a forthcoming study.

Acknowledgments

V.F. acknowledges financial support from the ESA PRODEX contract nr. 4000132054. The work of P.B. was supported by the Volkswagen Foundation under the special stipend No. 9B870. P.B. acknowledges the support within the grant No. AP14869395 of the Science Committee of the Ministry of Science and Higher Education of Kazakhstan (“Triune model of Galactic center dynamical evolution on cosmological time scale”). The work of P.B. was also supported under the special program of the NRF of Ukraine Leading and Young Scientists Research Support – “Astrophysical Relativistic Galactic Objects (ARGO): life cycle of active nucleus”, No. 2020.02/0346.

References

Appendix A: Impactors’ orbital elements

The orbital elements of the asteroid-type and comet-type impactors determined at the time of impact are given in Table A.1. However, it should be noted that the impactor does not approach on these orbits, as the Earth strongly perturbs the impactor’s orbit during a close encounter.

In order to determine the unperturbed orbit of the impactor, we conducted a backward integration in time, for instance, to the point of interception and then re-calculated the orbital elements. This way, the effect of the Earth’s perturbation is taken into account to correct the orbital elements. An illustration of orbital perturbations by the Earth’s gravitational field is given in Fig. A.1. There is no significant perturbation to the orbits of impactors on a head-on or rear-end collision trajectory (models A1, C1, A3, and C3). However, assuming a δ = 90° or a δ = 270° impact angle, the impactor’s trajectory is strongly perturbed, see the departure of fly path in Fig. A.1 for models A2, C2, A4, and C4.

Fig. A.1. Perturbation of impactor orbits during a close encounter with Earth. The impactor trajectories shown as dashed red lines (models A2, C2, A4, and C4) are strongly perturbed. However, the trajectories for head-on and rear-end impact (models A1, C1, A3, and C3) are not perturbed by the Earth.

It should be noted that the A1 and C1 models are special configurations and that in both cases the impactor is in its apocentre position at the impact event. In model A1, the Earth’s orbital speed is higher than that of the impactor, so the Earth will run into the impactor. In model C1, the impactor is travelling in a retrograde orbit, but the Earth’s orbital speed is still higher than that of the impactor when impacting.

Orbital elements determined at one year before impact are given in Table A.2. We emphasise that asteroid-type and one comet-type impactors resemble the Aten and Apollo group of NEOs which are indicated in Table A.2.

Figure A.2 shows the visualisation of the orbits of the hypothetical impactors whose orbital elements are presented in Table A.2. The orbits of planets of the solar system are also shown. The JPL Custom Orbital Visualization tool⁴ was used to create the orbital plots.

Fig. A.2. Visualisation of the orbits of the impactors (yellow): left for the asteroid-type, right for the comet-type impactors. The orbits of the inner planets are also shown. Red arrow indicates the rotation of the planets around the Sun, while black arrow indicates the fly path of the impactor.

Appendix B: Integrator order and precision

The adaptive timestep, dt⁽¹⁾, is recalculated at each step according to Nitadori & Makino (2008) as

(B.1)

where dt⁽⁰⁾ is the timestep at the previous step, and are the acceleration of the i-th body at the predictor and previous stage, and p is the order of the integrator.

We conducted a series of integrator precision measurements to determine the optimal order with acceptable position error. The impactor is modelled as a single body omitting the particle distribution generation. All major bodies in the Solar System are included in the integration. After 3 years of backward and forward integration, the impactor should be at the same position as it initially was, that is, at the impact point. The hypothetical impact event for the precision error measurements assumes an asteroid model with parameters of δ = 0, ϕ = 0, v_imp = 20 km s⁻¹. Due to the finite number representation used for computation (32 bit double precision) and the applied order of the Hermite scheme, the calculation resulted measurable position error.

By measuring this position error, we can determine the precision of the applied integrator. Figure B.1 shows our measurements for three different integrator orders. With the smallest η = 0.001 setting, the fourth-order scheme results in about 100 metre precision, while the sixth and eighth order gives about a metre precision. We note, however, that η = 0.001 results in a very small time steps, that is, extremely slow integration. With η = 0.05 the eighth-order scheme gives a position error below the size of the impactor (1000 m). However, using the fourth- or sixth-order scheme the precision grows several orders above the size of the impactor. Thus, to integrate with a relative position error less then the physical size of the impactor, the eighth-order scheme is selected with η = 0.05 throughout this study.

Fig. B.1. Measurements of position errors with fourth-, sixth-, and eighth-order Hermite schemes using different η parameters for the adaptive timestep calculation. The size of the impactor is also shown with a horizontal line. The optimal scheme with the desired precision (eighth order with η = 0.05) is indicated with red circle.

Appendix C: Size and shape of the fragment cloud

To describe the shape of the fragment cloud, a best fit ellipsoid is determined for its particles’ distribution (see an example in Fig C.1). This is done by calculating the cloud’s moment of inertia tensor. The routine def-triax.c is based on the well-known Jacobi eigenvalue algorithm⁵, which is an iterative method of rotation for the calculation of the eigenvalues and eigenvectors of a real symmetric matrix. For a review, we refer to the paper by Golub & van der Vorst (2000) on the topic of efficient and high performance eigenvalue computation algorithms.

Fig. C.1. Example of the best-fit triaxial ellipsoid model of a fragment cloud: a is the semi-major axis, and b,  c are the two semi-minor axes of the ellipsoid. The fragment distribution is taken from model C4.

For each snapshot, we determined the half mass radius r_hm of the particle distribution (simply sorting by the individual distances of the particles in the cloud from the centre of the mass of the whole system). In our further shape investigation, we used only the particles which we find inside the 4 × r_hm.

In our case, the matrix of inertia (MOI) in a centre-of-mass frame for the set of N particles by definition is a symmetric and real matrix:

(C.1)

Here m_i, x_i, y_i, and z_i stand for the i-th particle’s mass and Cartesian coordinates. A summation is performed for all particles found inside 4 × r_hm

After the Jacobi rotations (which in the end define for us a resulting first Euler angle θ), the MOI matrix contains only diagonal components:

(C.2)

In the next step, we find the maximum - I_max, middle - I_mid and minimum - I_min values from these diagonal elements. The axis ratio of the particle distribution ellipsoid can then easily be defined as:

(C.3)

The actual code is available from the authors upon an email request.

Appendix D: Impact number, analytical estimation

The number of fragments that hit Earth, N_imp, can be estimated by counting the fragments in a spherically symmetric homogeneous fragment cloud model, see Fig. D.1. Here, we re-introduce the method presented in Lubin & Cohen (2023). For simplicity, orbital dynamics (which can distort the spherical symmetry) and self-gravity (which can distort the homogeneity) are neglected in this model. The radius of the expanding cloud of fragments is given by R_clo(t, v₀) = R₀ + v(D/2)t at a given time t, where v(D/2) is the speed at which the cloud expands at the surface, and an D is the initial diameter of the impactor. Assuming a completely spherical fragment cloud, its volume will be V_clo(t, v₀) = (4/3)πR_clo(t, v₀)³ at a given time. The cylindrical section (whose diameter is equal to that of the Earth, represented by the white area in Fig. D.1) is , where R_⊕ is Earth’s radius. Assuming a homogeneous distribution of fragments in the cloud, the impact number (i.e. the number of fragments in the cylinder) is:

(D.1)

Fig. D.1. Fragment cloud before impact. The number of impacts can be estimated by counting the number of fragments that are inside the white cylinder.

with N_tot being the total number of fragments in the cloud. We note that due to the non-uniform expansion, the assumption of a homogeneous distribution of the fragment cloud fails for gravitationally non-self-interacting models with v₀ = 1, but can be plausible for gravitationally interacting cloud models.

Appendix E: Simulations on interception dates

The number of impacts on Earth as a function of Δt = t_int − t_imp with different values of expansion velocity, v₀ = 1,  2,  5, and 10 is presented in Fig. E.1. All asteroid-type (A1-A4) and comet-type (C1-C4) models are computed assuming impact direction δ = 0° , 90° , 180°, and 270°. Self-gravitating and non-self-gravitating models are also compared with the analytical predictions. Discussion is given in Sect. 3.2.

Fig. E.1. Number of fragments that hit Earth, Nimp, as a function of timp − tint for a set of impact models (assuming longitude of δ = 0° , 90° , 180° and 270°, see model details in Table 1). Left and right panels show asteroid-type (models A1-A4) and comet-type (models C1-C4) impactors, respectively. Models assume four cloud expansion velocities, v0 indicated with four distinct colours. The warm and cool colours represent the asteroid-type (vimp = 20 km s−1) and comet-type (vimp = 40 km s−1) impactors. The gravitationally interacting models are shown with symbols. The non-interacting impact cloud models are shown with thin solid and thin dashed lines, respectively. Thick solid lines represent the theoretical impact number given by equation D.1, for which case the effect of orbital dynamics on the cloud and the fragment’s self-gravity are neglected.

Appendix F: Simulations on rotation axis

Additional runs of non-self-gravitating fragment cloud models assuming v₀ = 2 which take into account the rotation of the impactor are presented in Figs F.1 and F.2. A discussion is given in Sect. 3.3.

Fig. F.1. Effect of the orientation of the rotational axis on the number of fragments that hit Earth, Nimp. The minimum rotational period is assumed (2.04 hours for a rubble pile composition corresponding to 0.42 m s−1 surface velocity) for an asteroid-type impactor. The integration length is 0.6 and 1 years for left and right panels, respectively. Models A1-A4 (impact angles of ϕ = 0 and δ = 0°  , 90° , 180° , and 270°) are investigated. The ratio of the minimum and maximum impact numbers, ΔNimp, are calculated for each model.

Fig. F.2. Same as Fig. F.1 for a comet-type impactor, i.e. computations for models C1-C4. In this case the minimum rotational period is 4.1 hours and the corresponding surface velocity is 0.21 m s−1.

All Tables

All Figures

	Fig. 1. Fragment clouds immediately before the impact to the impact for models C4 (left) and A4 (right). Two interception dates are shown: Δt = 1.6 yr and 2.2 yr. Note: the Z axis scale is independent of the X and Y axis scales and the fragment clouds are shifted arbitrarily for better visualisation.
In the text

	Fig. 2. Variation of the parameters of the ellipsoid that best fits the cloud of fragments as it approaches the Earth. Top panels: first Euler angle θ. Middle panels: size of the semi-major axis, a. Bottom panels: oblateness (ratios of the semi-minor axes b/a and c/b). Two interception dates are compared: Δt = 1.6 yr vs. 2.25 yr for models C4 and A4, Δt = 1.885 yr vs. 2.10 yr for model A1.
In the text

	Fig. A.1. Perturbation of impactor orbits during a close encounter with Earth. The impactor trajectories shown as dashed red lines (models A2, C2, A4, and C4) are strongly perturbed. However, the trajectories for head-on and rear-end impact (models A1, C1, A3, and C3) are not perturbed by the Earth.
In the text

	Fig. A.2. Visualisation of the orbits of the impactors (yellow): left for the asteroid-type, right for the comet-type impactors. The orbits of the inner planets are also shown. Red arrow indicates the rotation of the planets around the Sun, while black arrow indicates the fly path of the impactor.
In the text

	Fig. B.1. Measurements of position errors with fourth-, sixth-, and eighth-order Hermite schemes using different η parameters for the adaptive timestep calculation. The size of the impactor is also shown with a horizontal line. The optimal scheme with the desired precision (eighth order with η = 0.05) is indicated with red circle.
In the text

	Fig. C.1. Example of the best-fit triaxial ellipsoid model of a fragment cloud: a is the semi-major axis, and b,  c are the two semi-minor axes of the ellipsoid. The fragment distribution is taken from model C4.
In the text

	Fig. D.1. Fragment cloud before impact. The number of impacts can be estimated by counting the number of fragments that are inside the white cylinder.
In the text

Fig. E.1. Number of fragments that hit Earth, Nimp, as a function of timp − tint for a set of impact models (assuming longitude of δ = 0° , 90° , 180° and 270°, see model details in Table 1). Left and right panels show asteroid-type (models A1-A4) and comet-type (models C1-C4) impactors, respectively. Models assume four cloud expansion velocities, v0 indicated with four distinct colours. The warm and cool colours represent the asteroid-type (vimp = 20 km s−1) and comet-type (vimp = 40 km s−1) impactors. The gravitationally interacting models are shown with symbols. The non-interacting impact cloud models are shown with thin solid and thin dashed lines, respectively. Thick solid lines represent the theoretical impact number given by equation D.1, for which case the effect of orbital dynamics on the cloud and the fragment’s self-gravity are neglected.

In the text

Fig. F.1. Effect of the orientation of the rotational axis on the number of fragments that hit Earth, Nimp. The minimum rotational period is assumed (2.04 hours for a rubble pile composition corresponding to 0.42 m s−1 surface velocity) for an asteroid-type impactor. The integration length is 0.6 and 1 years for left and right panels, respectively. Models A1-A4 (impact angles of ϕ = 0 and δ = 0°  , 90° , 180° , and 270°) are investigated. The ratio of the minimum and maximum impact numbers, ΔNimp, are calculated for each model.

In the text

	Fig. F.2. Same as Fig. F.1 for a comet-type impactor, i.e. computations for models C1-C4. In this case the minimum rotational period is 4.1 hours and the corresponding surface velocity is 0.21 m s−1.
In the text