A systematic ranging technique for followups of NEOs detected with the Flyeye telescope
Abstract
When new objects are detected in the sky, an orbit determination needs to be performed immediately to find out their origin, to determine the probability of an Earth impact and possibly also to estimate the impact region on Earth. ESA’s Flyeye telescope is expected to revolutionize the effort of predicting potential asteroid or deep space debris impact hazards due to the expected increase of nearEarth object discoveries. As the observed orbit arc for such an object is very short, classical Gaussian orbit determination cannot be used. We adopt the systematic ranging technique to overcome the lack of information and predict a region of the sky where the body can most likely be found. We also provide a detection probability for followup observations and investigate potential followup telescopes for the Flyeye telescope.
1 Introduction
In recent years, many nearEarth objects (NEOs) have been detected, leading to a current (2019, Jan.) number of more than 19 300 known objects, listed on ESA’s NEO Coordination Centre website^{1}^{1}1http://neo.ssa.esa.int. Due to an increasing number of survey telescopes, the discovery rate has risen over the past decades, reaching a value of 2033 objects in 2017 and 1822 objects in 2018. The monthly NEOCC newsletter^{1} provides regular updates of the newly detected objects. The lower rate of 2018 can be explained by bad weather conditions in Hawaii and an unplanned downtime of PanSTARRS [1], showing that the discovery of NEOs shall not rely on a single observatory, but on several independent telescopes, spread all over the globe.
ESA will contribute to the task of searching for uncatalogued NEOs with the Flyeye telescope, a telescope splitting the large field of view ( deg by deg) by use of a prismatic faceted mirror onto 16 individual cameras [2]. The telescope will allow the coverage of almost a full hemisphere every two observing nights. Each area of the sky is visited four times, collecting information on the possible movement of objects brighter than mag. When a new object is found, immediate orbit determination needs to be done to filter potential impactors from passing NEOs and mainbelt asteroids, estimating the impact probability and a potential impact region on Earth. Examples of impactors are the asteroids 2008 TC3, 2014 AA and 2018 LA, which all were discovered by R. A. Kowalski as part of the Mt. Lemmon Survey (G96) and impacted less than a day after they were detected [3, 4, 5].
For a newly discovered object, the observed orbit arc is usually very short [6], with just four reliable observables: right ascension, declination, and the rate of change of right ascension and declination. For this reason, classical Gaussian orbit determination cannot be applied until further followup observations are obtained. However, without a good knowledge of the orbit, followup observations might point to wrong celestial coordinates and miss the object, leading to a vicious circle. In order to find the object a few or even many hours later, it is necessary to put some constraints on the region of the sky to observe. For this goal, we are implementing the technique of systematic ranging, described in Chapter 2, in our system, where a raster scan in the space of topocentric range and range rate is performed. By assigning a probability density to this space, we generate in Chapter 5 a sample cloud on the sky, which can be used by observers to select a telescope pointing. In Chapter 6 we calculate a promising pointing for a specific telescope that increases the detection probability. Furthermore we investigate the detection probabilities for a number of collaborating followup stations of the Flyeye telescope.
2 Systematic ranging
For a newly discovered object, we have a set of observations, the socalled tracklet, consisting of right ascensions , declinations and potentially apparent magnitudes . With this data, rates of right ascension and rates of declination can be determined, leading to the socalled attributable [6]. However, an orbit is fully defined by 6 parameters, which means two additional parameters have to be guessed.
We follow the systematic ranging technique as described in [7] and [8], using the topocentric range , the distance between the observer station and the observed object, and the range rate as intuitive 5th and 6th parameter. The lack of information is overcome by systematically scanning a (,)grid, in contrast to the Monte Carlo based statistical ranging [9], where and are randomly sampled. For each grid point with fixed and , a bestfit attributable for the first observation epoch and a corresponding orbit can be determined by minimizing the cost function [10]:
(1) 
where is the vector of the astrometric residuals of observations with the observed values and the computed values . The weight matrix is a matrix describing the uncertainties of the various measurements. We use a weight matrix with the inverse standard errors as diagonal elements. The standard deviations are given in Table 1 [8]. An iterative solution with the least squares fit method can be achieved by adding
(2) 
to A, where is a matrix of measurement equation coefficients. An appropriate initial attributable can be estimated by the first and last observations . Note, that we include a lighttime correction by subtracting from the observed epochs for any propagation computation with as speed of light.
Obs. code  [arcsec] 

568  0,15 
F51  0,20 
H01  0,30 
E10, F65, J04, W84  0,40 
291, 691, 950, G96  0,50 
W85, W86, W87  0,60 
H21  0,70 
G45, K91, K92, K93  0,80 
Q63, Q64, V37  0,80 
703  1.00 
For each grid point a weighted astrometric root mean square (RMS) error can be computed, indicating the quality of the orbit fit:
(3) 
An example can be seen in left plot of Figure 1 for the object with the temporary designation ZW04D95^{2}^{2}2https://www.minorplanetcenter.net. The black dotted line corresponds to the boundary between elliptic and hyperbolic orbits. The blue solid contour lines show the weighted astrometric RMS error of the bestfit solution. To validate our results we compare them with the contour map (red dashed contour lines) of JPL’s Scout website^{3}^{3}3https://cneos.jpl.nasa.gov/scout, provided by D. Farnocchia (personal communication, 2018). Both RMS errors are computed by using the same astrometric standard deviations . The two results match very well, except for approaching orbits close to Earth. In this region the gravity of Earth and Moon gets significant, yet for performance reasons we only use a Kepler propagator assuming unperturbed heliocentric orbits. Hence, this difference is expected. The perturbations are important for an accurate orbit determination, which is why we will include them in future work, though they can be neglected for the goal of this paper.
We compute a posterior probability density for topocentric range and range rate, by multiplying the error function of the normally distributed observation errors with a prior probability density function , as given by [8]:
(4)  
(5)  
(6) 
with a power law size distribution and ranging between and [8]. We select for our computations . For the probability density computation only heliocentric bound orbits with eccentricity are considered, while the probabilities for unbound orbits are set to zero. The probability density is plotted on the right side of Figure 1 as blue solid lines. By summing up the values, starting with the lowest probabilities, to 5% of the total grid sum, we determine a 95% confidence region of the ()space. This region is enclosed by the magenta dashed line.
3 Monte Carlo samples
To predict the distribution, where the object may be found at an arbitrary epoch, we generate Monte Carlo samples, following [8]:

Randomly drawing a topocentric range and range rate within the grid limits;

Randomly picking a posterior probability density between zero and the upper limit of the maximum computed value of the grid;

Computing the attributable and posterior probability density ;

Continuing if , otherwise restarting with step 1;

Adding Gaussian noise, defined by the covariance matrix , to ;
For a better resolution at small topocentric ranges, we sample the range in a logarithmic scale. In this case, the density function of Equation 4 has to be multiplied with . We generate 500 Monte Carlo samples in Figure 1, where as expected most of them are located inside the magenta dashed 95% confidence region. Starting from the generated topocentric ranges, range rates and computed attributables, the orbital elements can be derived and processed, e.g. to propagate the samples and check for impacts on Earth to estimate an impact probability.
4 Change of RMS error and 95% confidence region after followup observations
We investigate the change of RMS error and 95% confidence region after followup observations. To this aim, we repeat the computations of the object P10LiAR with different number of observations. The first tracklet consists of 3 observations, the initial discovery, by the station PanSTARRS 1 (F51), spanning a time period of min. The second tracklet includes 3 followups taken from Mauna Kea (568), leading to a total of 6 observations covering h. Finally we use 9 observations within h, including data from the Steward Observatory (691). The result can be seen in Figure 2, where the red dashdotted lines corresponds to 3 observations, the blue dashed lines to 6 observations and the green solid lines to 9 observations. The left plot shows the weighted astrometric RMS error and the right plot shows the 95% confidence region. One can see that the region with a good fit of the observations shrinks with more data taken over a larger time span and from independent data sources. As a consequence, the 95% confidence region decreases too.
5 Promising followup pointing
We want to find a promising right ascension and declination to point at for followup observations. To this aim, we propagate the samples to a certain epoch.
A very simple approach is to compute the arithmetic means of and of the propagated samples. This can be seen in Figure 3 as a red cross, where we propagated 500 random samples (black dots) of P10LfOv, ZW04D95 and ZTF027b to h after their last observation. To keep the computation general, we use the geocenter as reference for the computation instead of a specific station. An example field of view (FOV) of arcmin, corresponding to ESA’s Optical Ground Station (J04), is shown as dashed frame around the cross. While for P10LfOv the sample cloud is centered very well, the mean of ZW04D95 is already shifted to a position where only 2 out of 500 samples are located in the FOV. This shift is the result of a few large outliers, which occur for samples that already came close to Earth and hence moved to the other side of the sky. The shape of the ZTF027b samples is not a roundish cloud anymore but very stretched and quite a few outliers can be seen. This shape arises from the computed posterior probability distribution and the resulting samples being placed in orbits closer to Earth than the samples of the other objects, which is shown in Figure 5. As a result of the different ranges, the apparent velocities of the samples differ too, where P10LfOv has the lowest apparent velocities and ZTF027b the highest. Those high apparent velocities lead to the spread of the cloud. For ZTF027b, the mean computation would suggest a bad pointing far off the sample line. If the samples are propagated h instead of h, as shown in Figure 4 for ZW04D95, the cloud spreads too, leading to an even more unfavorable pointing recommendation.
Another approach is to compute the median of all and , which can be seen as a blue cross with solid frame in Figure 3 and Figure 4. Similar to the mean, the median centers the small cloud of P10LfOv very well, why one cannot see a notable difference of the FOVs. In contrast, the median focuses very well on the h and h propagated samples of ZW04D95. For ZTF027b, the median places the pointing recommendation on the sample line, yet it centers not necessarily on the densest region. However, since the samples are widely spread and the sample density inside the FOV is low anyways, not finding a promising pointing for objects like ZTF027b will not make a big difference for the later computation of the detection probability. To rediscover ZTF027b after h, an observation campaign or a piece of good luck is needed.
It can be seen that using the mean results in physically unrealistic solutions, the median gives better results. This is due to the fact that errors are not Gaussiandistributed, similar to what was shown in [11]. Therefore, for the rest of the paper, we will use the median to compute the promising pointing.
6 Followup pointing for chosen telescope
If a specific telescope is chosen for the followup observations, not all samples might be detectable as,

the topocentric range differs for each sample and the resulting apparent magnitude might be too faint for a specific telescope with a given limiting visual magnitude.

depending on the telescope latitude, longitude and followup observation epoch, some of the samples could have and , which are below the horizon.
Those undetectable samples are not considered in determining the median values, but only the detectable ones.
The absolute magnitude in the HG magnitude system [12] can be computed from the apparent magnitude by
(7)  
(8)  
(9)  
(10)  
(11)  
(12) 
and vice versa, where . The function is called reduced magnitude with the phase functions . is the phase angle and corresponds to the heliocentric distance of the sample. The corresponding constants are , , , , and . As we do not have observed phase curves for newly discovered NEOs, we assume a slope parameter for all objects.
Similar to the attributable fit of Equation 1, we compute the best fit from the photometric data by a least squares fit for each sample, where we consider the previously found best fit attributable. Finally the of the samples depend on the propagation epoch.
An example of the computation can be seen in the upper right plot of Figure 6, where we computed the absolute magnitude of ZW04D95 as function of and .
7 Detection probability for a chosen telescope
To compute the number of detected samples by a followup observation, we map the samples with their right ascension and declination on the field of view of a specific telescope by [13]:
(13)  
(14) 
where corresponds to the right ascension center of the followup observation and the declination center. If a sample lies within the boundaries of the FOV by checking and , the sample is found. Depending on the specific telescope and CCD properties, further mapping, e.g. due to a rotated CCD, must be implemented before the check.
Only the visible samples, which are used to determine the center of observation, have to be accounted for. Dividing the sum of found samples by the full sample number, including the too faint and below horizon objects, leads to an overall detection probability of a single observation by a specific telescope at a certain epoch. To get a reliable result, a large enough Monte Carlo sample sequence has to be used.
Obs. code  Telescope name  FOV [arcmin] 

J04  1.0 m Optical Ground Station at Teide  47 
Z84  0.8 m Schmidt telescope  17 
Z58  0.56 m TestBed Telescope  149 
309  8.2 m Very Large Telescope  7 
344  1.8 m BOAO telescope  15 
345  0.6 m SOAO telescope  15 
Y28  1.0 m OASI telescope  10 
679  2.12 m telescope at San Pedro Martir  6 
An example can be seen in Figure 6, where the promising pointing for the object ZW04D95 is shown in the lower plots. We chose ESA’s Optical Ground Station (J04) with a FOV of arcmin as followup observatory and h after the last observation as observation epoch. Instead of simply considering the horizon as observation limit, we set the limit to deg elevation (cyan solid line), where one would expect good observation data. The zenith of the station is indicated by the cyan X and the horizon by the cyan dashed line in the lower left plot. For demonstration reasons, the limiting magnitude is set to . The blue cross marks the computed promising pointing, surrounded by the projected FOV. Green dots show detectable objects while red diamonds are below the horizon and yellow triangles are too faint. Out of 500 samples, 356 observable ones match with the FOV projection, leading to a detection probability of .
8 Detection probability evolution
As we have already noticed in Chapter 5, the sample cloud spreads with time. To investigate this effect on the detection probability, we neglect the other influences on , i.e. the horizon and the limiting magnitude, which both depend on the observatory. The propagated right ascensions and declinations are once more computed from the geocenter, to obtain a general, nonstation dependent solution.
The detection probability of a certain object depends now on the time after the last observation and the field of view, which is shown in Figure 7 for the same 3 examples as in Chapter 5. The main difference between the observed asteroids is their apparent velocity: P10LfOv (green dashdotted) is slow, ZW04D95 (blue dashed) is intermediate and ZTF027b (red solid) is fast. As expected, shrinks with increasing time, as the sample cloud spreads, and with decreasing FOV. As observed before in Figure 3, we know that the sample cloud of P10LfOv is very narrow after h, which is why even after h the spread has not a large effect on the detection probability for FOVs of the order of arcmin. A clear difference can be seen for ZW04D95, where ranges from less than 0.1 to more than 0.9 in the given boundaries. The most challenging example is ZTF027b, where the detection probability decreases very quickly within 10 h even for large FOVs of several hundreds of arcmin to . For FOVs of the order of arcmin, a probability of less than 10% is already reached within a few hours.
As we expect the Flyeye telescope to find many new NEOs, we investigate the chances of successful followups of those objects. To this aim, a set of data points representing several followup telescopes, given in Table 2, is placed in Figure 7. We distinguish between stations of the northern hemisphere (black cross) and southern hemisphere (magenta X). The abscissa represents the longitude difference between the future Flyeye telescope longitude at Monte Mufara and the followup telescope longitude in hours. Therefore if an object similar to the given examples is detected by the Flyeye telescope, we can estimate the detection probability by certain stations at their location.
Among our chosen set of stations, the likelihood of observing P10LfOv and ZW04D95like objects is very high. In contrast, ZTF027blike objects are a bigger challenge as there is only a single considered station with , the TestBed Telescope (Z58). Two more observatories of the northern hemisphere have chances close to 50%. The other stations are clearly below 10%. Based on Figure 7 we can estimate the chance of observing an object and decide for a station with high enough detection probability to do followups. For objects like ZTF027b, Z58 will be one of the important followup observatories of the Flyeye telescope. However, the epoch of the real followup observation is not strictly limited to the longitude difference of the locations, but might be much earlier. This would lead to better for the stations.
Please note that the object might be below the horizon due to not ideal station latitude, which is only marginally considered by distinguishing between northern and southern hemisphere stations. As the Flyeye telescope will be located at Monte Mufara, Italy, telescopes of the northern hemisphere will have to do the early followups to keep track of the new NEOs.
9 Conclusions
We created a software, based on the systematic ranging technique, scanning the space of topocentric range and range rate of a tracklet for the bestfitting orbits. Using a prior probability density function, we get a probability distribution, which is used to create a set of Monte Carlo generated samples. Those samples can be propagated to an arbitrary epoch, leading to a right ascension and declination distribution, which can be used by observers of any station to select a preferred telescope pointing to find the object. Finally, using the field of view of a chosen telescope, the detection probability for a pointing to the median of the Monte Carlo samples. With this information, observers quickly get an impression of the success chance of the scheduled followup.
We also investigated the evolution of the weighted astrometric RMS error and 95% confidence regions for a typical NEO detection. Both shrink with additional data, showing the importance of followups to have a longer arc for orbit determination. From Figure 3 and Figure 4, one can see a large difference of the sample clouds. In particular NEOs close to Earth have elongated, spread clouds due to the high apparent velocities. Those objects are a challenge to followup, as their detection probability decreases quickly with time, shown in Figure 7. Expecting similar tracklets from the Flyeye telescope in the future leads to the conclusion that a number of stations, most notably the TestBed Telescope (Z58) in Cebreros, Spain, will have the important task to do the followups of NEOs close to Earth, as they might be lost otherwise. In addition, such close objects could be impactors, emphasizing the urgency of observations to estimate a more precise impact probability and predict the potential impact locations. To this aim, the software must be fully automatized and will have to run immediately after the discovery of a new object, notifying observers of the discovery of a critical NEO.
Acknowledgments
We thank Davide Farnocchia for assistance with the systematic ranging technique and providing his data for validation.
This research has made use of data and/or services provided by the International Astronomical Union’s Minor Planet Center.
References
 [1] Ramanjooloo, Y., Wainscoat, R. J., Weryk, R., Chambers, K. C., Magnier, E. A. (2018). The PanSTARRS search for NearEarth Objects. In American Astronomical Society, DPS Meeting #50, id304.08.
 [2] Cibin, L., Chiarini, M., Gregori, P, Bernardi, F., Ragazzoni, R., Sessler, G., Kugel, U. (2019, Jan.). The flyeye telescope, development and first factory tests results. Paper presented at the ESA NEO and Debris Detection Conference, Darmstadt, Germany.
 [3] Minor Planet Center (2008). MPEC 2008T50 : 2008 TC3. Online at https://minorplanetcenter.net/mpec/K08/K08T50.html (as of 14 Jan. 2019)
 [4] Minor Planet Center (2014). MPEC 2014A02 : 2014 AA. Online at https://minorplanetcenter.net/mpec/K14/K14A02.html (as of 14 Jan. 2019)
 [5] Minor Planet Center (2018). MPEC 2018L04 : 2018 LA. Online at https://minorplanetcenter.net/mpec/K18/K18L04.html (as of 14 Jan. 2019)
 [6] Milani, A., Gronchi, G. F., Vitturi, M. D. M., Knežević, Z. (2004). Orbit determination with very short arcs. I admissible regions. Celestial Mechanics and Dynamical Astronomy, 90(12), 5785.
 [7] Chesley, S. R. (2004). Very short arc orbit determination: The case of asteroid 2004 FU. In Proceedings of the International Astronomical Union, 2004(IAUC197), pp255–258.
 [8] Farnocchia, D., Chesley, S. R., Micheli, M. (2015). Systematic ranging and late warning asteroid impacts. Icarus, 258, 18–27.
 [9] Virtanen, J., Muinonen, K., Bowell, E. (2001). Statistical ranging of asteroid orbits. Icarus, 154(2), 412–431.
 [10] Jehn, R. (1989). ORBDET: A program for orbit determination (MAS Working paper Nr. 306), European Space Operations Center, Darmstadt, Germany.
 [11] Schmidt, A. K. (2019). Analysis of uncertainties and algorithmic errors in a doublestation meteor camera setup. (Unpublished master’s thesis). University of Oldenburg, Oldenburg, Germany.
 [12] Bowell, E., Hapke, B., Domingue, D., Lumme, K., Peltoniemi, J., Harris, A. W. (1989). Application of photometric models to asteroids. In Asteroids II (Eds. R. P. Binzel, T. Gehrels & M. S. Matthews), Univ. of Arizona Press, Tucson, USA, pp524–556.
 [13] Montenbruck, O., Pfleger, T. (2000). Astronomy on the personal computer. SpringerVerlag Berlin Heidelberg, Berlin Heidelberg, Germany, pp251–266.