Issue 
A&A
Volume 677, September 2023



Article Number  L6  
Number of page(s)  12  
Section  Letters to the Editor  
DOI  https://doi.org/10.1051/00046361/202347205  
Published online  04 September 2023 
Letter to the Editor
Mitigating potentially hazardous asteroid impacts revisited
^{1}
Konkoly Observatory, Research Centre for Astronomy and Earth Science, KonkolyThege Miklós út 1517, 1121 Budapest, Hungary
email: regaly@konkoly.hu
^{2}
CSFK, MTA Centre of Excellence, Budapest, Konkoly Thege Miklós út 1517, 1121 Budapest, Hungary
email: frohlich.viktoria@csfk.org
^{3}
Eötvös Loránd University, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary
^{4}
Zentrum für Astronomie der Universität Heidelberg, Astronomisches RechenInstitut, Mönchhofstr. 1214, 69120 Heidelberg, Germany
^{5}
Main Astronomical Observatory, National Academy of Sciences of Ukraine, 27 Akademika Zabolotnoho St, 03143 Kyiv, Ukraine
email: berczik@mao.kiev.ua
Received:
15
June
2023
Accepted:
3
August
2023
Contact. Potentially hazardous asteroids (PHA) in Earthcrossing orbits pose a constant threat to life on Earth. Several mitigation methods have been proposed, and the most feasible technique appears to be the disintegration of the impactor and the generation of a fragment cloud by explosive penetrators at interception. However, mitigation analyses tend to neglect the effect of orbital dynamics on the trajectory of fragments.
Aims. We aim to study the effect of orbital dynamics of the impactor’s cloud on the number of fragments that hit the Earth, assuming different interception dates. We investigate the effect of selfgravitational cohesion and the axial rotation of the impactor.
Methods. We computed the orbits of 10^{5} fragments with a highprecision direct Nbody integrator of the eighth order, running on GPUs. We considered orbital perturbations from all large bodies in the Solar System and the selfgravity of the cloud fragments.
Results. Using a series of numerical experiments, we show that orbital shear causes the fragment cloud to adopt the shape of a triaxial ellipsoid. The shape and alignment of the triaxial ellipsoid are strongly modulated by the cloud’s orbital trajectory and, hence, the impact crosssection of the cloud with respect to the Earth. Therefore, the number of fragments hitting the Earth is strongly influenced by the orbit of the impactor and the time of interception. A minimum number of impacts occur for a welldefined orientation of the impactor rotational axis, depending on the date of interception.
Conclusions. To minimise the lethal consequences of an PHA’s impact, a wellconstrained interception timing is necessary. A tooearly interception may not be ideal for PHAs in the Apollo or Aten groups. Thus, we find that the best time to intercept PHA is when it is at the pericentre of its orbit.
Key words: celestial mechanics / Earth / minor planets / asteroids: general / methods: numerical
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1. Introduction
The CretaceousPaleogene extinction event wiped out about 75% of all species (Raup & Sepkoski 1982): all nonavian 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).
NearEarth objects (NEOs) are asteroids or comets with a perihelion distance, q, less than 1.3 au. NearEarth asteroids (NEAs), are divided into four groups^{1} (Atira, Aten, Apollo, and Amor) based on their semimajor axis and eccentricity. Asteroids with a diameter of more than 140 m in the Apollo and Aten groups having Earthcrossing 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 socalled 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^{−1} 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 highprecision Nbody simulations is necessary.
2. Numerical simulations
A fully gravitationally interacting manybody system is used to model the impact event. The impactor is resolved as 10^{5} 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 highprecision direct Nbody integrator. We performed two sets of simulations: (1) a gravitationally noninteracting 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 GPUaided code HIPERION^{2}, utilising an eighthorder Hermite scheme described in Nitadori & Makino (2008). Based on the predictorcorrector 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^{3} 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 comettype 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.
Parameters of the asteroid and comet type impactors.
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 headon or a rearend collision with Earth, respectively. The calculated orbital elements of the impactor models are presented in Appendix >A.
2.2. Backwardforward 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 backwardforward 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^{5} 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:
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_{0} = 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,
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 bestfit 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 semimajor axis of the bestfit 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 nontrivial 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 comettype 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 bestfit 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 semimajor axis of the cloud grows to 2.592 × 10^{4} km and 5.281 × 10^{4} 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 semimajor axis of the ellipsoid grows monotonically, and the associated increase in volume is V/V_{0} = 1.303 × 10^{13} and 3.579 × 10^{13}. 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 semimajor axis, a. Bottom panels: oblateness (ratios of the semiminor 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 asteroidtype 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 semimajor axis of the cloud at impact is a = 4.023 × 10^{5} km for the late and a = 1.3086 × 10^{5} 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^{5} km and a = 2.422 × 10^{5} 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 selfgravity
The number of impacts N_{imp} as a function of t_{int} − t_{imp} assuming different expansion velocities (1 ≤ v_{0} ≤ 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 selfgravity is also investigated.
The analytical estimation and the numerical measurement of N_{imp} are very similar for comettype impactors approaching on a hyperbolic trajectory (models C3 and C4), assuming v_{0} ≥ 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 selfgravitating models with v_{0} = 1, while they slightly overestimate the numerical simulations for nonselfgravitating models. For an asteroidtype 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 comettype 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_{0} 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_{0} > 2 models.
Generally, selfgravity of the fragments has a tendency to increase the number of impacts. This effect is strong for models assuming v_{0} = 1, while it is very weak for a relatively fast cloud expansion, v_{0} > 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 selfgravity of the cloud is neglected. The selfgravity of the fragment is much weaker for comettype impactors because they are less massive than asteroidtype impactors.
3.3. Orientation of axial rotation matters
To investigate the effect of the impactor’s axial rotation, we run additional nonselfgravitating fragment cloud models that account for an additional velocity component due to rotation based on Eq. (2). The expansion velocity is v_{0} = 2, in which case the nonselfgravitating 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 asteroidtype and comettype 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 welldefined minimum number of impacts for each scenario. A comparison of the lefthand and righthand 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 asteroidtype (v_{imp} = 20 km s^{−1}) and comettype (v_{imp} = 40 km s^{−1}) 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.

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.

To predict the number of impacts on Earth the orbital dynamics must be correctly calculated. Selfgravity 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 comettype impactors if the cloud expands faster than the local escape velocity.

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.

The impact number is affected by the orientation of the spin axis of the impactor prior to disassembly. For a welldefined 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.
See details at NASA’s Near Earth Objects Basics page: https://cneos.jpl.nasa.gov/about/neo_groups.html
High Precision Integrator for Nbody, HIPERION code can be accessed at https://konkoly.hu/staff/regaly/research/hiperion.html
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
 Alvarez, L. W., Alvarez, W., Asaro, F., & Michel, H. V. 1980, Science, 208, 1095 [NASA ADS] [CrossRef] [Google Scholar]
 Baranau, V., & Tallarek, U. 2017, J. Chem. Phys., 147, 224503 [NASA ADS] [CrossRef] [Google Scholar]
 Collins, G. S., Melosh, H. J., & Marcus, R. A. 2005, Meteorit. Planet. Sci., 40, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Fastovsky, D. E., & Sheehan, P. M. 2005, GSA Today, 15, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Ginsburg, A., Sipőcz, B. M., Brasseur, C. E., et al. 2019, AJ, 157, 98 [Google Scholar]
 Golub, G. H., & van der Vorst, H. A. 2000, J. Comput. Appl. Math., 123, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Kaplinger, B., Wie, B., & Dearborn, D. 2010, in Proceedings of the 60th International Astronautical Congress, Daejeon, Republic of Korea, 2009 [Google Scholar]
 Kelly, A., & Dachille, F. 1953, Target: Earth: The Role of Large Meteors in Earth Science (Pensacola: Pensacola Engraving Company) [Google Scholar]
 King, P. K., Syal, M. B., Dearborn, D. S., et al. 2021, Acta Astron., 188, 367 [NASA ADS] [CrossRef] [Google Scholar]
 Lubin, P., & Cohen, A. N. 2023, AdSpR, 71, 1827 [NASA ADS] [Google Scholar]
 Nitadori, K., & Makino, J. 2008, New Astron., 13, 498 [NASA ADS] [CrossRef] [Google Scholar]
 Pravec, P., & Harris, A. W. 2000, Icarus, 148, 12 [NASA ADS] [CrossRef] [Google Scholar]
 Raup, D. M., & Sepkoski, J. J. 1982, Science, 215, 1501 [NASA ADS] [CrossRef] [Google Scholar]
 Renne, P. R., Deino, A. L., Hilgen, F. J., et al. 2013, Science, 339, 684 [NASA ADS] [CrossRef] [Google Scholar]
 Sanchez Cuartielles, J. P., Vasile, M., & Radice, G. 2008, in 59th International Astronautical Congress, Glasgow, Scotland [Google Scholar]
 Sekanina, Z. 2003, ApJ, 597, 1237 [CrossRef] [Google Scholar]
 Smit, J. 1980, Nature, 285, 198 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Impactors’ orbital elements
The orbital elements of the asteroidtype and comettype 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.
Orbital elements of asteroidtype (A1A4) and comettype (C1C4) impactors determined at the impact event.
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 recalculated 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 headon or rearend 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 headon and rearend 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 asteroidtype and one comettype impactors resemble the Aten and Apollo group of NEOs which are indicated in Table A.2.
Orbital elements of asteroidtype (A1A4) and comettype (C1C4) impactors determined at the interception date.
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^{4} was used to create the orbital plots.
Fig. A.2. Visualisation of the orbits of the impactors (yellow): left for the asteroidtype, right for the comettype 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^{(1)}, is recalculated at each step according to Nitadori & Makino (2008) as
where dt^{(0)} is the timestep at the previous step, and are the acceleration of the ith 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^{−1}. 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 fourthorder 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 eighthorder scheme gives a position error below the size of the impactor (1000 m). However, using the fourth or sixthorder 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 eighthorder scheme is selected with η = 0.05 throughout this study.
Fig. B.1. Measurements of position errors with fourth, sixth, and eighthorder 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 deftriax.c is based on the wellknown Jacobi eigenvalue algorithm^{5}, 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 bestfit triaxial ellipsoid model of a fragment cloud: a is the semimajor axis, and b, c are the two semiminor 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 centreofmass frame for the set of N particles by definition is a symmetric and real matrix:
Here m_{i}, x_{i}, y_{i}, and z_{i} stand for the ith 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:
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:
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 reintroduce the method presented in Lubin & Cohen (2023). For simplicity, orbital dynamics (which can distort the spherical symmetry) and selfgravity (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_{0}) = R_{0} + 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_{0}) = (4/3)πR_{clo}(t, v_{0})^{3} 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:
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 nonuniform expansion, the assumption of a homogeneous distribution of the fragment cloud fails for gravitationally nonselfinteracting models with v_{0} = 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_{0} = 1, 2, 5, and 10 is presented in Fig. E.1. All asteroidtype (A1A4) and comettype (C1C4) models are computed assuming impact direction δ = 0° , 90° , 180°, and 270°. Selfgravitating and nonselfgravitating models are also compared with the analytical predictions. Discussion is given in Sect. 3.2.
Fig. E.1. Number of fragments that hit Earth, N_{imp}, as a function of t_{imp} − t_{int} 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 asteroidtype (models A1A4) and comettype (models C1C4) impactors, respectively. Models assume four cloud expansion velocities, v_{0} indicated with four distinct colours. The warm and cool colours represent the asteroidtype (v_{imp} = 20 km s^{−1}) and comettype (v_{imp} = 40 km s^{−1}) impactors. The gravitationally interacting models are shown with symbols. The noninteracting 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 selfgravity are neglected. 
Appendix F: Simulations on rotation axis
Additional runs of nonselfgravitating fragment cloud models assuming v_{0} = 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, N_{imp}. 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 asteroidtype impactor. The integration length is 0.6 and 1 years for left and right panels, respectively. Models A1A4 (impact angles of ϕ = 0 and δ = 0° , 90° , 180° , and 270°) are investigated. The ratio of the minimum and maximum impact numbers, ΔN_{imp}, are calculated for each model. 
Fig. F.2. Same as Fig. F.1 for a comettype impactor, i.e. computations for models C1C4. In this case the minimum rotational period is 4.1 hours and the corresponding surface velocity is 0.21 m s^{−1}. 
All Tables
Orbital elements of asteroidtype (A1A4) and comettype (C1C4) impactors determined at the impact event.
Orbital elements of asteroidtype (A1A4) and comettype (C1C4) impactors determined at the interception date.
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 semimajor axis, a. Bottom panels: oblateness (ratios of the semiminor 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 headon and rearend 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 asteroidtype, right for the comettype 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 eighthorder 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 bestfit triaxial ellipsoid model of a fragment cloud: a is the semimajor axis, and b, c are the two semiminor 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, N_{imp}, as a function of t_{imp} − t_{int} 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 asteroidtype (models A1A4) and comettype (models C1C4) impactors, respectively. Models assume four cloud expansion velocities, v_{0} indicated with four distinct colours. The warm and cool colours represent the asteroidtype (v_{imp} = 20 km s^{−1}) and comettype (v_{imp} = 40 km s^{−1}) impactors. The gravitationally interacting models are shown with symbols. The noninteracting 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 selfgravity are neglected. 

In the text 
Fig. F.1. Effect of the orientation of the rotational axis on the number of fragments that hit Earth, N_{imp}. 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 asteroidtype impactor. The integration length is 0.6 and 1 years for left and right panels, respectively. Models A1A4 (impact angles of ϕ = 0 and δ = 0° , 90° , 180° , and 270°) are investigated. The ratio of the minimum and maximum impact numbers, ΔN_{imp}, are calculated for each model. 

In the text 
Fig. F.2. Same as Fig. F.1 for a comettype impactor, i.e. computations for models C1C4. 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.