The Manifold Of Variations: impact location of short-term impactors

The interest in the problem of small asteroids observed shortly before a deep close approach or an impact with the Earth has grown a lot in recent years. Since the observational dataset of such objects is very limited, they deserve dedicated orbit determination and hazard assessment methods. The currently available systems are based on the systematic ranging, a technique providing a 2-dimensional manifold of orbits compatible with the observations, the so-called Manifold Of Variations. In this paper we first review the Manifold Of Variations method, to then show how this set of virtual asteroids can be used to predict the impact location of short-term impactors, and compare the results with those of already existent methods.


Introduction
Astronomers observe the sky every night to search for new asteroids or for already known objects. The Minor Planet Center (MPC) collects the observations coming from all over the world and then tries to compute orbits and to determine the nature of the observed objects. New discoveries which could be Near-Earth Asteroids (NEAs) are posted in the NEO Confirmation Page (NEOCP 1 ), which thus contains observational data of unconfirmed objects. They could be real asteroids as well as non-real objects, and cannot be officially designated until additional observations are enough to compute a reliable orbit and confirm the discovery.
Some asteroids with an Earth-crossing orbit may impact our planet and the goal of impact monitoring is to identify potentially hazardous cases and solicit follow-up observations. Two automated and independent systems continually scan the catalogue of known NEAs with this purpose, namely clomon-2 2 and Sentry 3 , which are respectively operational at the University of Pisa/SpaceDyS and at the NASA JPL. Nevertheless, also objects waiting for confirmation in the NEOCP could be on a collision trajectory with the Earth, sometimes with the impact occurring just a few hours after the discovery. This is exactly what happened for asteroids 2008 TC 3 , 2014 AA, 2018 LA, and 2019 MO, all discovered less than one day prior to striking the Earth and prior to being officially designated by the MPC. Thus being able to perform a reliable hazard assessment also in these cases is a fundamental issue, but needs dedicated techniques due to the very different nature of the problem. Indeed when an asteroid is first observed, the available astrometric observations are few and cover a short time interval. This amount of information is usually not enough to allow the determination of a well-constrained orbit and in fact the orbit determination process shows some kind of degeneracies. The few observations constrain the position and motion of the object on the celestial sphere, but leave practically unknown the topocentric range and range-rate. As a consequence the set of orbits compatible with the observations forms a region in the orbital elements space which has a two-dimensional structure, so that every one-dimensional representation of the region such as the Line Of Variations (LOV, Milani et al. (2005b)) would be unreliable.
Two systems are now publicly operational and dedicated to the orbit determination and hazard assessment of unconfirmed objects: Scout at the NASA JPL (Farnocchia et al., 2015) and NEOScan at SpaceDyS (Spoto et al., 2018). They are both based on the systematic ranging, an orbit determination method which explores a subset of admissible values for the range and range-rate. NEOScan 4 makes use of the Admissible Region (AR, Milani et al. (2004)) as a starting point to explore the range and range-rate space. Then the short arc orbit determination process ends with the computation of the Manifold Of Variations (MOV, Tommei (2006)), a 2-dimensional compact manifold of orbits parameterized over the AR. A finite sampling of the MOV is thus a suitable representation of the confidence region, because it accounts for its two-dimensional structure, and thus the resulting set of virtual asteroids can be used for the short-term hazard assessment of such objects (Spoto et al., 2018;. A part of this activity is the prediction of the impact location of a potential impactor, especially when the associated impact probability is high. A method to predict the impact corridor of an asteroid has been developed in Dimare et al. (2020), by which the impact region is given by semilinear boundaries on the impact surface at a given altitude above the Earth and corresponding to different confidence levels. The algorithm is conceived to be a continuation of the impact monitoring algorithm at the basis of the clomon-2 system, since the semilinear method requires a nominal orbit obtained by full differential corrections and an impacting orbit 5 , as provided by the LOV method (Milani et al., 2005a). When the observational arc is very short and the object is waiting for confirmation in the NEOCP, the semilinear method could become inapplicable because very often a reliable nominal orbit simply does not exist. In this paper we propose a method which uses the MOV to predict the impact location of an imminent impactor, and we then test this technique using the data available for the four impacted asteroids so far, namely 2008 TC 3 , 2014 AA, 2018 LA, and 2019 MO.
In Section 2 we give a brief recap about the AR and the MOV method. In Section 3 we introduce the impact surface at a given height over the Earth and the impact map, to then describe how to exploit the MOV sampling for impact location predictions. Section 4 contains the results of our method applied to the impacted asteroids 2008 TC 3 , 2014 AA, 2018 LA, and 2019 MO. When possible, we also compare the impact regions with those computed with the semilinear method (Dimare et al., 2020) and with a Monte Carlo simulation. Lastly, Section 5 is dedicated to the conclusions.

The Manifold Of Variations method
Suppose we have a short arc of observations, typically 3 to 5 observations over a time span shorter than one hour. In most cases the arc is too short to allow the determination of a full orbit and it is called a Too Short Arc (TSA, Milani et al. (2007)). When in presence of a TSA, either preliminary orbit methods fail or the differential corrections do not converge to a nominal orbit. Nevertheless, as anticipated in the introduction, the little information contained in the arc can be summarised in the attributable, the vector A := (α, δ,α,δ) ∈ S 1 × − π 2 , π 2 × R 2 formed by the angular position and velocity of the object at the mean observational time. Note that if there is at least one measurement of the apparent magnitude, the attributable could also contain an average value for this quantity. The topocentric range ρ and range-rateρ are thus not known, otherwise we would have had a full orbit. The AR comes into play here, to provide a set of possible values of (ρ,ρ) by imposing general conditions on the object's orbit. It can be shown that the AR is a compact subset of R ≥0 × R, which can have at most two connected components. For the precise definition of the AR and the proof of its properties the reader can refer to Milani et al. (2004), here we limit ourselves to the general idea. We essentially impose that the object is a Solar System body and that it is sufficiently large to be source of meteorites. Indeed short-term impactors are usually very small asteroids, with a few meters in diameter, thus the main interest in such objects is not for planetary defence, but rather to observe them during the last part of their impact trajectory and possibly to recover meteorites on ground, as it happened for asteroid 2008 TC 3 .
The AR is sampled by a finite number of points and with different techniques. In case of a TSA we explore the whole AR with rectangular grids: first we cover the entire region and compute the corresponding sample of orbits, then we identify the subregion corresponding to the best-fitting orbits with respect to the observations, and lastly we cover this subregion with a second grid. The two-step procedure is shown in Figure 1. Even if this is not the most common case, it may happen that the short arc of observations allows the computation of a preliminary orbit and even of a full orbit. In this case we consider the nominal solution as reliable only if the arc curvature is significant (Milani et al., 2007). Given its importance for orbit determination purposes, we give more weight to the geodesic curvature with respect to the along-track acceleration, imposing that the signal-to-noise ratio of the geodesic curvature is greater than 3. In this case we switch to a different sampling method to exploit the additional strong information coming with the reliable nominal orbit. Indeed we consider the marginal covariance of the orbit in the range and range-rate space, whose probability density function has level curves which are concentric ellipses around the nominal range and range-rate. We select a maximum confidence threshold and sample the AR by following these ellipses up to this confidence level. This cobweb technique is shown in Figure 2. Full details about the sampling of the AR in the various cases can be found in Spoto et al. (2018). Left. First step, with the rectangular grid covering the entire AR. Right. Second grid, covering the subregion of good orbits identified in the first step. In both plots the sample points are marked in blue when χ ≤ 2 and in green when 2 < χ < 5. The orange cross marks the orbit with the minimum χ value. Now we describe how to obtain a sampling of orbits, namely the MOV, once a sampling of the AR is given. We start with some notation. The target function is where x are the orbital elements, m is the number of observations used for the fit, and ξ is the vector of the normalised observed-computed debiased astrometric residuals. Let A 0 be the attributable computed from the arc of observations. The AR sampling is a set K of admissible values of the range and range-rate. For each ρ 0 = (ρ 0 ,ρ 0 ) ∈ K we consider the full orbit (A 0 , ρ 0 ,ρ 0 ) and fit only the first four components, that is the attributable part, with a constrained differential corrections procedure.
We call K ⊆ K the set of values of (ρ,ρ) such that the constrained differential corrections converge, giving an orbit on M. To estimate the goodness of the fit to the observations, for each orbit x ∈ M we compute where Q * is the minimum value of the target function over K . We always consider MOV orbits having χ < 5, thus corresponding at most to the 5σ confidence level. For hazard assessment, this guarantees to find impact possibilities with a probability > 10 −3 , the so-called completeness level of the impact monitoring system (Del Vigna et al., 2019a).
Hereafter we summarise the main steps of the Manifold Of Variations method, as it is implemented in the software system NEOScan.
(i) The MPC NEO Confirmation Page is scanned every two minutes to look for new objects or for new observations to add to previous detections.
(ii) Computation of the attributable and sampling of the Admissible Region with rectangular grids or with the cobweb techinque, depending on the existence of a reliable nominal solution.
(iii) Computation of the Manifold Of Variations by constrained differential corrections, obtaining a set of virtual asteroids with a two-dimensional structure.
(iv) Propagation of the virtual asteroids for 30 days and projection of the propagated MOV sampling on the Modified Target Plane 6 .
(v) Searching for impacting solutions and, if any, computation of the impact probability. The computation of this quantity is done with a propagation of the probability density function from the normalised residuals space to the sampling space and by integrating the resulting density over the set of impacting sample points (Del Vigna, 2020).

Impact location prediction
Before describing how to exploit the MOV sampling for impact predictions, we emphasise a key difference in the treatment of short-term impactors with respect to the long-term hazard monitoring of known asteroids. The MOV method does not use interpolation of the sample to find impacting orbits as it is done for the LOV method (Milani et al., 2005a), but just considers the impacting orbits of the sampling. For instance, if none of the MOV sample points impacts the Earth, we assign a null impact probability. Indeed we adopt a dense sampling, thus if a set of impacting orbits is not detected by the sampling it has a too small probability to be interesting. This is a crucial difference with respect to the LOV method at the basis of long-term impact monitoring. Many of the NEAs in the catalogue are big objects and thus also impact events with probability of few parts per million are worth detecting. Since nearby LOV orbits are separated by the chaotic dynamics introduced by close approaches, the prediction of impacts with such low probabilities, especially if far in time with respect to the time of the observations, would require an extremely high number of sample points along the LOV. Hence if we limited the analysis to the sample points only, the computations would be too heavy. As a consequence, interpolation of the LOV sampling is essential 7 , but with the MOV method we are allowed to consider just the impacting orbits of the sampling.
On the basis of the above comment, the short-term hazard analysis performed with the MOV method ends with the identification of a subset V of impacting orbits among the MOV sampling (see step (v) at the end of the previous section). The idea of our location prediction method is actually very simple: once the set V is given, we just propagate the orbits of V until the impact and compute the geodetic coordinates of the impacting points.
For such predictions we consider the WGS 84 model (NIMA -National Imagery and Mapping Agency, 2000), for which the Earth surface is approximated by a geocentric oblate ellipsoid with semimajor axis equal to 6378.137 km and flattening parameter f defined by 1/f = 298.257223563. The eccentricity e of the ellipsoid can be computed as e 2 = f (2 − f ).
Definition 3.1. Let h ≥ 0. The impact surface S h at altitude h above ground is the set of points at altitude h above the WGS 84 Earth ellipsoid.
The impact region will be a subset of the impact surface S h for a given value of the altitude h and points on S h are given by means of the geodetic coordinates. To compute this region we use the impact map introduced in Dimare et al. (2020), that is This map takes an impacting orbit x ∈ V, computes the time t(x) at which the orbit reaches the surface S h , applies the integral flow of the system to propagate the orbit to the time t(x), and lastly converts the state vector into the geodetic coordinates on S h . Thus our MOV impact region is the set F h (V) ⊆ S h , that is the set of impacting MOV orbits propagated to the surface S h .
3.1. Semilinear boundaries for the impact corridor prediction. We briefly recap the semilinear method applied to the problem of the impact location prediction, described in Dimare et al. (2020).
Let X be the orbital elements space, let Y ⊆ R 2 be the target space, and let F : W → Y be a differentiable function defined on an open subset W ⊆ X , to be thought of as the prediction function. The dimension 7 Actually in some extremely non-linear cases this is not enough, and LOV sampling densification techniques have been developed to overcome the problem and to guarantee that the preset completeness level is reached .
N of X is 6 if we consider the six orbital elements, or is > 6 if some dynamical parameter is included 8 . Furthermore we consider a nominal orbit x * ∈ X provided with its N × N covariance matrix Γ X . The uncertainty of x * can be described through the confidence region where σ > 0 is the confidence parameter. The prediction function F maps the orbital elements space onto the target space, and thus maps the confidence region Z X (σ) to Z Y (σ) := F (Z X (σ)). The semilinear method provides an approximation of the boundary ∂Z Y (σ) of the non-linear prediction region (Milani and Valsecchi, 1999).
To explain the construction of the semilinear boundaries we start with the notion of linear prediction. The inverse of Γ X , that is C X = Γ −1 X , is the normal matrix and defines the linear confidence ellipsoid according to the covariance propagation law. When F is strongly non-linear, this linear prediction region is a poor approximation of Z Y (σ) and here the semilinear method comes into play. The boundary ellipse ∂Z Y lin (σ) is the image by the linear map DF x * of an ellipse E X (σ) ⊆ X which lies on the surface ∂Z X lin (σ). The semilinear boundary K(σ) is the non-linear image in the target space of the ellipse E X (σ), that is If the curve K(σ) is simple (no self-intersections) and closed, then by the Jordan curve theorem it is the boundary of a connected subset Z(σ) ⊆ Y, which we use as an approximation of Z Y (σ).
Remark 3.2. The differential DF x0 is clearly not injective. More precisely, each point of Y has a (N − 2)dimensional fiber in X , and thus in principle the selection of the ellipse E X (σ) could be done in infinitely many ways. The choice foreseen in the semilinear method consists in selecting as E X (σ) the ellipse resulting from the intersection between a suitable regression subspace in X and the confidence ellipsoid Z X lin (σ). See Milani and Valsecchi (1999) for all the details.
The application of this method to the impact corridor problem is described in Dimare et al. (2020) and works as follows. We have a nominal orbit x 0 of an asteroid with a non-negligible chance of impacting the Earth, and thus a virtual impactor representative x imp as provided by the LOV method. In this problem the prediction map is the impact map F h : W → S h defined as above, where W ⊆ X is an open neighbourhood of x imp . The semilinear method can now be applied, with the following adaptation. For the linear prediction on the impact surface and the selection of the ellipse E X (σ) we use the differential (DF h ) ximp of the impact map at the virtual impactor representative orbit x imp . The result of this method are curves on the surface S h corresponding to different values of σ. Note that these curves need not be closed in this context, because in general not all the orbits on E X (σ) impact.  Tholen and Farnocchia (2020) for asteroid Apophis, in both cases succeeding in ruling out the most relevant virtual impactors.

Numerical tests
We test the MOV method for the impact location prediction of four asteroids, namely 2008 TC 3 , 2014 AA, 2018 LA and 2019 MO. When we have a reliable nominal solution we also compare these results with those of the semilinear method and, additionally, to an observational Monte Carlo simulation. The latter amounts to adding noise to each observation to then compute a new orbit based upon the revised observations and use this orbit as a virtual asteroid. When the observed arc is short the uncertainty in the orbit determination is typically so large that the true uncertainty is not represented by the confidence ellipsoid. In this case the observational Monte Carlo is the best approach to follow, definitely better than an orbital Monte Carlo, which works well when the orbital uncertainty is fairly good. 4.1. Asteroid 2008 TC 3 . The small asteroid 2008 TC 3 was discovered by Richard A. Kowalski at the Catalina Sky Survey on October 6, 2008. The preliminary orbit computations done at the MPC immediately revealed that the object was going to impact the Earth within 21 hours. Thus the astronomical community made a great effort to observe it and the currently available dataset contains nearly 900 observations, though not all of them are of high quality, and need to be properly treated for a precise estimate of the trajectory of 2008 TC 3 (Farnocchia et al., 2017).
When asteroid 2008 TC 3 was recognised to be an impactor it only had few observations and for our analysis we consider its first seven observations. They are enough to compute a reliable nominal solution, so that our algorithm samples the AR with the cobweb technique. The application of the MOV method results in an impact probability of 99.7%, which means that the vast majority of the MOV orbits impacts the Earth. We then propagate the MOV sampling to the Earth surface, i.e. we set h = 0, and obtain the result shown in Figure 3. Since the impact region is very extended, we limit our analysis to the MOV orbits with χ < 3.
In the left figure we show the geodetic coordinates of the impacting points, with different shades of gray depending on the χ value. In the right figure we show the results of a Monte Carlo simulation with 10,000 sample points, to compare with the MOV impact region. Moreover, thanks to the existence of a nominal orbit we can also apply the semilinear method, thus in both plots of Figure 3 we also draw the semilinear boundaries on ground corresponding to the confidence levels 1, 2, and 3. The result is a strong agreement of the three methods.
We notice that the semilinear boundaries enclose a region which is larger than that obtained with the MOV method. This happens also in most of the other examples which we present in the next sections. We can explain this behaviour as follows. Recall that the basic idea of the semilinear method is to select a curve in the orbital elements space to propagate non-linearly to approximate the boundary of the nonlinear prediction region. The choice of this curve is made by considering the boundary of the marginal uncertainty on a suitable space, which is equivalent to considering the corresponding regression subspace (see Remark 3.2). The marginal covariance is the largest projected uncertainty on a given space thus, so to say, is the most conservative choice. 4.2. Asteroid 2014 AA. Asteroid 2014 AA was discovered by R. Kowalski at the Catalina Sky Survey on January 1, 2014 at 06:18 UTC. Similarly to 2008 TC 3 , also 2014 AA impacted the Earth just 21 hours after its first detection. But very differently from 2008 TC 3 , due to the particular night in which 2014 AA was spotted, it was not recognised to be an impactor. As a consequence the astrometric dataset is very limited, containing just 7 observations. Also in this case we can fit a nominal orbit and thus the MOV method starts with the cobweb sampling of the AR. All the MOV orbits are impacting, resulting in an impact probability of 100%. The impact regions on ground of 2014 AA computed with the MOV or with the Monte Carlo method are shown in Figure 4.   Figure 3. Impact location prediction on ground for asteroid 2008 TC 3 . Left. Comparison between the impact region computed with the MOV orbits having χ < 3 and the semilinear boundaries corresponding to the confidence levels 1, 2, and 3. Right. Comparison between a Monte Carlo run and the same semilinear boundaries.
Again we also show the semilinear boundaries for comparison, this time corresponding to the confidence levels σ = 1, 3, and 5.  Figure 4. Impact region on ground of asteroid 2014 AA computed with the MOV method and by the semilinear method. Left. Comparison between the impact region computed with the MOV orbits having χ < 5 and the semilinear boundaries corresponding to the confidence levels 1, 3, and 5. Right. Comparison between a Monte Carlo run and the same semilinear boundaries.
The plot deserves some comment to avoid misunderstadings, due to the particular shape of the impact region. It is known that semilinear boundaries are not necessarily simple curves. As a consequence, if this is the case, the curve cannot be the boundary of the non-linear prediction and indeed the impact regions extend outside the drawn boundaries. This behaviour is confirmed by the fact that, in the vicinity of the torsion, the boundaries with lower σ extends outside the ones with higher σ. With a heuristic approach we can state that the impact region extends outside the region delimited by the boundary, on the side of the lower confidence level boundaries. As a consequence this result is somewhat unsatisfactory, because we can just guess the actual impact area. Of course this issue disappears as soon as we start with a sampling of the whole confidence region and not with just the sampling of a curve. As we can see, the impact region computed with the MOV gives a clear representation of the impact area and also shows why the semilinear boundaries are twisted. Indeed the confidence region folds on itself before being projected on the impact surface S h , and this is a consequence of the non-linear effects due to the ongoing close approach. This behaviour is also confiermed by the Monte Carlo simulation shown in the right plot of Figure 4.  The astrometric dataset contains 14 observations, which are enough to compute a full orbit and to apply the cobweb sampling. The impact probability computed with the MOV method gives the certainty of the impact and all the MOV orbits impact the Earth. Figure 5 shows the MOV impact region at h = 28.7 km, together with the fireball event. Figure 6 shows the comparison between the MOV impact region or the Monte Carlo impact region against the semilinear buondaries computed at the height of the fireball event. Also in this case we have full agreement among the methods. As pointed out in the introduction, the MOV method can be applied for impact location predictions also when the amount of information provided by the observational arc is not enough to constrain a full orbit. In this case the semilinear method of Dimare et al. (2020) cannot be used, because it requires a nominal orbit and a virtual impactor representative, as provided through the LOV method. To show an example of such computation, we applied the MOV method to 2018 LA with its first 12 observations. In this case the observations are not enough to compute a reliable nominal orbit, thus the AR is sampled with rectangular grids. The second grid samples a subregion of the AR, with the orbits having χ < 5 spanning a quite narrow strip (see Figure 7). The propagation of the impacting MOV orbits gives the result shown in Figure 8, with the impact region being quite elongated due to the poor constraint on the asteroid orbit.
4.4. Asteroid 2019 MO. Asteroid 2019 MO was discovered from the ATLAS Mauna Loa observatory on June 22, 2019, less than 12 hours before impacting the Earth between Jamaica and the south American coast. Also in this case we have a fireball event reported in Table 1, which took place at 25 km of altitude.  Figure 6. Impact location prediction at h = 28.7 km for asteroid 2018 LA, with all the 14 observations. Left. Comparison between the impact region computed with the MOV orbits having χ < 5 and the semilinear boundaries corresponding to the confidence levels 1, 3, and 5. Right. Comparison between a Monte Carlo run and the same semilinear boundaries. The star marks the location of the fireball event. Figure 7. Second step of the Admissible Region grid sampling for asteroid 2018 LA with 12 observations. The sample points are marked in blue when χ ≤ 2 and in green when 2 < χ < 5 (note that the blue points are contained in the narrow region delimited by the green points). The red circles mark the impacting orbits. Figure 9 shows the MOV impact region at h = 25 km extending above the south American coast, together with the fireball location. Figure 10 shows the comparison between the MOV method, the semilinear boundaries and a Monte Carlo run with 10,000 sample points, which are fully compatible. 2018LA χ <1 1 < χ <3 3 < χ <5 Fireball location Figure 8. Impact region at altitude h = 28.7 km of asteroid 2018 LA computed with the MOV method and using the first 12 observations. The star marks the location of the fireball event. -

Conclusions
Very small asteroids can be only observed during a deep close approach with the Earth and it may be the case that an impact occurs a few days after the discovery. In this paper we considered the problem of predicting the impact location of such objects by exploiting the MOV, a set of virtual asteroids representing the orbital uncertainty. Once a MOV sampling is available it suffices to propagate each orbit for a given  Comparison between the impact region computed with the MOV orbits having χ < 5 and the semilinear boundaries corresponding to the confidence levels 1, 3, and 5. Right. Comparison between a Monte Carlo run and the same semilinear boundaries. The star marks the location of the fireball event.
amount of time and, for the impacting orbits, to compute the geodetic coordinates of the impacting points. The impact region is thus given by a set of points on the impact surface at a certain height over the Earth. The advantage of the MOV method with respect to already existing techniques, like for example the semilinear method of Dimare et al. (2020), is that it can be used even when it is impossible to fit a nominal orbit to the few available astrometric observations (see the example presented in Section 4.3 and shown in Figure 8).
We tested our method using the data available for the four impacted asteroids so far, namely 2008TC 3 , 2014AA, 2018LA, and 2019MO. Since these data are also enough to constrain a full orbit, we compared the results of the MOV method with the semilinear boundaries and with the outcome of an observational Monte Carlo simultation, getting a very strong consistency among the three methods.