A systematic ranging technique for follow-ups of NEOs detected with the Flyeye telescope

TL;DR

This paper presents a systematic ranging method for rapidly determining the probable sky region of newly detected near-Earth objects by ESA's Flyeye telescope, aiding immediate follow-up observations and impact risk assessment.

Contribution

It introduces a systematic ranging approach tailored for short-arc orbit determination of NEOs detected by Flyeye, enhancing impact prediction capabilities.

Findings

Effective prediction of sky regions for follow-up observations.

Quantitative detection probability estimates.

Assessment of suitable follow-up telescopes.

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 near-Earth 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 follow-up observations and investigate potential follow-up telescopes for the Flyeye telescope.

Figures

Tables

Table 1. Table 1: Astrometric standard deviation for a number of important observatories as used in [ 8 ] . We use σ α , δ = 1 ​ arcsec subscript 𝜎 𝛼 𝛿 1 arcsec \sigma_{\alpha,\delta}=1\,\mathrm{arcsec} for non-listed observatories.

Obs. code	${\sigma }_{\alpha ,\delta }$ [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

Table 2. Table 2: Set of follow-up telescopes with their field of view (FOV).

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 Test-Bed 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

Equations

$${\mathbf{\nu }}^{T}W\mathbf{\nu }=\text{minimum}$$

$$\Delta \mathbf{A}=-(B^{T}WB{)}^{-1}{B}^{T}W\mathbf{\nu }$$

$$\mathrm{RMS}=\sqrt{\frac{1}{2n}{\mathbf{\nu }}^{T}W\mathbf{\nu }}$$

$$f_{\mathrm{post}}\propto f_{\mathrm{err}}{f}_{\mathrm{prior}}$$

$$f_{\mathrm{err}}=\exp \left(-\frac{1}{2}{\mathbf{\nu }}^{T}W\mathbf{\nu }\right)$$

$$f_{\mathrm{prior}}={\varrho }^{2}f(H)={\varrho }^{2-5\eta }$$

$$\displaystyle H=H(\varphi)+2.5\log\big{(}(1-G)\Phi_{1}(\varphi)+G\Phi_{2}(\varphi)\big{)}$$

$$H(\phi )=m-5{\log }_{10}(d\varrho )$$

$${\Phi }_i=W{\Phi }_{i\mathrm{S}}+(1-W){\Phi }_{i\mathrm{L}}$$

$$W=\exp \left(-90.56\tan (\phi /2{)}^2\right)$$

$${\Phi }_{i\mathrm{S}}=1-\frac{C_i\sin (\phi )}{0.119+1.341\sin (\phi )-0.754\sin (\phi {)}^2}$$

$${\Phi }_{i\mathrm{L}}=\exp \left(-A_i\tan (\phi /2{)}^{B_i}\right)$$

$$X=-\frac{\cos (\delta )\sin (\alpha -{\alpha }_{\mathrm{f}})}{\cos ({\delta }_{\mathrm{f}})\cos (\delta )\cos (\alpha -{\alpha }_{\mathrm{f}})+\sin ({\delta }_{\mathrm{f}})\sin (\delta )}$$

$$Y=-\frac{\sin ({\delta }_{\mathrm{f}})\cos (\delta )\cos (\alpha -{\alpha }_{\mathrm{f}})-\cos ({\delta }_{\mathrm{f}})\sin (\delta )}{\cos ({\delta }_{\mathrm{f}})\cos (\delta )\cos (\alpha -{\alpha }_{\mathrm{f}})+\sin ({\delta }_{\mathrm{f}})\sin (\delta )}$$

Full text

1 Introduction

In recent years, many near-Earth 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 website111http://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 newsletter1 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 Pan-STARRS [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 ($6.7$ deg by $6.7$ 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 $21.5$ mag. When a new object is found, immediate orbit determination needs to be done to filter potential impactors from passing NEOs and main-belt 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 follow-up observations are obtained. However, without a good knowledge of the orbit, follow-up 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 follow-up stations of the Flyeye telescope.

2 Systematic ranging

For a newly discovered object, we have a set of $n$ observations, the so-called tracklet, consisting of right ascensions $\alpha$, declinations $\delta$ and potentially apparent magnitudes $m$. With this data, rates of right ascension $\dot{\alpha}$ and rates of declination $\dot{\delta}$ can be determined, leading to the so-called attributable $\mathbf{A}=(\alpha,\delta,\dot{\alpha},\dot{\delta})$ [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 $\varrho$, the distance between the observer station and the observed object, and the range rate $\dot{\varrho}$ as intuitive 5th and 6th parameter. The lack of information is overcome by systematically scanning a ($\varrho$,$\dot{\varrho}$)-grid, in contrast to the Monte Carlo based statistical ranging [9], where $\varrho$ and $\dot{\varrho}$ are randomly sampled. For each grid point with fixed $\varrho_{i}$ and $\dot{\varrho}_{j}$, a best-fit attributable $\mathbf{A}_{ij}(\varrho_{i},\dot{\varrho}_{j})$ for the first observation epoch and a corresponding orbit can be determined by minimizing the cost function [10]:

[TABLE]

where $\boldsymbol{\nu}=(\alpha_{1}^{O}-\alpha_{1}^{C},\delta_{1}^{O}-\delta_{1}^{C},\ldots,\delta_{n}^{O}-\delta_{n}^{C})$ is the vector of the astrometric residuals of $n$ observations with the observed values $O$ and the computed values $C$. The weight matrix $W$ is a $(2n\times 2n)$ matrix describing the uncertainties of the various measurements. We use a weight matrix with the inverse standard errors $1/\sigma^{2}$ 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

[TABLE]

to A, where $B=\partial\boldsymbol{\nu}/\partial\mathbf{A}$ is a $(2n\times 4)$ matrix of measurement equation coefficients. An appropriate initial attributable can be estimated by the first and last observations $\mathbf{A}_{\mathrm{ini}}=(\alpha_{1},\delta_{1},(\alpha_{n}-\alpha_{1})/(t_{n}-t_{1}),(\delta_{n}-\delta_{1})/(t_{n}-t_{1}))$. Note, that we include a light-time correction by subtracting $\varrho/c$ from the observed epochs for any propagation computation with $c$ as speed of light.

For each grid point a weighted astrometric root mean square (RMS) error can be computed, indicating the quality of the orbit fit:

[TABLE]

An example can be seen in left plot of Figure 1 for the object with the temporary designation ZW04D95222https://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 best-fit solution. To validate our results we compare them with the contour map (red dashed contour lines) of JPL’s Scout website333https://cneos.jpl.nasa.gov/scout, provided by D. Farnocchia (personal communication, 2018). Both RMS errors are computed by using the same astrometric standard deviations $\sigma_{\alpha,\delta}$. 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 $f_{\mathrm{post}}$ for topocentric range and range rate, by multiplying the error function $f_{\mathrm{err}}$ of the normally distributed observation errors with a prior probability density function $f_{\mathrm{prior}}$, as given by [8]:

[TABLE]

with a power law size distribution $f(H)\propto\varrho^{-5\eta}$ and $\eta$ ranging between $0.35$ and $0.47$ [8]. We select for our computations $\eta=0.41$. For the probability density computation only heliocentric bound orbits with eccentricity $e<1$ 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 $f_{post}$ values, starting with the lowest probabilities, to 5% of the total grid sum, we determine a 95% confidence region of the ($\varrho,\dot{\varrho}$)-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 $\varrho_{k}$ and range rate $\dot{\varrho}_{k}$ within the grid limits; 2. 2.

Randomly picking a posterior probability density $\gamma_{k}$ between zero and the upper limit of the maximum computed value of the grid; 3. 3.

Computing the attributable $\mathbf{A}_{k}(\varrho_{k},\dot{\varrho}_{k})$ and posterior probability density $f_{\mathrm{post},k}(\varrho_{k},\dot{\varrho}_{k})$; 4. 4.

Continuing if $f_{\mathrm{post},k}\geq\gamma_{k}$, otherwise restarting with step 1; 5. 5.

Adding Gaussian noise, defined by the covariance matrix $\Gamma_{k}=({B_{k}}^{T}WB_{k})^{-1}$, to $\mathbf{A}_{k}(\varrho_{k},\dot{\varrho}_{k})$;

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 $\varrho$. 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 follow-up observations

We investigate the change of RMS error and 95% confidence region after follow-up 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 Pan-STARRS 1 (F51), spanning a time period of $37$ min. The second tracklet includes 3 follow-ups taken from Mauna Kea (568), leading to a total of 6 observations covering $47$ h. Finally we use 9 observations within $49$ h, including data from the Steward Observatory (691). The result can be seen in Figure 2, where the red dash-dotted 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 follow-up pointing

We want to find a promising right ascension $\alpha_{\mathrm{f}}$ and declination $\delta_{\mathrm{f}}$ to point at for follow-up observations. To this aim, we propagate the samples to a certain epoch.

A very simple approach is to compute the arithmetic means of $\alpha$ and $\delta$ 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 $6$ 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 $47$ 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 $12$ h instead of $6$ 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 $\alpha$ and $\delta$, 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 $6$ h and $12$ 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 $6$ 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 Gaussian-distributed, 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 Follow-up pointing for chosen telescope

If a specific telescope is chosen for the follow-up observations, not all samples might be detectable as,

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

depending on the telescope latitude, longitude and follow-up observation epoch, some of the samples could have $\alpha$ and $\delta$, which are below the horizon.

Those undetectable samples are not considered in determining the median values, but only the detectable ones.

The absolute magnitude $H$ in the H-G magnitude system [12] can be computed from the apparent magnitude $m$ by

[TABLE]

and vice versa, where $i=1,2$. The function $H(\varphi)$ is called reduced magnitude with the phase functions $\Phi_{i}$. $\varphi$ is the phase angle and $d$ corresponds to the heliocentric distance of the sample. The corresponding constants are $A_{1}=3.332$, $B_{1}=0.631$, $C_{1}=0.986$, $A_{2}=1.862$, $B_{2}=1.218$ and $C_{2}=0.238$. As we do not have observed phase curves for newly discovered NEOs, we assume a slope parameter $G=0.15$ for all objects.

Similar to the attributable fit of Equation 1, we compute the best fit $H_{k}$ from the photometric data by a least squares fit for each sample, where we consider the previously found best fit attributable. Finally the $m_{k}$ of the samples depend on the propagation epoch.

An example of the $H$ computation can be seen in the upper right plot of Figure 6, where we computed the absolute magnitude of ZW04D95 as function of $\varrho$ and $\dot{\varrho}$.

7 Detection probability for a chosen telescope

To compute the number of detected samples by a follow-up observation, we map the samples with their right ascension $\alpha$ and declination $\delta$ on the field of view of a specific telescope by [13]:

[TABLE]

where $\alpha_{\mathrm{f}}$ corresponds to the right ascension center of the follow-up observation and $\delta_{\mathrm{f}}$ the declination center. If a sample lies within the boundaries of the FOV by checking $|X|\leq\arctan(\mathrm{FOV}/2)$ and $|Y|\leq\arctan(\mathrm{FOV}/2)$, 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 $p$ 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.

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 $47$ arcmin as follow-up observatory and $20$ h after the last observation as observation epoch. Instead of simply considering the horizon as observation limit, we set the limit to $30$ 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 $m=18.5$. 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 $p(20\;\mathrm{h},\mathrm{J04})\approx 0.71$.

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 $p$, 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, non-station 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 dash-dotted) is slow, ZW04D95 (blue dashed) is intermediate and ZTF027b (red solid) is fast. As expected, $p$ 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 $6$ h, which is why even after $12$ h the spread has not a large effect on the detection probability for FOVs of the order of $1$ arcmin. A clear difference can be seen for ZW04D95, where $p$ 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 $p<0.5$. For FOVs of the order of $10$ 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 follow-ups of those objects. To this aim, a set of data points representing several follow-up 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 follow-up 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 ZW04D95-like objects is very high. In contrast, ZTF027b-like objects are a bigger challenge as there is only a single considered station with $p>0.9$, the Test-Bed 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 follow-ups. For objects like ZTF027b, Z58 will be one of the important follow-up observatories of the Flyeye telescope. However, the epoch of the real follow-up observation is not strictly limited to the longitude difference of the locations, but might be much earlier. This would lead to better $p$ 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 follow-ups 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 best-fitting 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 follow-up.

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 follow-ups 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 follow-up, 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 Test-Bed Telescope (Z58) in Cebreros, Spain, will have the important task to do the follow-ups 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.