Short arc orbit determination and imminent impactors in the Gaia era

Short-arc orbit determination is crucial when an asteroid is first discovered. In these cases usually the observations are so few that the differential correction procedure may not converge. We have developed an initial orbit computation method, based on the systematic ranging, an orbit determination techniques which systematically explores a raster in the topocentric range and range-rate space region inside the admissible region. We obtain a fully rigorous computation of the probability for the asteroid that could impact the Earth within few days from the discovery, without any a priori assumption. We test our method on the two past impactors 2008 TC3 and 2014 AA, on some very well known cases, and on two particular objects observed by the ESA Gaia mission.


Introduction
Short-arc orbit determination is a very important step when an asteroid is first discovered. In these cases the timing is essential, because we are interested in a rapid follow-up of a possible imminent impactor, which is an asteroid impacting the Earth shortly after its discovery, within the same apparition (interval of observability). The observations are so few that the standard differential correction procedure (Milani & Gronchi 2010) to find an orbit by a least-squares minimization fails, and other methods need to be used to extract information on the orbit of the object.
Several initial orbit computation methods have been developed in the last 25 years. For instance, Muinonen & Bowell (1993) define a Gaussian probability density on the orbital elements space using the Bayesian inversion theory. In particular, they determine asteroid orbital elements from optical astrometric observations using both a priori and a posteriori densities, the latter computed with a Monte Carlo method.
The few observations in the short arc constrain the position of the object in the sky, but they leave almost unknown the distance from the observer (topocentric range) and the radial velocity (topocentric range-rate). Thus, ranging methods have been developed over the years to replace or refine the Monte Carlo approach in the short arc orbit determination. There are two alternative approaches to the ranging methods: the statistical and the systematic ones. The original statistical ranging method , Muinonen et al. (2001)) starts from the selection of a pair of astrometric observations. Then, the the topocentric ranges at the epoch of the observations are randomly sampled. Candidate orbital elements are included in the sample of accepted elements if the χ 2 value between the observed and computed observations is within a pre-defined threshold. Oszkiewicz et al. (2009) improve the statistical ranging using the Markov-Chain Monte Carlo (MCMC) to sample the phase space. The MCMC orbital ranging method is based on a bi-variate Gaussian proposal PDF for the topocentric ranges. Then, Muinonen et al. (2016) develop a random-walk ranging method in which the orbital-element space is uniformly sampled, up to a χ 2 value, with the use of the MCMC method. The weights of each set of orbital elements are based on a posteriori probability density value and the MCMC rejection rate. They have developed this method for the ESA Gaia mission, in the framework of Gaia alerts on potentially new discovered objects by Gaia (see Tanga et al. (2016)).
On the other side, Chesley (2005) and Farnocchia et al. (2015b) introduce the so-called systematic ranging, which systematically explores a raster in the topocentric range and range-rate space (ρ, 9 ρ). This technique permits to describe the asteroid orbital elements as a function of range and range-rate. Then, the systematic ranging also allows one to determine the subset of the sampling orbits which leads to an impact with the Earth.
In this paper we describe a new approach to the systematic ranging, based on the knowledge of the Admissible Region (AR) (Milani et al. 2004), and a new method to scan the re-gion. The process has some main advantages if compared to the other methods described above.
1. Our grid is more efficient, for two main reasons: -we discard all the objects that are not in the AR, saving CPU time and making the systematic ranging more accurate in finding the region in the (ρ, 9 ρ) space of the possible orbital solutions; -we use two different grids depending on the boundary of the AR. The first grid is larger and less dense, the second one is based on a refinement using the value of the post-fit χ 2 of each point in the first grid (see Section 2.1).
2. The computation of the probability for the potential impactors is a rigorous probability propagation from the astrometric error model, without any assumption of a priori probability density function on the range/range-rate space (see Section 3).
In Section 2.1 we summarize the features of the AR, and describe the sampling of the (ρ, 9 ρ) space needed for the creation of the Manifold Of Variations (MOV), an extension of the Line Of Variations (LOV, Milani et al. (2005b); Milani & Gronchi (2010)). In Section 2.2 we present the so-called spider-web (Tommei 2006), to be used when a nominal orbit is available. Section 3 describes how to compute the impact probability when an impactor is found. In Section 4 we apply our method to the well known cases 2008 TC 3 and 2014 AA, and we show other examples to explore the capabilities of the new method. Then in Section 5 we test our method on new objects discovered by Gaia in the framework of the Gaia alerts, and we show the points of strength of this approach applied to Gaia observations. Section 7 contains our conclusions and comments.
2 Sampling of the topocentric range and range-rate space

Admissible Region and systematic ranging
Even if the observations are too scarce, we can anyway compute the right ascension α, the declination δ, and their time derivatives 9 α and 9 δ, by fitting both angular coordinates as a function of time with a polynomial model. This four quantities can be assembled together to form the attributable (Milani & Knežević 2005): at a chosen timet, which could be the time of the first observation or the mean of the observation times. The information contained in the attributable leaves completely unknown the topocentric distance ρ and the radial velocity 9 ρ. We would have a full description of the topocentric position and velocity of the asteroid in the attributable elements (α, δ, 9 α, 9 δ, ρ, 9 ρ), if ρ and 9 ρ were known. From now on, we will use x = (A, ρ), where ρ = (ρ, 9 ρ), to describe the attributable elements.
Given an attributable, we can define the AR as the set of all the possible couples (ρ, 9 ρ) satisfying the following conditions (see Milani et al. (2004) for more mathematical details): 1. the object belongs to the Solar System, and it is not a too long period comet. We consider only the objects for which the value of the heliocentric energy is less than −k 2 /(2a max ), where a max = 100 au and k = 0.01720209895 is the Gauss' constant; 2. the corresponding object is not a satellite of the Earth, i.e. the orbit of the object has a non-negative geocentric energy when inside the sphere of influence of the Earth, whose radius is The AR is a compact set, and it can have at most two connected components, which means that it could be represented as the union of no more than two disjoint regions in the (ρ, 9 ρ) space. The AR has usually one component, and the case with two components indicates the possibility for the object to be distant (perihelion q > 28 au). As shown in Milani et al. (2004), the number of connected components depends on the number of the roots of a polynomial resulting from condition number 1. This polynomial cannot have more than three distinct real positive roots: the AR has two connected components if the roots are three, and it has one component if there is only one root. It is worth noting that the region defined by condition number 1 could contain points with arbitrarily small values of ρ. The boundary of the region given by condition 2 turns out to have two different shapes: it can be formed just by the line of geocentric energy equal to 0 (if it is entirely contained in the region 0 < ρ < R S I ), or by a segment of the straight vertical line ρ = R S I and two arcs of the zero-curve of the geocentric energy (for 0 < ρ < R S I ). We also discard the orbits corresponding to meteors too small to be source of meteorites, using the condition H ≤ H max , where H max = 34.5 is the shooting star limit (Milani et al. 2004), and H is the absolute magnitude. Given all these conditions, we can sample the AR with a finite number of points.
If a nominal solution does not exist, we make use of the systematic ranging to scan the AR. We sample the AR in two different ways, depending on the number of connected components, and the values of the roots (r 1 , r 2 and r 3 , in ascending order). Table 1 summarizes the conditions and the grids used in the (ρ, 9 ρ) space. In particular, we compute a rectangular grid in the range/range-rate space, with range as in Table 1, and range-rate controlled by the AR equations (Milani et al. 2004). Nevertheless, since the AR has a shape dictated by a polynomial equation and it is not a rectangle, we check the value of the heliocentric energy for each grid point, and we discard those not satisfying condition 1. Orbits not satisfying condition 2 are discarded as well, except when we compute the probability for the asteroid to be a satellite of the Earth 1 . ρ) space with respect to the values of the roots, and the connected components of the AR. The sampling in 9 ρ is always uniform.
Roots AR Grid Sampling in ρ components r 1 < √ 10 au 1 50 × 50 Unif. in log 10 (ρ) The target function is defined by where x = (A, ρ) are the fit parameters, m is the number of observations used in the least squares fit, ξ is the vector of the observed-computed debiased astrometric residuals 2 , and W is the weight matrix. The choice of the weights for each observatory is fundamental, and it has to take into account the debiasing of the star catalog systematic errors, unless the astrometric reduction has already been performed with an essentially bias-free star catalog, e.g. the Gaia DR1 (Lindegren, L. et al. 2016). Given a subset K of the AR, we define the Manifold Of Variations M as the set of the points (A * (ρ 0 ), ρ 0 ) such that ρ 0 ∈ K and A * (ρ 0 ) is the local minimum of the function Q| ρ=ρ 0 . In addition, the value of the minimum RMS of the residuals is less than a given threshold Σ. In general, the MOV is a 2-dimensional manifold, such that the differential of the map from the sampling space to M has rank 2.
In the case of the systematic ranging, K is the AR, scanned with a regular semi-logarithmic or uniform grid. For each sample point ρ 0 = (ρ 0 , 9 ρ 0 ) we fix ρ = ρ 0 and 9 ρ = 9 ρ 0 in the target function, and then we search A * (ρ 0 ) by means of an iterative procedure, the doubly constrained differential corrections. The normal equation is We indicate as K the subset of K on which the doubly constrained differential corrections converges. In this way, the sampling of the MOV is done over K . For each point x on the MOV, we also compute a χ value where Q * is the minimum value of the target function: Q(x * ) if a reliable nominal solution exists, or the minimum value of Q(x) over K otherwise.
When a nominal solution does not exist, the systematic ranging is performed by a two-step procedure. The first iteration is to compute a grid following the conditions in Table 1. Once we obtain a first preliminary grid, we densify for a higher resolution. We select the minimum and the maximum value of ρ and 9 ρ among all the values of the points for which we have a convergence of the 4-dimension differential correction, and the value of χ is less than 5. We also compute the score with respect to the first grid, and we use the value to select the grid for the second step. The score gives us a first insight into the nature of the object, even when the asteroid is not a potential impactor. We define the score as a probability of the object to belong to different classes (NEO, MBO, DO, and SO), where each class is defined by the following conditions: • NEO: Near Earth Object, an object with q < 1.3 au, where q is the perihelion distance (in au).
• MBO: Main Belt Object, belonging either to the Main Belt or to the Jupiter Trojans. In particular we choose the conditions 1.7 au < a < 4.5 au e < 0.4 or 4.5 au < a < 5.5 au e < 0.3 where a is the semi-major axis (in au), and e is the eccentricity.
• SO: Scattered Object, not belonging to any of the previous classes.
If the object is close to be a NEO (the NEO-score is more than 50%), we use a uniform grid in log 10 (ρ), otherwise we use a uniform grid in ρ. Then we compute again the Manifold Of Variations on the new and denser 100 × 100 grid.

Spider web
Let us suppose that a nominal orbit has been obtained by unconstrained differential corrections, starting from a preliminary orbit as first guess (for instance using the Gauss' method, (Milani & Gronchi 2010)). Then, we can use the nominal solution as the center of the subset of the MOV we are interested in, and we can adopt a different procedure to scan the AR. If a nominal orbit exists, and the value of the geodesic curvature signal-to-noise ratio (SNR) (Milani et al. 2008) is greater than 3, instead of using a grid (Section 2.1), we compute a spider web sampling in a neighborhood of the nominal solution (Tommei 2006). This is obtained by following the level curves of the quadratic approximation of the target function used to minimize the RMS of the observational residuals. The advantage of the use of the cobweb is that, firstly, it is faster than the systematic ranging, and secondly it is more accurate in the cases for which we have already a reliable nominal solution.
Let x * be the nominal solution with its uncertainty, represented by the 6 × 6 covariance matrix Γ. In a neighborhood of x * , the target function can be well approximated by means of the quadratic form defined by the normal matrix C = Γ −1 . The matrix C is positive definite, hence the level curves of the target function are concentric 5-dimensional ellipsoids in the 6-dimensional orbital elements space. The level curves on the (ρ, 9 ρ) space are represented by the marginal ellipsoids, defined by the normal matrix C ρρ = Γ −1 ρρ , where Γ ρρ is the restriction of Γ to the (ρ, 9 ρ) space. To sample these curves we choose the maximum value σ max = 5 for the confidence parameter. Then, for each level curve within the confidence level σ max , we select the points corresponding to some fixed directions. We initially create a regular grid of points in the space of polar elliptic coordinates (R, θ), where 0 ≤ θ < 2π and 0 ≤ R ≤ σ max 3 . Then we apply to each point of the grid the following transformation, depending on the covariance matrix of the nominal orbit and on the orbit itself: where λ 1 > λ 2 are the eigenvalues of the 2 × 2 matrix Γ ρρ , v 1 is the eigenvector corresponding to the greatest eigenvalue λ 1 , and ρ * and 9 ρ * are the range and range-rate values of the nominal solution. Figure 1 shows an example of spider web sampling on the (ρ, 9 ρ) plane.
Figure 1: The figure shows an example of spider web around the nominal solution ρ * = (ρ * , 9 ρ * ). The points follow concentric ellipses, corresponding to different values of the parameter R. For each fixed direction (identified by θ), there is one point on each level curve.

Probability density computation
We obtain a probability distribution on the sampling space to be used for several applications, such as the computation of the impact probability or the score. We begin assuming that the residuals are a Gaussian random variable Ξ, with zero mean and covariance Γ ξ = W −1 . Hence the probability density function on the residuals space is 3 We assume that the nominal solution corresponds to the point (0, 0).
A possible approach to propagate the density (4) to the sampling space uses the Bayesian theory to combine the density coming from the residuals with a prior distribution. The a posteriori probability density function for (ρ, 9 ρ) is given in Muinonen & Bowell (1993) as where p prior is a prior distribution on the sampled space.
Hereinafter we report some possible choices for the prior probability.
• Jeffreys' prior. It has been used for the first time in Muinonen et al. (2001). It takes into account the partial derivatives of the vector of the residuals with respect to the coordinates (ρ, 9 ρ). Jeffreys' prior tends to favor orbits where the object is close to the observer, because of the sensitivity of the residuals for small topocentric distances. Granvik et al. (2009) makes versions of the code (which utilizes Jeffreys' prior) publicly available for the first time.
• Prior based on a population model. This approach requires the computation of the prior distribution as posterior of another prior, which is selected as proportional to ρ 2 by geometric considerations on the Cartesian space of position and velocity.
• Uniform distribution. Uniform distribution in the (ρ, 9 ρ) space. Farnocchia et al. (2015b) gives a detailed description of all these different possible choices, and they also analyze how the impact probabilities change according to different prior distributions. They conclude that the uniform distribution is a good choice for an a priori probability density function, because it represents a good compromise between a simple approach and the identification of potential impactors.
Hereinafter we propose a new method to propagate the probability density function p Ξ (ξ) back to the sampling space. This method is a rigorous propagation of the density function according to the probability theory, and it does not use any a priori assumption. The main idea is that we propagate the probability density function from the residual space to the orbital element space (more precisely, on the Manifold of Variations), and then to the sampling space, according to the Gaussian random variable transformation law. The main difference in using this approach instead of the uniform distribution, is that we compute the Jacobian determinant of the transformation and it is not equal to 1, which is the value chosen by Farnocchia et al. (2015b).
We define the following spaces: • S is the space of the sampling variables. It changes depending on the case we are considering: S = R + × R if the sampling is uniform in ρ, S = R 2 if the sampling is uniform in log 10 ρ, and S = R + × S 1 in the cobweb case.
• K is the subset of the points of the AR such that the doubly constrained differential corrections give a point on the MOV; • M is the MOV, a 2-dimensional manifold in the sixdimensional orbital elements space X; • R m is the residuals space. The residuals are a function of the fit parameters: ξ = F(x), with F : X → R m differentiable, and we define the manifold of possible residuals as V = F(X) (Milani & Gronchi 2010, Section 5.7).
Without loss of generality, we can assume that where I m is the m × m identity matrix. As explained in Milani & Gronchi (2010, Section 5.7), this is obtained by using the normalized residuals in place of the true residuals; biases due to star catalog can also be removed while forming the normalized residuals. With this technique, the probability density function becomes normalized to a standard normal distribution. Thus from now on, we will use ξ to indicate the normalized residuals, and the function F maps the orbital elements space to the normalized residuals space. Then we consider the following chain of maps (defined in A) and we use their Jacobian matrices to compute the probability density function on S . Let s be the variable of the sampling space S , and let S be the corresponding random variable. We use the compact notation χ 2 (s) to indicate χ 2 (x(ρ(s))). The probability density function of S is where M µ is the 2 × 2 matrix associated to f µ (considered as tangent map to the surface M), and M σ is the 2 × 2 Jacobian matrix of f σ . The derivation of (6) is given in A (explicit expressions for the Jacobian determinants are provided in A.2 and A.3, respectively). It is worth noting that we limit our analysis to Solar System orbits (condition number 1 of Section 2.1), because interstellar objects are very rare. As a consequence, we use a Bayesian theory with a population limited to the Solar System, and all the probability computations we describe are actually conditional probabilities to the AR, which is even a subset of the Solar System.

Impact probability computation
Each point on the MOV can be thought as an orbit compatible with the observations, and we call it Virtual Asteroid (VA). We propagate the VAs into the future, currently for 30 days from the date of the observations, and we search for Virtual Impactors (VIs), which are connected sets of initial conditions leading to an impact (Milani et al. 2005a). If a VI has been found on the Modified Target Plane (MTP, Milani & Valsecchi (1999)), it is associated to a subset V ⊆ S of the sampling space, and hence its probability is given by If for a given object we find impacting solutions, we assign to the object an impact flag, which is an integer number related to the computation of the impact probability. It depends on the impact probability and on the arc curvature, as shown in Table 2. An arc has significant curvature if χ 2 > 10, where χ is the chi-value of the geodesic curvature and the acceleration (as defined in Milani et al. (2007)). The impact flag can take the integer values from 0 to 4: 0 indicates a negligible chance of collision with the Earth, whereas the maximum value 4 expresses an elevated impact risk (≥ 1%). It is conceived as a simple and direct communication tool to assess the importance of collision predictions, and to give the priority for the follow-up activities.
Table 2: The table shows the conditions on the Impact Probability (IP) and on the arc quality to assign the impact flag to an object.
IP > 10 −2 and no significant curvature 4 IP > 10 −2 and significant curvature The goal for a system dedicated to imminent impactors is to detect all the possible VIs down to a probability level of about 10 −3 , called completeness level (Del Vigna et al. 2018). To reach the completeness level of 10 −3 for our system, we can neglect the terms which are of smaller order of magnitude in equation (7). This allows us to consider the VAs with a χvalue less than 5. If χ = 5 then exp(−χ 2 /2) 10 −5.4 , and the corresponding term is negligible. Moreover, the choice χ < 5 is valid for the score computation, because we are interested in a score accuracy of about 10 −2 .

Results
We are setting up a service dedicated to the scan of the Minor Planet Center NEO Confirmation Page 4 (NEOCP). The goal is to identify asteroids as NEOs, MBOs or distant objects to be confirmed or removed from the NEOCP, and to give early warning of imminent impactors, to trigger immediately follow-up observations. The software used to produce these results is a new version of the OrbFit Software version 5.0 5 .
The service involves the following steps, based on the algorithm presented in Sections 2.2 and 3.
• Scanning of the NEOCP every 2 minutes. New cases or old cases just updated are immediately run.
• Computation and sampling of the AR using a 2dimensional representation in the (ρ, 9 ρ) plane with a either grid or a spider web.
• Computation of the MOV, obtaining a set of VAs.
• Propagation of the VAs in the future (currently for 30 days).
• Projection on the Modified Target Plane, searching for VIs.
• If VIs exist, computation of the Impact Probability.
• Computation of the score.
The time required to run one target strictly depend on the characteristics of the object, but usually it is between 15 and 20 minutes. When predicting possible imminent impacts, one of the most important requirements to fulfil is to minimize the number of unjustified alarms. We mark as non-significant cases the objects for which there are less than 3 observations or the arc length is less than 30 minutes (see also Sec. 5), unless there exists a nominal solution with a geodesic curvature SNR greater than 1. The classification of a case as non-significant does not mean we skip the computation: we anyway perform all the steps of the algorithm, and assign the score and the impact flag. Nevertheless, being non-significant automatically decreases the priority of the object in case of an alarm.
Unfortunately, these techniques are not enough to remove all the spurious cases. They usually occur when the astrometry is either known to be erroneous or noisy, or anyway not reliable. We cannot solve this problem, and we acknowledge that the astrometric error models based on large number statistic are not enough to distinguish erroneous and accurate astrometry in a small sample (see comments in Section 7).
We test our algorithm on the two well known cases of NEAs that have impacted the Earth few hours after the discovery, namely 2008 TC 3 and 2014 AA. We have already pointed out that the choice of the weights is very important in these cases: to be able to compare the results with Farnocchia et al. (2015b), we choose the same weights. Furthermore, we also select some cases among the objects that won't impact to also show the importance of the score computation.

Graphical representation of the results
We use plots showing the AR and its sampling to present our results. Hereinafter we describe the color code present in our figures. Concerning the AR, we make use of the following lines.
• The red solid line represents the level curve of the heliocentric energy equal to −k 2 /(2a max ). Namely, it is the outer boundary of the AR, corresponding to the boundary of the region defined by condition 1 in Section 2.1.
• The green dashed line shows where the geocentric energy is equal to 0, also taking into account the condition about the radius of the Earth sphere of influence, as discussed in Section 2.1 • The magenta dashed line (which is parallel to the rangerate axis) represents the shooting star limit condition.
• The magenta solid lines (which are parallel to the rangerate axis) represent different values of the absolute magnitude.
We now provide a description of the colors used for the sampling points. No point is marked if the 4-dimensional differential corrections does not converge, because the point does not belong to the MOV.
• The dots are black if χ > 5.
• In case a VI has been found, we mark with red circles the points representing possible impacting orbits.
• The orange star represents the point with the minimum χ 2 value. 3 2008 TC 3 has been discovered by Richard A. Kowalski at the Catalina Sky Survey on October 7, 2008. The object was spotted 19 hours before the impact, and it is the first body to be observed and tracked prior to falling on the Earth. After the discovery, hundreds of astrometric observations were submitted to the Minor Planet Center (MPC) and these observations allowed the computation of the orbit and the prediction of the impact. We use the first tracklet composed by 4 observations, and then the first two tracklets (7 observations) to ascertain whether we could predict the impact. We compute a uniform densified grid in log 10 (ρ) (Fig. 2, top panel) where we consider only the first 4 observations, while we are able to compute a reliable nominal orbit, and the consequent spider web using 7 observations (Fig. 2, down  panel). Table 3 shows that with 4 observations and using the grid we are able to predict a possible impact of the object with the Earth with an impact probability of 3.6%, and the score of the object to be classified as a NEA is 100%. This would have produced an alert for the observers that could have immediately followed-up the object. With 7 observations we can confirm the certainty that the asteroid is a NEA (score = 100%), and the impact probability grows to 99.7%.

Asteroid 2014 AA
2014 AA has been discovered by Richard A. Kowalski at the Catalina Sky Survey on the new year's eve of 2014. The object was discovered 21 hours before the impact, but it has not been followed-up as 2008 TC 3 because of the exceptional night in which it has been spotted. We use first the first tracklet composed by 3 observations, and then the whole set of 7 observations to test whether we could have predicted the impact with our method.
These two examples have several analogies: we compute a uniform densified grid in log 10 (ρ) with the first tracklet, that is only 3 observations (Fig. 3, top panel), and we are able to  compute a reliable orbit and the consequent spider web only with 7 observations (Fig. 3, down panel). Table 3 shows that using the first tracklet only, we are able to predict a possible impact with the Earth with an impact probability of 3.0%, and the NEA score of the object is 100%. Again, this would have produced an alert for the observers that could have im- mediately followed-up the object. As in the previous case, with the second tracklet we confirm both that the asteroid is a NEA (score = 100%) and the collision, since the impact probability grows up reaching the value of 100%.

Asteroid 2014 QF 433
The previous examples show how the systematic ranging is capable of identifying imminent impactors. Although this is one of the most important applications of this technique, the systematic ranging is also essential in the first short arc orbit determination process.
Asteroid 2014 QF 433 has been discovered by F51 -Pan-STARRS 1, Haleakala on 2014 August 26. The first four observations have been posted on the NEO Confirmation Page, with the temporary designation TVPS7NV. It has been on the NEOCP until 2014 September 5. On that day (with 18 observations) it has been confirmed to be a distant object by the Minor Planet Center. Figure 4 shows the results of the systematic ranging on this asteroid, with the four discovery observations only, and 51 minutes of arc length. In this case the AR has two connected components, indicating the possibility for the object to be distant. The values of the three positive roots of equation of degree 6 are r 1 = 1.103 au, r 2 = 40.072 au, and r 3 = 59.786 au. The attributable is A = (α, δ, 9 α, 9 δ) = (5.7358902, −0.3008327, −3.35275 · 10 −4 , −9.94065 · 10 −5 ), with α and δ in radians and 9 α and 9 δ in radians per day. The two plots in Figure 4 clearly show that the object is distant, since almost all the grid points corresponding to the MOV lie in the second connected component. As a consequence, the cumulative score for the Distant and Scattered classes is 99%.
As a further validation, we take the orbital elements of this asteroid from the AstDyS database 6 , and we compute the range and the range-rate at the epoch of the attributable. The result is shown in the down panel of Figure 4: the black star represents the orbit of 2014 QF 433 computed with all the available observations, and it is in perfect agreement with the systematic ranging sampling. 6 Asteroid Dynamics Site, available at http://hamilton.dm.unipi. it/astdys/

Asteroid 2017 AE 21
The case of 2017 AE 21 shows the importance of the score computation. An object is worthy of attention and has to be followed-up even though it is not an impactor: for instance, it could be a potential NEA.
Asteroid 2017 AE 21 has been discovered by F51 -Pan-STARRS 1, Haleakala on 2017, January 3. It appeared on the NEOCP as a tracklet of 3 observations spanning 30 minutes, with the temporary designation P10yBuc. It has been confirmed to be a NEA on 2017 January 24, when it had 5 observations. With the first tracklet, our system produces an impact flag of 2, indicating a modest impact risk, with an impact probability IP = 2 · 10 −3 . Moreover, the NEO score is 92%, encouraging some follow-up to confirm. The top panel of Figure 5 shows the result when using the first tracklet only. We do not have any reliable nominal orbit to use, and as a consequence we adopt the grid sampling. The portion of the grid corresponding to low χ values (blue points) is quite wide, indicating a great uncertainty in the orbit determination, and the uncertainty region contains also impacting solutions.
With just two additional observations, the differential corrections still fail in computing a reliable nominal orbit, but now the good portion of the grid is located in a small subregion of the AR (see Figure 5, down panel). In this case the uncertainty region does not contain impacting orbits, thus we get an impact flag of 0, with IP = 0, whereas the NEO score grows to 100%. As a consequence, the low probability VI has been contradicted by the new observation, but the follow-up suggestion coming from the high 98% NEO score of the first run was reliable.

NEOCP object P10vxCt
As we stated in the introduction of this section, noisy astrometry can be the cause for unjustified alarms. In fact, if an object has a single tracklet of few observations and one of them is erroneous, the arc usually shows a significant curvature, implying that the object seems very close and fast moving. Most likely, it could be classified as an immediate impactor with very high impact probability.
Object P10vxCt has been spotted by F51 -Pan-STARRS, Haleakala on 2016 June 8. The first time it appeared on the NEOCP it had a tracklet with 3 observations spanning about 44 minutes (Table 4, above). It has never been confirmed, but it is anyway an important example to show the risk posed by noisy astrometric data. With the first 3 observations, our system computes a nominal solution compatible with a very close orbit, resulting in a spider web sampling over a small subset of the AR (see Fig. 6, top panel). A very large fraction of the MOV orbits are impacting solutions, and it results in an impact probability of 99.2% and impact flag 4, considering the significance of the curvature. The second batch of observations consists of 4 positions: 3 of them are a remeasurement of the first tracklet obtained from the discovery images of the object, plus an additional observation. With this new astrometry, the impact has been ruled out and the object has been removed from the NEOCP.  To show the role of the remeasurements in the impact removal, we consider the 3 remeasured observations only (see Table 4, below). The second observation in the first tracklet was badly determined, being off by about 3 arcsec from the corresponding one in the second batch. The effect of this shift can be seen in the curvature parameters κ (geodesic curvature) and 9 η (acceleration). For the first tracklet we have κ 1 = (0.0010073 ± 0.0001015) and 9 η 1 = (0.0003218 ± 0.0001013), whereas for the remeasured tracklet κ 2 = (0.0000649 ± 0.0006749) and 9 η 2 = (−0.0000430 ± 0.0006750), both significantly lower than the first ones (see Fig. 7 for a graphic representation of the two arcs). Moreover, as we can see from the curvature uncertainties, both curvature components are not significantly different from 0. As a consequence, with the remeasured observations only, the impact solution is sharply downgraded: a nominal solution cannot be computed anymore, resulting in a grid sampling of the AR, and the impact orbits are a very small fraction of the MOV orbits (see Fig. 6, right panel). Thus the impact probability lowers to about IP = 7.5 · 10 −5 , with an impact flag of 1. Providing remeasured observations is not the only way to solve the problem caused by bad astrometry. The second observation is not as good as the other two, and let us suppose this information were provided along with the observation itself. In this case, we could have properly down-weighted the second observation to take into account the additional information, and the case would have been solved. To prove this claim, we assign a formal uncertainty of 3 arcsec to both the right ascension and the declination of the second observation. With this choice, the impact solution still remain, but with an impact probability IP = 4.4 · 10 −4 . Until this additional metadata will be provided together with the observations, cases  Table 4, for NEOCP object P10vxCt. The red line represents the originally submitted tracklet, whereas the blue one the remeasured tracklet. The higher curvature of the first arc with respect to the second is clear.
like the one presented here can be solved only by a manual intervention after all the computations (remeasurement) or by a fast follow-up (see Section 7 for general comments on this issue).

The ESA Gaia mission and short arc orbit determination
The ESA Gaia mission, currently surveying the sky from the Sun-Earth L 2 Lagrangian Point, is providing astrometry of stars and asteroids, at the sub-milliarcsec accuracy (Prusti 2012) down to magnitude V = 20.7. The spin of the satellite is 6 h, and it operates in a continuous scanning mode. It has two lines of sight, separated by an angle of 106.5 • in its scanning direction. The continuous mode of observation implies targets are not pointed at, but are rather passing in front of Gaia fields of view. Such crossings are called transits. Over 5 years of nominal mission duration, the objects observed by Gaia will have a coverage of 80 − 100 observations for an average direction (60 − 70 for the ecliptic (Tanga et al. 2016)) Gaia focal plane is a large Giga-pixel array of 106 CCDs. The CCDs are organized as follows.
• The first two CCD strips are devoted to the source detection. This is the instrument called Sky Mapper (SM).
• Other CCD strips are devoted to low resolution spectrophotometry (red and blue photometry, RP/BP) and high resolution spectroscopy (RVS), which is not considered for asteroid studies.
Each Solar System Object (SSO) transit is composed, at most, by 10 astrometric observations (AF and SM instruments), distributed over 50 seconds. The Data Processing and Analysis Consortium (DPAC) has, as part of its activities, the source identification and further processing done on the ground. In this context, the Coordination Unit 4 (CU4) performs the analysis of objects deserving a specific treatment, as Solar System object (PI: P. Tanga). The DPAC CU4 has implemented two pipelines for Solar System processing (Tanga et al. 2007;Mignard et al. 2007): • SSO-ST is the Solar System short-term processing, devoted to provide a first, approximate orbit for the recovery of objects potentially discovered by Gaia. A groundbased follow-up network (Gaia-FUN-SSO, Thuillot et al. (2014)) is currently operating, realizing follow-up observations of Gaia potential discoveries from the ground; • SSO-LT is the Solar System long-term processing, which runs for the data releases, performing a more sophisticated data reduction with the best possible instrument calibration and astrometric solution.

The short-term processing
The Solar System short-term processing is based on few transits of the object (at least 3), covering a time span of few hours (at least six). At now, the astrometry for the alerts is based on a preliminary calibration, which is not the same used for the long-term processing. As a result, the error model required by these observations is not different from the one already used for the ground-based cases with the best astrometry. We uniformly weight (0.1 arcsec) both right ascension and declination.
If the detected source is not successfully linked with a known solar system object, then it is potentially a discovery. It is thus crucial to predict a possible sample of orbits for the ground-based follow-up network to certify the discovery. In the framework of the CU4 data treatment, this is done by random-walk statistical ranging, which has been developed by Muinonen et al. (2016). This is the Java code which is currently producing the orbital data for Gaia alerts in the SSO-ST pipeline.
We take the opportunity of the availability of the method presented in Sec. 2 to validate independently the results of the SSO-ST pipeline.
The impact probability computation, which is a key feature when we look for possible impactors, is not so essential when we have to deal with Gaia short-term observations. We expect indeed that among the objects that Gaia will discover, there will be a large fraction of Main Belt Asteroids, few Near Earth Asteroids, and it would be very surprising if it could discover even an imminent impactors (Carry 2014). On the other side, the use of the double grid or of the cobweb is essential in this case.

Results
The same service presented in Section 4 to scan the NEO Confirmation Page, has been ran on possible Gaia discoveries. The graphical representation of the results is identical to the one given in Section 4.1. We also use the same version of the OrbFit software cited in Sec. 4, adapted to deal with a different input given by Gaia observations. All the alerts are available at https://gaiafunsso. imcce.fr/, and hereinafter we present two particular objects, that at the time seemed to be possible Gaia discoveries, but then their observations have been linked to ground-based observations submitted earlier in time.
The first object that could have been discovered by Gaia on 2016, December 29, is a Main Belt Asteroid. It has been identified as g0T015. It has four Gaia transits which cover a time span of 16 hours.
We apply the systematic ranging (using the two grids) on the observations, and the results are summarized in Fig. 8. The object has been then followed-up by the Observatoire de Haute-Provence (OHP) for two consecutive nights (2017, January 3-4), and the observations have been reported to the Minor Planet Center. The object has been recognized as a Main Belt Object, and it obtained the provisional designation 2017 AD 17 . It could have been the first Gaia discovery, if it hasn't been linked to few observations from F51 Pan-STARSS submitted earlier (March 2014).
To be sure that the observed object was the same potentially discovered by Gaia, we consider as initial guess the elements corresponding to the point with the minimum value of χ 2 , and we perform an orbit determination process including the outlier rejection procedure (Milani & Gronchi 2010).  Figure 9 shows that follow-up observations from the OHP, and Gaia observations match together, according to the orbit selected as first guess. The residuals for the OHP (blue points) are way larger than the ones obtained from Gaia observations (red points), as expected. The weights used in this case are uniform in right ascension (RA) and declination (DEC), and correspond to 0.1 arcsec. The result also shows a high correlation (close to 1) between RA and DEC, which is typical for Gaia observations due to its scanning law.
The second Gaia object is g1 j0D7, again a MBA. It has been observed by Gaia on 2017, September 2. It has 7 transits, which cover a time span of ∼ 22 hours. The object has followup observations from the Abastuman Observatory (MPC code 119) the night between September 10 and September 11. It obtained the MPC preliminary designation 2017 RW 16 , but then it has been linked to the asteroid 2006 UL 1 89 discovered by the Catalina Sky Survey on June 2005. Figure 10 shows the result of the systematic ranging when we use only Gaia observations. Again, we choose the point with the minimum value of the χ 2 as starting point, and we use it as preliminary point. We then perform a differential corrections least-squares fit, and we obtain the residuals (see Fig. 11), as described in the previous case.  Figure 11: Residuals in right ascension and declination, in arcsec, obtained as a result of the orbit determination applied to the whole set of observations: Gaia (red points) and ground based (blue points) for the Gaia object g1 j0D7.
Then, when the accuracy for the short term will improve, we will also be able to run the systematic ranging on less than three transits, given the correct error model to the observations. Here it is worth noting that the concept of short arc strongly depends on the concept of curvature (see Sections 2.2 and 4), which is related not only to the time interval, but also to the accuracy of the observations themselves.
The score computation represents a key feature in the Gaia frame (as already pointed out in Sec. 5.1), but it will also be very useful in future applications, like the ESA Euclid mission. Euclid is an ESA mission with the aim of mapping the geometry of the dark universe down to V AB =24.5, with a launch scheduled for 2020. While conducting its primary goal survey, Euclid will also observe asteroids during his whole lifetime, and a Solar System Working Group has been created within the Euclid consortium.
Euclid will observe at solar elongation ∼ 91 • , and each SSO will be imaged 16 times (more precisely, 95% of them 12 times, and 50−60% of them 16 times) over 67 minutes (Carry 2017). With a Hubble-like angular resolution, astrometry will be quite good, at least at the same level of the best observatories on ground (like Mauna Kea), or as the short term accuracy of the Gaia alerts as it is now. The estimate accuracy is around 100 mas. While the southern sky will be repeatedly covered to the same depth by LSST (LSST Science Collaborations & LSST Project 2009), only Euclid will systematically cover high declinations, with a strong potential for discoveries. For each, having the score will be crucial to select the asteroid for which trigger or not follow-up observations.

Conclusions
One of the main issue in the impact hazard assessment for imminent impactors is given by the computation of the impact probability. The main results of this article are a new algorithm to propagate the probability density function from the space of the astrometric residuals to the Manifold Of Variation, a geometric device to sample the set of possible orbits, available even after a very short observed arc. In previous works, this computation was supported with an a priori number density of asteroids. Our computation is complete, rigorous, and uses no a priori hypothesis.
Does this new algorithm solve the problem of assessing the risk of imminent impacts from a freshly discovered asteroid, with observations limited to 1-2 tracklets? By using the AR and one of our grid sampling, we have shown how to approximate a probability integral on the portion of the MOV leading to an imminent impact, if it is found. However, to accept this integral as Impact Probability we need to check three conditions.
First, the probability density on the space of residuals needs to be based upon a probabilistic model of the astrometric errors, taking into account the past performances of the observatories. Second, the observations used in the computation must be "typical" of the observatory: even the best astronomical program produces a comparatively small subset of "faulty" observations, with errors much larger than the usual ones. Third, we should assume that the small sample of observations has statistical properties, such as mean and standard deviation (STD), close to the ones of the full distribution.
The first hypothesis is reasonable, in that a lot of work has been developed in the last 20 years to produce astromet-ric error models for asteroid observations (see Carpino et al. (2003), Chesley et al. (2010), Baer et al. (2011), Farnocchia et al. (2015a). These models are not perfect, but they represent an increasingly reliable source of statistical information. The second hypothesis is not trivial: the current format for asteroid observations does not contain sufficient metadata to discriminate the "weak observations" from the good ones. The full adoption of the new Astrometric Data Exchange Standard (ADES 7 ), approved by the IAU in 2015, will provide information such as SNR, timing uncertainty, and so on, allowing to adapt the weighting of the individual observations. The example of P10vxCt shows how just one lower quality observation can completely spoil the orbit results, generating a false impact alarm. This can be avoided either with remeasuring by the observer or by the orbit computer, provided such down-weighting is supported by the metadata.
The third hypothesis is the most troublesome. Assuming that the probability density of an astrometric error model is a perfect statistical description, then by the law of large numbers a large enough sample of N observations shall have approximately the same statistical properties of the model, with the differences going to zero for N → +∞ (law of large numbers). Unfortunately, N = 3, 4, 5 is not large enough for the law of large numbers to apply. For instance, a tracklet with N = 3 observations can have all the observations in one coordinate with errors > 2.5 STD: this statistical fluke would be very rare, occurring in a little more than 1 tracklet over 1 million. Still, if a large asteroid survey submits to the MPC more than one million tracklets per year, such a fluke may occur about once a year, whereas the discovery of imminent impactors is currently more rare (2 in 10 years). Detection of a rare astronomical event cannot be a priori discriminated from rare statistical events.
The tests on real cases discussed in this paper, and many more from the NEOCP, convinced us that our algorithm computes a reliable impact probability when the impact actually occurs. Nevertheless, we cannot show that our algorithm is immune from "false" alarms. They are not false in the sense of a wrong computation, or even worse a malicious disinformation, they are statistical flukes which cannot be avoided because of lack of information (hypothesis 2) and the need to use statistics on a small sample (hypothesis 3). The question is what should be done to mitigate the damage by these false alarms, given that we cannot avoid disseminating them: otherwise, how could we disseminate the alarm in the true case?
The only answer is to have a follow-up chain which does not waste resources: the discoverers could themselves either remeasure or follow-up on the short term, like 1 hour after discovery, the cases announced as possible impactors. Other telescopes should be available to perform follow-up, to avoid improper use of survey telescopes for a less demanding task. The ideal solution should be the availability of a Wide Survey, capable of covering the entire dark sky every night and of detecting, e.g., an asteroid with absolute magnitude H = 28 at 0.03 au distance (near opposition). Then the same asteroid 7 It is available at http://minorplanetcenter.net/iau/info/ IAU2015_ADES.pdf. would be recovered by the survey the next day, before the impact, and without the need for auxiliary follow-up. acknowledgements This work has made use of data from the European Space Agency (ESA) mission Gaia (http://www.cosmos.esa. int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, http://www.cosmos.esa.int/ web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.
Part of this work is based on observations made at Observatoire de Haute Provence (CNRS), France. The authors would also like to acknowledge all the observers involved in the detection of the asteroids discovered by Gaia: for g0T015 they are E. Saquet, M. Dennefeld, D. Gravallon, and for g1j0D7 they are R. Y. Inasaridze, V. Ayvazian, T. Mdzinarishvili, Y. Krugly.
F. Spoto acknowledges support by the CNES post-doctoral program, and A. Del Vigna acknowledges support by the company SpaceDyS.

A Probability density function computation
In this appendix we give the mathematical details for the derivation of equation (6) for the probability density function on the sampling space S .
namely the normal matrix of the constrained differential corrections leading to x 0 , computed at convergence. Hence having used (2), that is mQ(x) = mQ * + χ 2 (x). Concerning the choice of x 0 we proceed as follows: if a reliable nominal solution x * exists, we set x 0 = x * ; if not, we select as x 0 the sample orbit having the minimum value of the target function. This choice is also coherent with the χ computation given by (2).

A.2 From the MOV to the Admissible Region
We can now compute the determinant of the map f µ . Let M µ be the 2 × 2 matrix representing the tangent map between K and M. f µ is a differentiable function, with Jacobian matrix We now consider ρ 0 ∈ K such that f µ (ρ 0 ) = x 0 . The matrix B(ρ 0 ) has rank 2, thus M is smooth in the neighborhood of each of its points. We can linearize in ρ 0 , obtaining the tangent map B(ρ 0 ) = D f µ (ρ 0 ) : which is a linear map between two 2-dimensional spaces. We use a rotation matrix R in the orbital elements space, such that