The Hazardous km-sized NEOs of the Next Thousands of Years

The orbits of the NEOs are propagated using the JPL small-body integrator which is based on an N-body model that includes Sun, planets, Pluto, Moon, and small-body perturbers (Farnocchia et al. 2015). When the Yarkovsky effect was detected from astrometric data (Farnocchia et al. 2013b), we added it to the force model. The ephemeris models used in the integration are DE441 (Park et al. 2021) for the planets and SB441-N16 for the largest main-belt bodies. ⁵

Figure 1 shows the propagation of the orbit of 2015 FP332, which reveals the relevant dynamical effects to the long-term dynamics of near-Earth objects. 2015 FP332 is in a Lidov-Kozai cycle (Lidov 1962; Kozai 1962), in which periods of high eccentricity are exchanged with periods of high inclination. In this case, both longitude of the node and argument of perihelion drift secularly.

Zoom In Zoom Out Reset image size

Planetary encounters can cause the exponential growth of the distance between initially neighboring trajectories, a necessary condition for chaos (Tancredi 1998). Neighboring trajectories of near-Earth asteroids can diverge in timescales ranging from decades, such as 99942 Apophis (Farnocchia et al. 2013a); to hundreds of years, such as 29075 (1950 DA) (Farnocchia & Chesley 2014); to tens of thousands of years, such as 433 Eros (Michel et al. 1996). In this process the linear approximation of the state uncertainty can quickly become inaccurate.

Figure 1 shows the propagation of the multiple samples or virtual asteroids of the orbit of 2015 FP332. In this example, the nominal trajectory of 2015 FP332 experiences a very close Venus encounter that causes the rapid increase in semimajor axis. Once each initially neighboring virtual asteroid diverges to a different trajectory it experiences a unique sequence of close encounters. This effect motivates the use the MOID to estimate long-term probabilities of collision. The resulting dynamics under these encounters are very nonlinear, and the orbits of near-Earth objects in these timescales become stochastic. For this reason, we sample the uncertainty in the orbits of NEOs and propagate them in a Monte Carlo simulation. The detection of potential impactors of small probabilities is out of the scope of this work, in which the main metric of interest is the MOID. For this reason we run a limited number of Monte Carlo samples (N = 21), which allows us to distinguish the main dynamical effects as well as the uncertainty in mean anomaly.

The presence of close encounters is expected if the near-Earth object has a small MOID with any of the planets. Thus, tracking the evolution of the MOID is not only relevant for the evaluation of the probability of collision with Earth but to understand when the dynamics are subject to nonlinear stochastic variations. The evolution of the orbit of 2015 FP332 in Figure 1 shows the effect of a low-MOID period in the long-term prediction. A Venus low-MOID enables close approaches that cause the rapid expansion of the distribution of orbits and mean anomaly to become unknown.

The MOID is the result of the optimization of the relative distance between two bodies over their respective fast angles. There are multiple algorithms to compute the MOID in the literature, including analytical methods (Gronchi 2005) and numerical methods such as Wiśniowski & Rickman (2013) or Hedo et al. (2018), which is used in this work. The MOID, a function of the osculating orbit elements, is then computed when postprocessing the numerically integrated trajectory.

Depending on the dynamical effects on the asteroid described in the previous section we find a variety of MOID long-term dynamics trajectories. Figure 2 shows a few examples of the dynamics of the propagation of the MOID for four km-sized NEOs. 2021 UY9 represents the simplest case, in which the MOID does not become small throughout the simulation time therefore making Earth impacts impossible. The uncertainty in the orbit of 2021 UY9 remains small throughout the propagation. The case of 29075 (1950 DA) is the very common case for NEOs, in which the MOID drifts secularly until a future zero crossing that lasts a short period of time, of about a century. The example of 136618 (1994 CN2) is similar to the case of 2015 FP332 in Figure 1, in which the date in which the MOID becomes small for the first time is uncertain. After a low-MOID period the trajectory becomes more uncertain. The last example, 2329 Orthos (1976 WA), illustrates the scenario of an extended period of time with a low Earth MOID. This is caused by the combination of two effects, a large amplitude of short-period oscillations and a favorable phasing of the secular cycle.

Zoom In Zoom Out Reset image size

The examples of Figure 2 are also representative of the growth in uncertainty of the MOID. In the two top cases the uncertainty remains small for thousands of years. The km-sized NEA population have well defined orbits and their Earth MOIDs remain well known for hundreds of years. Previous works focus on mapping the orbit covariance into a confidence region of the MOID (Gronchi & Tommei 2007). However, the uncertainty in the orbit can become far from Gaussian in long-term orbit propagations. Thus, we use statistics of the Monte Carlo propagation as indicators of the spread of the distribution as well as the confidence in our predictions.

The complexity in long-term MOID dynamics that we showed in the previous section motivates the development of a systematic method to quantify the long-term Earth impact hazard of NEOs. We propose a novel metric to characterize the potential impact hazard that consists in an estimated probability of collision between planet and NEO. The probability of collision is a problem primarily studied in two major timescales. The fundamental problem of impact hazard assumes the position of the asteroid within its orbit is reasonably well determined and it is possible to precisely determine the geometry of the subsequent close encounters. In the case of potentially hazardous asteroids, this analysis can typically be completed for one or two centuries (Chamberlin et al. 2001; Roa et al. 2021). In these timescales the uncertainty in the orbit of many NEOs starts to grow large enough that the position within the orbit can become unknown. This effect motivates the statistical assumption of a uniformly distributed mean anomaly.

Traditional impact probability theories assume that the orbit elements of the two objects involved are constant and have one intersection point (Öpik 1951; Wetherill 1967). Then, the probability depends on the timing of the orbits, which is when the mean anomalies are assumed uniformly distributed. This timing probability, here P_MA, is the probability that both bodies are in the right time at the right place, i.e., the range of mean anomalies that corresponds to a collision or flyby within a small distance.

There are a few options for the timing probability in the literature. In the most simplified case, we can assume the planet's orbit to be circular with Öpik's formula (Öpik 1951), in which the probability is function of the Keplerian elements of the asteroid a-e-i. Wetherill (1967) then derived an expression for an elliptic orbit of the planet, with the problem of being singular at zero inclination. In this work we use this expression as rederived recently in JeongAhn & Malhotra (2017) for regular nontangential encounters. In particular, we use the extended expression for the case in which the two objects do not exactly intersect. That means that the MOID is a positive value between 0 and the distance threshold for the close encounters of interest d. Thus, the probability that two objects have a close encounter with closest approach distance smaller than d is:

$$\begin{eqnarray}&&{P}_{\mathrm{MA}}=\displaystyle \frac{2{Ud}}{{T}_{p}{T}_{\mathrm{NEO}}| {{\boldsymbol{v}}}_{p}\times {{\boldsymbol{v}}}_{\mathrm{NEO}}| }\sqrt{1-\displaystyle \frac{{\mathrm{MOID}}^{2}}{{d}^{2}}}\end{eqnarray} \tag{ 1 }$$

where v _p and v _NEO are the velocities of the planet and the asteroid at the point that defines the MOID, U is the relative velocity at the same point, and T_p , T_NEO are the respective orbit periods. The square root term of Equation (1) adjusts the probability for a nonzero MOID. If MOID >d, the probability is assumed to be zero. This expression can be averaged for a MOID uniformly distributed between 0 and d. However, in this work we do not need to make this assumption as we keep track of the MOID throughout our long-term propagation and the distribution can be far from uniform in the range 0 < MOID ≤d.

Once we allow the orbit of the NEO to be time-varying, we can obtain the probability of collision as the combination of two terms: the probability that there is an intersection between the planet and near-Earth object and P_MA. If we investigate a potential Earth collision, the condition is that the Earth MOID is smaller than the combined radii of the two bodies considering gravitational focusing as required. The gravitational focusing factor virtually extends the radius of the planet to account for trajectories that lead to a collision due to the planet's gravity, and is a function of the incoming velocity of the asteroid V_∞ and the mass and radius of the planet M_p , R_p :

$$\begin{eqnarray}&&\gamma =\sqrt{1+\displaystyle \frac{2{{GM}}_{p}}{{R}_{p}{V}_{\infty }^{2}}}.\end{eqnarray} \tag{ 2 }$$

This approach has been used in the past to obtain the probability of collision for asteroids under the Lidov-Kozai cycle (Vokrouhlický et al. 2012; Pokorný & Vokrouhlický 2013). In that case, the generalized probability of collision is obtained as the sum over all the crossing configurations (noted with *) of the fraction of time that the NEO spends withing the distance threshold times the timing probability:

$$\begin{eqnarray}&&P=\displaystyle \sum _{* }{\left(\displaystyle \frac{{\rm{\Delta }}{t}_{\mathrm{MOID}\lt d}}{{T}_{\sec }}\right)}^{* }{P}_{\mathrm{MA}}(d,{K}_{p}^{* },{K}_{\mathrm{NEO}}^{* })\end{eqnarray} \tag{ 3 }$$

where K_P , K_NEO are the Keplerian elements of the planet and the NEO and Δt_MOID<d is the amount of time that the NEO spends within the distance threshold d. Vokrouhlický et al. (2012) and Pokorný & Vokrouhlický (2013) model the long-term asteroid dynamics with an analytical solution of the Lidov-Kozai cycle of the Jupiter perturbation, which defines the secular period T_sec as the Lidov-Kozai period. As a result, the fraction Δt and intersection configurations are computed analytically. As we show in the previous section, defining the times in which the MOID is small can be a complex problem under a wide range of dynamical contributions. In this work we propagate the orbit numerically to find the low-MOID periods. Because there is not a small discrete number of crossings along the long-term dynamics of the NEO, we estimate the probability as the average throughout the propagation time T using

$$\begin{eqnarray}&&P=\displaystyle \frac{1}{T}{\int }_{T}\kappa {P}_{{\rm{M}}{\rm{A}}}(d,{K}_{p},{K}_{{\rm{N}}{\rm{E}}{\rm{O}}})dt.\end{eqnarray} \tag{ 4 }$$

This integral is computed numerically using the numerically integrated trajectories. The factor κ is introduced so we can null the contribution of the trajectory in which the position of the object is deterministic within its orbit, i.e., when the uncertainty in mean anomaly is small. κ is set to 0 before the first date in which we find that the standard deviation in mean anomalies is larger than 10°, and set to 1 elsewhere. This distinction allows us to rule out the associated risk of objects that currently have a very low MOID but their position is properly constrained for the duration of their visit to the planet's vicinity. In addition, we check the close encounters that were recorded in the propagation before we can use the analytical expression for P_MA in Equation (1).

The orbits of the known km-sized NEO population are propagated starting from their orbit solutions in JPL's Small-body Database as of 2022 September 15. ⁶ For each NEO we find the first date at which a low-MOID period is found between all of the Monte Carlo samples, with a threshold defined as MOID < 0.01 au (235 Earth Radii or 3.89 Lunar Distances). At this threshold, the incoming velocity V_∞ required for a collision is of 0.05 km s⁻¹ or less, from solving Equation (2). From a statistical point of view, this relative velocity is extremely unlikely (Farnocchia & Chodas 2021; Harris & Chodas 2021).

The first question we answer is how many km-sized NEOs currently have a MOID < 0.01 au, and how will this number evolve in the next 1000 yr. As of the time in which the JPL's SBDB was queried, there are 40 NEOs that fulfill this condition. The evolution of this estimated number of NEOs is shown in Figure 3. As the uncertainty in the orbits of the NEOs grows into the future, only some of the MC samples may have MOID < 0.01 au. This phenomenon is shown in more detail in Figure 4, which shows the estimated range of km-sized NEOs with MOID < 0.01 au based in the Monte Carlo samples. The uncertainty remains very small (±1 body) throughout the next 500 yr. By the end of the millennium, this number is in the range of 26–72 km-sized NEOs. As mentioned earlier, none of these objects pose a collision threat to Earth in the next 100 yr.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Individual results of the MOID propagation are shown in Figure 4. We sorted the km-sized NEOs by the date in which they meet the MOID < 0.01 au condition. As defined by the length of their low-MOID periods, we show NEOs that are expected to be continuously in the vicinity of Earth as opposite to the ones that are for a brief period of time. We observe that even if the number of NEOs will never exceed an average value of 40–45 NEOs, the total number of unique low-MOID NEOs in the next 1000 yr is of almost 150.

We list the km-sized NEOs that frequently experience MOID < 0.01 au in Table 1 based in the fraction of the next 1000 yr that they meet the condition. Their implicit probability of collision is assessed in the next section. There are 4 objects whose MOID remains lower than 0.01 au throughout this millennium: 7482 (1994 PC1), 68950 (2002 QF15), 164121 (2003 YT1), and 144332 (2004 DV24). In the second and third case, the mean anomaly remains well defined throughout the millennium.

The propagation of the MOID of a few of the top-ranked NEOs is shown in Figure 5. Most remarkably, we can see how low the MOID of 7482 (1994 PC1) persists throughout the next 1000 yr. In 68950 (2002 QF15) and 164121 (2003 YT1) we observe a secular drift in the MOID. In the case of 68950 (2002 QF15), this secular drift predicts a near-zero MOID around year 2500. In the case of 164121 (2003 YT1), the MOID is increasing at a relatively slow rate. The last example, 143651 (2003 QO104), shows a large amplitude of the MOID around zero, which motivates additional analysis of the long-term hazard.

Zoom In Zoom Out Reset image size

In the previous section we inspected the necessary condition for very close encounters to occur: a MOID <0.01 au. The next step is to estimate the collision probability by making assumptions on their mean anomalies. The method to compute this timing probability once the MOID is low was described in Section 2.3. To study the potential impact hazard we set a smaller close approach threshold (1 LD) and take into account the deterministic parts of the NEO position during the orbit propagation. In addition, we study the list of close approaches generated in the Monte Carlo experiment to validate our predictions.

The analytical expressions for the probability of collision assume uniformly distributed mean anomalies of the bodies. The initial conditions of the propagation start from a well defined mean anomaly of the NEOs. Thus, we need to track the evolution of the uncertainty in mean anomaly to know when we can start using the analytical estimates. Using our Monte Carlo experiments we compute the standard deviation in mean anomaly separation from the nominal trajectory.

We propagate the orbits of the km-sized NEO population for 1000 yr and study when the MOID is smaller than a Lunar Distance (1 LD). When the standard deviation in mean anomaly is large, we estimate the probability of close encounters. In Figure 6 we list the km-sized NEOs showing the dates in which we found a low-MOID and sorted by their estimated probability of close encounters. The low-MOID regions are color coded with the standard deviation in mean anomaly. The combination of this information highlights the future periods of time in which the position of the NEOs is unknown.

Zoom In Zoom Out Reset image size

Among the 40 km-sized NEOs currently with an Earth MOID < 0.01 au, we find that their mean anomalies remain well defined typically for at least 200 yr, and in some cases for thousands of years. On the other hand, there are a few examples of growth in mean anomaly uncertainty after 2200, such as 35396 (1997 XF11). Because the MOID becomes greater than 1 LD by the time the uncertainty in mean anomaly is large, the estimated probability for this NEO is zero.

The objects with the highest estimated probability are shown in Figure 6 and listed in Table 2. The asteroid with the largest estimated probability of a deep close encounter is 7482 (1994 PC1). This result is to be expected as, in Section 3.1 and Figure 5, we show that 7482 (1991 PC1) has a continuous low MOID. In this analysis we find that 7482 spends about 98% of this millennium with an Earth MOID <1 LD. During this unusually lasting MOID <1 LD period the position is well determined until approximately year 2500.

The propagation of the MOID for the other km-sized NEOs on top of the list is shown in Figure 7. We find that either the Earth MOID of these bodies is secularly drifting to zero, or that the current low-MOID oscillates around zero for a longer period of time. The latter case was observed for 7482 (1994 PC1), but additionally 4179 Toutatis (1989 AC) and 314082 Dryope (2005 CZ36) are in similar situations. Figure 7 shows that deep encounters are expected for these bodies, both in a low-MOID format and as the result of the Monte Carlo experiment.

Zoom In Zoom Out Reset image size

The fact that the position is well determined allows us to determine the geometry of the subsequent close encounters, until the uncertainty grows too large. This assumption leaves a brief period of time between a very well constrained position and the range of validity of the uniformly distributed mean anomaly assumption. For this reason we check if there were actually such very close encounters among the low-MOID NEOs that we found earlier. In general, no close encounters within the Lunar Distance were found in the deterministic part of the trajectories.

There are a few exceptions that should be mentioned: 4179 Toutatis (1989 AC), 220839 (2004 VA) (Both in Figure 7), 20236 (1998 BZ7), 214869 (2007 PA8), and 175114 (2004 QQ) experience close encounters right before or right after the date in which S( MA )>10°. In all of these cases, the MOID tends to zero around the dates in which a deep encounter is expected. In some cases, uncertainty in the position grows largely due to preceding close encounters. In general, we find that the Monte Carlo experiment agrees in finding deep encounters. Thus, the current method is successful in identifying their potential for very deep encounters within the next millennium.

7482 (1994 PC1) has been highlighted in every section of this work because of its remarkable MOID evolution. Its Earth MOID is currently 6.09 × 10⁻⁴ au (0.237 LD), has been near zero for centuries and will remain very low for at least another 1000 yr as shown in Section 3.1. This condition is the reason why it is ranked as the most hazardous NEO in the list of Table 2.

The orbit elements of 7482 (1994 PC1) are shown for reference in Appendix. During the period in which the Earth MOID remains small there are close encounters that cause significant variations semimajor axis and eccentricity. However, arguments of node and perihelion follow a secular drift. In its current orbit within the inner solar system, there is not a large amplitude of short-period oscillations that could disperse the distributions further. However, it is important to highlight that after 500 yr the mean anomalies become uncertain.

Figure 8 shows the sequence of close encounters that are recorded in the Monte Carlo numerical propagation. It appears that the uncertainty in the encounter of 2525 is large enough that the range of possible closest approach distances is between 0 and 0.04 au. Right after the 2525 encounter the standard deviation in mean anomaly increases beyond 10°, and we start estimating its probability of collision using the methods of Section 2.3. Encounters below the Lunar Distance were found after this period, which is consistent with the higher probability that we previously estimated.

Zoom In Zoom Out Reset image size

The km-sized NEO with the shortest deterministic horizon is 143651 (2003 QO104), which was also previously introduced in Figure 6. The orbit solution of 143651 (2003 QO104) has an observation arc of decades, including light-curve observations (Birtwhistle 2009) and radar astrometry (Warner et al. 2009). Thus, we believe that the rapid increase in uncertainty is a dynamical effect of its orbit.

Among the list of NEOs in Table 2 with nonzero estimated probability of having encounters below the Lunar Distance, 143651 (2003 QO104) has the slowest close encounters. These relative velocities imply larger scatter during close encounters, including a rapid increase in mean anomaly uncertainty. As shown in Figure 9, there is a close encounter in 2220 after which the sequence of encounters becomes unique for each Monte Carlo run.

Zoom In Zoom Out Reset image size

Figure 14 shows the evolution of the orbital elements. By the end of the millennium there is a wide variety of orbits in which 143651 (2003 QO104), product of an undetermined sequence of both close and slow Earth close encounters.

The binary asteroid 66391 Moshup (1999 KW4) has been the object multiple studies relative to its binary system condition. It consists of a primary and secondary of respectively 1.317 and 0.59 km of diameter (Ostro et al. 2006; Scheirich et al. 2021). Its rotation states suggest that it is a product of YORP spin-up and disruption (Scheeres et al. 2006; Davis & Scheeres 2020), and its orbit is expanding in time due to the BYORP effect (Scheirich et al. 2021).

The heliocentric orbit of 66391 Moshup (1999 KW4) is in resonance with the Earth, as it experiences resonant close encounters every 17 or 18 yr, as shown in Figure 10. The apparition of 2019 allowed observations from multiple observatories, during the 0.0346au encounter (Scheirich et al. 2021). The next close encounter will be in 2036 May, with a closest distance of 0.0155au, much closer than the first radar observations obtained using the Goldstone and Arecibo radar systems in May of 2001 (Ostro et al. 2006). When the MOID becomes small, which is expected to happen slightly before year 3000, many close encounters below the Lunar Distance are recorded in our Monte Carlo analysis. These will cause a large scattering of the orbit as shown in Figure 15. Because of how relatively late in the millennium the MOID < 1LD condition is held, 66391 Moshup (1999 KW4) is not ranked higher in the list of Table 2.

Asteroid 29075 (1950 AD) is representative example of impact probability studies. After it was discovered and tracked for 17 days in 1950, it was lost for 50 yr until rediscovered on 2000 December 31. Giorgini et al. (2002) found a close approach in 2880 with the possibility of an impact. Farnocchia & Chesley (2014) modeled the Yarkovsky effect on 29075 (1950 AD) and estimated an impact probability of 2.5 × 10⁻⁴.

Zoom In Zoom Out Reset image size

The example of 29075 (1950 AD) is paradigmatic in MOID evolution of NEOs. As shown in Figure 2, its Earth MOID is secularly drifting to zero. During the decades that this condition is maintained, the probability of experiencing a deep encounter is nonzero. 29075 (1950 AD) was found among the list of km-sized NEOs for which we estimated this probability. As we show in Figure 11, the encounters of 2860 and 2880 will occur although with an uncertain closest approach distance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The km-sized 2022 AP7 is one of the largest PHAs recently discovered (Sheppard et al. 2022). The orbit of 2022 AP7 is in near resonance with the orbit of the Earth, meaning that even if its MOID will become small in the next hundreds of years, almost no close encounters are expected in this period of time. The only likely exception is a close encounter in 2363 which will probably to be at a closest approach distance larger than 0.05 au. An interesting finding is that 2022 AP7 comes from a sequence of resonant encounters every 5 yr during the 19th century.