© ESO 2018

1 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) of finding 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) defined 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 were 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: statistical and systematic methods.

The original statistical ranging method (Virtanen et al. 2001; Muinonen et al. 2001) starts from the selection of a pair of astrometric observations. Then, 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 χ² value between the observed and computed observations is within a pre-defined threshold. Oszkiewicz et al. (2009) improved the statistical ranging using the Markov-chain Monte Carlo (MCMC) to sample the phase space. The MCMC orbital ranging method is based on a bivariate Gaussian proposal PDF for the topocentric ranges. Then, Muinonen et al. (2016) have developed a random-walk ranging method in which the orbital-element space is uniformly sampled, up to a χ² 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. These authors have developed this method for the European Space Agency (ESA) Gaia mission, in the framework of Gaia alerts on potentially new discovered objects by Gaia (see Tanga et al. 2016).

On the other hand, Chesley (2005) and Farnocchia et al. (2015b) introduced the so-called systematic ranging, which systematically explores a raster in the topocentric range and range-rate space . This technique enables the description of 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 that lead to an impact with the Earth.

In this paper we describe a new approach to the systematic ranging, based on knowledge of the admissible region (AR; Milani et al. 2004), and a new method to scan the region. The process has some main advantages to other methods described above as follows:

• 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 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 is based on a refinement using the value of the post-fit χ² of each point in the first grid (see Sect. 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 Sect. 3).

In Sect. 2.1 we summarize the features of the AR, and describe the sampling of the space needed for the creation of the manifold of variations (MOV), which is an extension of the line of variations (LOV; Milani et al. 2005b; Milani & Gronchi 2010). In Sect. 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 Sect. 4 we apply our method to the well-known cases 2008 TC₃ and 2014 AA, and we show other examples to explore the capabilities of the new method. Then in Sect. 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

2.1 Admissible region and systematic ranging

Even though the observations are too scarce, we are able to compute the right ascension α, the declination δ, and their time derivatives and , by fitting both angular coordinates as a function of time with a polynomial model. These four quantities could be assembled together to form the attributable (Milani & Knežević 2005): (1)

at a chosen time , which could be the time of the first observation or the mean of the observation times. The information contained in the attributable leaves the topocentric distance ρ and radial velocity completely unknown. We would have a full description of the topocentric position and velocity of the asteroid in the attributable elements , if ρ and were known. From then on, we use , where , to describe the attributable elements.

Given an attributable, we define the AR as the set of all the possible couples 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 only consider the objects for which the value of the heliocentric energy is less than − k²∕(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 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 space. The AR usually has 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 there are three roots and 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_SI), or by a segment of the straight vertical line ρ = R_SI and two arcs of the zero curve of the geocentric energy (for 0 < ρ < R_SI). We also discard the orbits corresponding to meteors that are too small to be sources 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 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₁, r₂ and r₃, in ascending order). Table 1 summarizes the conditions and grids used in the 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¹.

The target function is defined by

where are the fit parameters, m is the number of observations used in the least squaresfit, ξ is the vector of the observed-computed debiased astrometric residuals², 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, for example the Gaia DR1 (Lindegren et al. 2016).

Given a subset K of the AR, we define the MOV as the set of the points such that ρ₀ ∈ K and is the local minimum of the function . In addition, the value of the minimum RMS of the residuals is less than a given threshold Σ. In general, the MOV is a two-dimensional manifold, such that the differential of the map from the sampling space to 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 we fix ρ = ρ₀ and in the target function, and then we search by means of an iterative procedure, the doubly constrained differential corrections. The normal equation is

where

We indicateas K′ the subset of K on which the doubly constrained differential corrections converges. In this way, the sampling of the MOV is performed over K′.

For eachpoint x on the MOV, we also compute a χ value (2)

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 among all the values of the points for which we have a convergence of the four-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 that an object belongs to various 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

where a is the semimajor axis (in au) and e is the eccentricity.

• DO: distant object, characterized by q > 28 au (e.g., a Kuiper belt object, hereafter KBO).

• SO: scattered object, not belonging to any of the previous classes.

If the object is close to being a NEO (the NEO-score is more than 50%), we use a uniform grid in log ₁₀(ρ), otherwise we use a uniform grid in ρ. Then we compute again the MOV on thenew and denser 100 × 100 grid.

Fig. 1 Example of spider web around the nominal solution . 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.

2.2 Spider web

Let us suppose that a nominal orbit was 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 (S/N) (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 a reliable nominal solution already.

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 = Γ⁻¹. The matrix C is positive definite, hence the level curves of the target function are concentric five-dimensional ellipsoids in the six-dimensional orbital elements space. The level curves on the space are represented by the marginal ellipsoids, defined by the normal matrix

where Γ_ρρ is the restriction of Γ to the 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³. Then we apply the following transformation to each point of the grid, depending on the covariance matrix of the nominal orbit and on the orbit itself: (3)

where λ₁ > λ₂ are the eigenvalues of the 2 × 2 matrix Γ_ρρ, v₁ is the eigenvector corresponding to the greatest eigenvalue λ₁, and ρ^* and are the range and range-rate values of the nominal solution. Figure 1 shows an example of spider web sampling on the plane.

3 Probability density computation

We obtain aprobability 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⁻¹. Hence the probability density function on the residuals space is (4)

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 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. This prior was 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 . Jeffreys’ prior tends to favor orbits in which the object is close to the observer, because of the sensitivity of the residuals for small topocentric distances. Granvik et al. (2009) made versions of the code, which uses 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 ρ² by geometric considerations on the Cartesian space of position and velocity.

• Uniform distribution. Uniform distribution in the 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 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 MOV), 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: if the sampling is uniform in ρ, if the sampling is uniform in log₁₀ρ, and 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.

• is the MOV, a two-dimensional manifold in the six-dimensional orbital elements space .

• is the residuals space. The residuals are a function of the fit parameters: ξ = F(x), with differentiable, and we define the manifold of possible residuals as (Milani & Gronchi 2010, Sect. 5.7).

Without loss of generality, we can assume that (5)

where I_m is the m × m identity matrix. As explained in Milani & Gronchi (2010, Sect. 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 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 χ²(s) to indicate χ²(x(ρ(s))). The probability density function of S is (6)

where M_μ is the 2 × 2 matrix associated with f_μ (considered as tangent map to the surface ), 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 Sect. 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.

3.1 Impact probability computation

Each point on the MOV can be thought as orbit compatible with the observations, and we call each point a 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 with a subset of the sampling space, and hence its probability is given by (7)

If for a given object we find impacting solutions, we assign an impact flag, which is an integer number related to the computation of the impact probability, to the object. The impact flag depends on the impact probability and on the arc curvature, as shown in Table 2. An arc has significant curvature if χ² > 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%). The impact flag is conceived as a simple and direct communication tool to assess the importance of collision predictions and give priority to the follow-up activities.

The goal for a system dedicated to imminent impactors is to detect all the possible VIs down to a probability level of about 10⁻³, called completeness level (Del Vigna et al. in prep.,). To reach the completeness level of 10⁻³ for our system, we neglect the terms that are of smaller orders of magnitude in Eq. (7). This allows us to consider the VAs with a χ-value less than 5. If χ = 5 then exp (−χ²∕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⁻².

4 Results

We are setting up a service dedicated to scanning the Minor Planet Center NEO Confirmation Page⁴ (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 follow-up observations immediately. The software used to produce these results is a new version of the OrbFit Software version 5.0⁵.

The service involves the following steps, based on the algorithm presented in Sects. 2.2 and 3.

• Scanning of the NEOCP every two minutes. New cases or old cases just updated are immediately run.

• Computation and sampling of the AR using a two-dimensional representation in the 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 MTP, searching for VIs.

• If VIs exist, computation of the IP.

• Computation of the score.

The time required to run one target strictly depends 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 fulfill is to minimize the number of unjustified alarms. We denote as nonsignificant cases the objects for which there are less than three observations or the arc length is less than 30 minutes (see also Sect. 5), unless there exists a nominal solution with a geodesic curvature S/N greater than 1. The classification of a case as nonsignificant does not mean we skip the computation. We perform all the steps of the algorithm in any case and assign the score and impact flag. Nevertheless, being nonsignificant 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 sufficient to distinguish erroneous and accurate astrometry in a small sample (see comments in Sect. 7).

We test our algorithm on the two well-known cases of NEAs that have impacted the Earth a few hours after the discovery, namely 2008 TC₃ and 2014 AA. We already pointed out that the choice of the weights is very important in these cases. We choose the same weights in order to be able to compare the results with Farnocchia et al. (2015b),. Furthermore, we also select some cases among the objects that would not be impacted to show the importance of the score computation as well.

4.1 Graphical representation of the results

We present our results with plots showing the AR and its sampling. 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²∕(2a_max). Namely, it is the outer boundary of the AR, corresponding to the boundary of the region defined by condition (1) in Sect. 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 Sect. 2.1

• The magenta dashed line (which is parallel to the range-rate axis) represents the shooting star limit condition.

• The magenta solid lines (which are parallel to the range-rate axis) represent different values of the absolute magnitude.

We now provide a description of the colors used for the sampling points. No point is denoted if the four-dimensional differential corrections does not converge because the point does not belong to the MOV.

• The dots are indicated in blue if χ ≤ 2, and green if 2 < χ ≤ 5.

• The dots are indicated in black if χ > 5.

• In case a VI has been found, we show the points representing possible impacting orbits with red circles.

• The orange star represents the point with the minimum χ² value.

4.2 Asteroid 2008 TC₃

2008 TC₃ was 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 four observations, and then the first two tracklets (seven observations) to ascertain whether we could predict the impact.

We compute a uniform densified grid in log₁₀(ρ) (Fig. 2, top panel) where we consider only the first four observations, while we are able to compute a reliable nominal orbit, and the consequent spider web using seven observations (Fig. 2, bottom panel). Table 3 shows that with four 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 seven observations we can confirm the certainty that the asteroid is a NEA (score = 100%), and the impact probability increases to 99.7%.

Fig. 2 Grid sample of the space for the first 4 observations(top panel) and spider web for the first 7 observations(bottom panel) of 2008 TC3.

Fig. 3 Grid sample of the space for the first 3 observations(top panel) and spider web for the whole set of 7 observations(bottom panel) of 2014 AA.

4.3 Asteroid 2014 AA

2014 AA was 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₃ because of the exceptional night in which it has been spotted. We initially used the first tracklet composed of three observations, and then the whole set of seven 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 ₁₀(ρ) with the first tracklet, which is only 3 observations(Fig. 3, top panel), and we are able to compute a reliable orbit and the consequent spider web only with seven observations (Fig. 3, bottom 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 immediately 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 increases to 100%.

Fig. 4 Grid sample of the space for the first 4 observations of 2014 QF433 (top panel) and an enlargement of the second component (bottom panel). The black star in the second component represents the orbit of the object computed with all the available observations.

4.4 Asteroid 2014 QF₄₃₃

The previous examples show how systematic ranging is capable of identifying imminent impactors. Although this is one of the most important applications of this technique, systematic ranging is also essential in the first short arc orbit determination process.

Asteroid 2014 QF₄₃₃ was discovered by F51 - Pan-STARRS 1, Haleakala on August 26, 2014. The first four observations were posted on the NEO confirmation page, with the temporary designation TVPS7NV. This asteroid has been on the NEOCP since September 5, 2014. On that day (with 18 observations) it was confirmed to be a distant object by the Minor Planet Center.

Figure 4 shows the results of systematic ranging on this asteroid with only four discovery observations and 51 minutes of arc length. In this case the AR has two connected components, indicating the possibility that object is distant. The values of the three positive roots of equation of degree 6 are r₁ = 1.103 au, r₂ = 40.072 au, and r₃ = 59.786 au. The attributable is (8)

with α and δ in radians and and in radians per day. The two plots in Fig. 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⁶, and we compute the range and the range-rate at the epoch of the attributable. The result is shown in the bottom panel of Fig. 4: the black star represents the orbit of 2014 QF₄₃₃ computed with all the available observations and is in perfect agreement with the systematic ranging sampling.

4.5 Asteroid 2017 AE₂₁

The case of 2017 AE₂₁ shows the importance of score computation. This object is worthy of attention needs follow up even though it is not an impactor; for instance, it could be a potential NEA.

Asteroid 2017 AE₂₁ was discovered by F51 - Pan-STARRS 1, Haleakala on January 3, 2017. It appeared on the NEOCP as a tracklet of three observations spanning 30 minutes with the temporary designation P10yBuc. This object was confirmed to be a NEA on January 24, 2017, when it had five observations. With the first tracklet, our system produces an impact flag of 2, indicating a modest impact risk and an impact probability IP = 2 × 10⁻³. Moreover, the NEO score is 92%, which encouraged some follow up observations of the object. The top panel of Fig. 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 vey wide, indicating a great uncertainty in the orbit determination, and the uncertainty region also contains impacting solutions.

With just two additional observations, the differential corrections still fail to compute a reliable nominal orbit, but now the good portion of the grid is located in a small subregion of the AR (see Fig. 5, bottom panel). In this case the uncertainty region does not contain impacting orbits, thus we get an impact flag of 0, where IP = 0, whereas the NEO score increases to 100%. As a consequence, the new observations contradicts the low probability VI, but the follow-up suggestion coming from the high 98% NEO score of the first run is reliable.

4.6 NEOCP object P10vxCt

As we stated in the introduction of this section, noisy astrometry can be the cause of unjustified alarms. In fact, if an object has a single tracklet of few observations and one of these 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 was spotted by F51 - Pan-STARRS, Haleakala on June 8, 2016. The first time it appeared on the NEOCP it had a tracklet with three observations spanning about 44 minutes (Table 4, above). It has never been confirmed, but it is in any case an important example to show the risk posed by noisy astrometric data. With the first three 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 percentage of the MOV orbits are solutions with possible impacts 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 four positions, three of which 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 was ruled out and the object was removed from the NEOCP.

To show the role of the remeasurements in the impact removal, we only consider the three remeasured observations (see Table 4, below). The second observation in the first tracklet was badly determined because it was off by about 3 arcsec from the corresponding observation in the second batch. The effect of this shift can be seen in the curvature parameters κ (geodesic curvature) and (acceleration). For the first tracklet we have (9)

while for the remeasured tracklet (10)

which are both significantly lower than the values obtained for the first, not remeasured, tracklet (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, the impact solution is sharply downgraded with the remeasured observations alone; 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⁻⁵, 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⁻⁴. Until this additional metadata will be provided together with the observations, cases 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 Sect. 7 for general comments on this issue).

Fig. 5 Grid sample of the space for the first 3 observations(top panel) and for the whole set of 5 observations(bottom panel) of 2017 AE21. In both cases we do not have any reliable nominal orbit to use, and as a consequence we adopt the grid sampling.

Fig. 6 Grid sample of the space for the first 3 observations(top panel) and for their remeasurement (bottom panel) of P10vxCt.

Fig. 7 Two batches of observations reported in Table 4, for NEOCP object P10vxCt. The red line represents the originally submitted tracklet, whereas the blue line the remeasured tracklet. The higher curvature of the first arc with respect to the second is clear.

5 ESA Gaia mission and short arc orbit determination

The ESA Gaia mission, which is currently surveying the sky from the Sun-Earth L₂ Lagrangian point, is providing astrometry of stars and asteroids, at the submilliarcsec 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. The satellite has two lines of sight that are 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 five 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)

The Gaia focal plane is a large Giga-pixel array of 106 CCDs. The CCDs are organized as follows:

• The firsttwo CCD strips are devoted to source detection. This is the instrument called Sky Mapper (SM).

• The following nine strips are astrometric CCDs, i.e. the astrometric field (AF).

• Other CCD strips are devoted to low resolution spectro-photometry, called red and blue photometry (RP and 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, source identification and further processing performed on the ground. In this context, the Coordination Unit 4 (CU4) performs the analysis of objects deserving a specific treatment as SSOs (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, providing a first, approximate orbit for the recovery of objects potentially discovered by Gaia. A ground-based 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.

5.1 Short-term processing

The solar system short-term processing is based on few transits of the object (at least three), 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 calibration used for long-term processing. As a result, the error model required by these observations is not different from that 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 performed by random-walk statistical ranging, which has been developed by Muinonen et al. (2016). This Java code currently produces the orbital data for Gaia alerts in the SSO-ST pipeline.

We take the opportunity of the availability of the method presented in Sect. 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 percentage of Main Belt Asteroids, few Near Earth Asteroids, and it would be very surprising if it could even discover imminent impactors (Carry 2014). On the other side, the use of the double grid or of the cobweb is essential in this case.

5.2 Results

The same service presented in Sect. 4 to scan the NEO Confirmation Page, was ran on possible Gaia discoveries. The graphical representation of the results is identical to that given in Sect. 4.1. We also use the same version of the OrbFit software cited in Sect. 4, which we adapted to deal with the different input given by Gaia observations⁷.

We present two particular objects, which at the time seemed to be possible Gaia discoveries, but whose observations have been linked to ground-based observations submitted earlier in time.

The first object that could have been discovered by Gaia on December 29, 2016 is a Main Belt Asteroid. It has been identified as g0T015. It has four Gaia transits that cover a time span of ≃ 16 hours.

We apply systematic ranging (using the two grids) on the observations and the results are summarized in Fig. 8. The object has then been followed up by the Observatoire de Haute-Provence (OHP) for two consecutive nights (January 3-4, 2017), and the observations have been reported to the Minor Planet Center. The object has been recognized as a MBO and obtained the provisional designation 2017 AD₁₇. This could have been the first Gaia discovery, if it had not 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 an initial guess the elements corresponding to the point with the minimum value of χ², 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, according to the orbit selected as first guess. The residuals for the OHP (blue points) are much larger than those 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 because of its scanning law.

The second Gaia object is g1j0D7, again a MBA. It has been observed by Gaia on September 2, 2017. It has seven transits, which cover a time span of ~22 hours. The object has follow-up observations from the Abastuman Observatory (MPC code 119) during the night between September 10 and September 11. This object obtained the MPC preliminary designation 2017 RW₁₆, but was then linked to the asteroid 2006 UL₁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 χ² 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.

Fig. 8 Admissible region and grid sample of the space for the Gaia object g0T015. Top panel: Grid sample of the space for the Gaia object g0T015. Bottom panel: Zoom-in of the second grid is shown.

Fig. 9 Residuals in right ascension and declination. Values in arcsec were 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 g0T015 (2017 AD17).

6 Future perspectives

Gaia alerts runs daily and they need a large effort to collect follow-up observations from ground. With the method and the softwaredescribed we validated the already existing Java code written for the alert pipeline. Moreover, our approach can also be used to confirm the discoveries by computing an orbit using both Gaia and ground-based observations.

Then, when the accuracy for the short term improves, we will also be able to run the systematic ranging on less than three transits, given the correct error model to the observations. It is worth noting that the concept of short arc strongly depends on the concept of curvature (see Sects. 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 Sect. 5.1), but it will also be very useful in future applications, such as 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; the Euclid launch is scheduled for 2020. While conducting its primary goal survey, Euclid also observes 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 over 67 minutes (Carry 2018 more precisely, 95% of them 12 times, and 50−60% of them 16 times). With a Hubble-like angular resolution, astrometry will be very good, at least at the same level of the best observatories on ground (likeMauna Kea), or at the level of 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 and have a strong potential for discoveries. For each, having the score will be crucial to select the asteroids the will trigger follow-up observations or not.

Fig. 10 Admissible region and grid sample of the space for the Gaia object g1j0D7. Top panel: Grid sample of the space for the Gaia object g1j0D7. Bottom panel: Zoom-in of the second grid is shown.

Fig. 11 Residuals in right ascension and declination. Values in arcsec were 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 g1j0D7.

7 Conclusions

One of the main issues 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 MOV, which is a geometric device to sample the set of possible orbits that are 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, for which observations are limited to 1–2 tracklets? By using the AR and one of our grid sampling methods, 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, we need to check three conditions to accept this integral as IP.

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 the so-called faulty observations that have errors much larger than the usual observations. Third, we should assume that the small sample of observations has statistical properties, such as mean and standard deviation (STD), close to those 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 astrometric 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 weak observations from the good observations. The full adoption of the new Astrometric Data Exchange Standard (ADES⁸), approved by the IAU in 2015, will provide information such as S/N, timing uncertainty, and so on, allowing adaptation of the weighting of 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 observationsshall 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.5STD; 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 so-called false alarms. These are not false in the sense of an incorrect computation, or even worse a malicious disinformation; rather they are statistical flukes that 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 these false alarms; otherwise, how could we disseminate a true alarm when we see one?

The only answer is to have a follow-up chain that does not waste resources. The discoverers could either remeasure or follow up in the short term, such as one hour after discovery, when the cases are 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 theentire dark sky every night and of detecting, for example, an asteroid with absolute magnitude H = 28 at 0.03 au distance (near opposition). Then the same asteroid 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 postdoctoral program, and A. Del Vigna acknowledges support by the company SpaceDyS.

Appendix A Probability density function computation

In this appendix we give the mathematical details for the derivation of Eq. (6) for the probability density function on the sampling space S.

A.1 From the residuals space to the MOV

The first step of the procedure is the classical propagation of the probability density function from the normalized residuals space to theorbital elements space. From Eq. (5) we recall that we start from the following density:

We consider a point , the corresponding image ξ₀ = F(x₀) ∈ V and the tangent application . Later we discuss the choice of the point x₀. As proved in (Milani & Gronchi 2010, Sect. 5.7), we obtain a probability density function on the orbital elements space given by

where is the covariance matrix resulting from the transformation of Gaussian random variables under a linear transformation. Moreover, the normal matrix of the variable X is

namely the normal matrix of the constrained differential corrections leading to x ₀, computed at convergence. Hence

having used Eq. (2), that is mQ(x) = mQ^* + χ²(x). Concerning the choice of x₀ we proceed as follows: if a reliable nominal solution x^* exists, we set x₀ = x^*; if not, we select as x₀ the sample orbit having the minimum value of the target function. This choice is also coherent with the χ computation given by Eq. (2).

A.2 From the MOV to the AR

We can now compute the determinant of the map f_μ. Let M_μ be the 2 × 2 matrix representing the tangent map between K′ and . f_μ is a differentiable function, with Jacobian matrix

We now consider ρ₀ ∈ K′ such that f_μ(ρ₀) = x₀. The matrix B(ρ₀) has rank 2, thus is smooth in the neighborhood of each of its points. We can linearize in ρ₀, obtaining the tangent map

which is a linear map between two two-dimensional spaces. We use a rotation matrix R in the orbital elements space, such that

This means that x^″ parametrizes . In these coordinates the linearized map has a simple structure:

with A(ρ₀) an invertible 2 × 2 matrix. By using that R is orthogonal, the following relation holds:

and hence (A.1)

The next step is to explicitly compute the matrix . Hereinafter we neglect terms containing the second derivatives of the residuals multiplied by the residuals themselves. The value is the attributable, which minimizes the target function . That is, is a zero of the function

The function F is continuously differentiable, and we have

where we used the approximation assumed at the beginning. The matrix is invertible, otherwise the doubly constrained differential corrections would fail, and the minimum point could not be reached. By applying the implicit function theorem, there exists a neighborhood U of ρ₀, a neighborhood W of , a continuously differentiable function f : U → W such that, for all ρ ∈ U holds

and (A.2)

We already computed , so we proceed with the other derivative, i.e.,

By using the Eq. (A.2) we obtain (A.3)

A.3 From the AR to the sampling space

The last step is the computation of the determinant of the map f_σ, and this depends on S, for which we have different possibilities. We call M_σ the Jacobian matrix associated with f_σ.

If the sampling is uniform in ρ, then f_σ is the identity map, and therefore detM_σ = 1. If the sampling is uniform in log₁₀(ρ),

we have

and hence

If we are in the cobweb case, the f_σ is given by Eq. (3). Its Jacobian matrix is (A.4)

where . We have omitted the dependency of M_σ on ρ₀ not to have a heavy notation. After some manipulation, the determinant is

References

All Tables

All Figures

	Fig. 1 Example of spider web around the nominal solution . 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.
In the text

	Fig. 2 Grid sample of the space for the first 4 observations(top panel) and spider web for the first 7 observations(bottom panel) of 2008 TC3.
In the text

	Fig. 3 Grid sample of the space for the first 3 observations(top panel) and spider web for the whole set of 7 observations(bottom panel) of 2014 AA.
In the text

	Fig. 4 Grid sample of the space for the first 4 observations of 2014 QF433 (top panel) and an enlargement of the second component (bottom panel). The black star in the second component represents the orbit of the object computed with all the available observations.
In the text

	Fig. 5 Grid sample of the space for the first 3 observations(top panel) and for the whole set of 5 observations(bottom panel) of 2017 AE21. In both cases we do not have any reliable nominal orbit to use, and as a consequence we adopt the grid sampling.
In the text

	Fig. 6 Grid sample of the space for the first 3 observations(top panel) and for their remeasurement (bottom panel) of P10vxCt.
In the text

	Fig. 7 Two batches of observations reported in Table 4, for NEOCP object P10vxCt. The red line represents the originally submitted tracklet, whereas the blue line the remeasured tracklet. The higher curvature of the first arc with respect to the second is clear.
In the text

	Fig. 8 Admissible region and grid sample of the space for the Gaia object g0T015. Top panel: Grid sample of the space for the Gaia object g0T015. Bottom panel: Zoom-in of the second grid is shown.
In the text

	Fig. 9 Residuals in right ascension and declination. Values in arcsec were 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 g0T015 (2017 AD17).
In the text

	Fig. 10 Admissible region and grid sample of the space for the Gaia object g1j0D7. Top panel: Grid sample of the space for the Gaia object g1j0D7. Bottom panel: Zoom-in of the second grid is shown.
In the text

	Fig. 11 Residuals in right ascension and declination. Values in arcsec were 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 g1j0D7.
In the text