© ESO, 2015

1. Introduction

Highlighted in the early 20th century, the Yarkovsky effect is a non-gravitational force due to an anisotropic emission of thermal radiation. It mainly depends on the size, the density, the spin, and the thermal characteristics of the asteroid. The effect has a diurnal component related to the asteroid rotation, which is generally dominant, and a seasonal component related to the orbital motion (Vokrouhlický et al. 2000).

The Yarkovsky effect leads to an acceleration on the asteroid that modifies its semi-major axis as a secular drift. It has important consequences for the long-term evolution of asteroids. In particular, the Yarkovsky effect could explain the transportation of asteroids from the main belt region to the inner solar system (Bottke et al. 2006). It also becomes necessary for accurate orbit determination or prediction of trajectory of near-Earth asteroids (NEA), in particular, in the context of close approach or impact predictions (Farnocchia et al. 2013a; Farnocchia & Chesley 2014; Chesley et al. 2014).

The first detection of Yarkovsky acceleration has been performed thanks to radar measurements of the asteroid (6489) Golevka (Chesley et al. 2003). Vokrouhlický et al. (2008) measured the Yarkovsky effect using only astrometric observations of the asteroid (152563) 1992BF.

More recently, Nugent et al. (2012b) detected a measurable secular drift for 54 NEAs and Farnocchia et al. (2013b) for 21 NEAs using more restrictive criteria. Nugent et al. (2012b) used the Orbfit package¹ for orbit determination and Farnocchia et al. (2013b) used JPL’s Comet and Asteroid Orbit Determination Package and cross-checked their results with Orbfit.

In this paper, we perform an independent study of the detection of the drift in semi-major axis of NEAs. We develop our own dynamical model allowing the orbit determination and the measurement of a secular drift in semi-major axis (Sect. 2). We aim to detect a secular drift for the 10953 NEAs² (1629 numbered and 9324 unnumbered) from ground-based observations. We first verify the robustness and quality of the detection with statistical methods (Sects. 3.1 and 3.2) and finally, by introducing filters to avoid unreliable detections, we detect a consistent measurement for 46 NEAs (Sect. 3). As a long period of observations appears necessary for an accurate detection, we analyse how data mining through precovery observations can help to refine the measurement of a secular drift (Sect. 4). Finally, we highlight the expected improvement with the Gaia mission for the detection of Yarkovsky acceleration (Sect. 5).

2. Orbit determination

To perform the orbit determination and to detect the drift, we build a dynamical model of an asteroid’s motion taking the main perturbations into account (Sect. 2.1), and we add the Yarkovsky effect as a transverse force depending on the orbital elements and on a drift in semi-major axis (Sect. 2.2). The observations and their treatment are presented in Sect. 2.3. Finally, the fitting process, used to determine the state-vector of the asteroid (i.e. the position and the velocity of the asteroid at a specific date) and the value of the drift ȧ, is presented in Sect. 2.4.

2.1. Dynamical model

To measure a drift in semi-major axis of asteroids, we developed a numerical integration of the motion of an asteroid (NIMA). The equations of the motion are numerically integrated with a Gauss-Radau method (Everhart 1985). The dynamical model includes:

• the perturbations of the planets, Pluto and the Moon, whose positions are given by JPL Ephemeris DE421 (Folkner et al. 2009);

• the perturbations of the ten main asteroids whose positions have been preliminary computed by NIMA and turned into Chebychev ephemeris to make quick computation. The masses of the main asteroids are provided in Table 1;

• the corrections of relativistic effects of Sun, Jupiter, and the Earth as described in Moyer (2000);

• the Yarkovsky effect (see Sect. 2.2).

All solar system objects in the dynamical model are considered point masses. The equations of variation are also numerically integrated with the equations of motion to fit the model to observations by a Levenberg-Marquardt algorithm (see Sect. 2.4). To be consistent with the planetary ephemeris, we use the masses of the asteroids provided by Folkner et al. (2009) and given in Table 1. The masses of some big asteroids are not well determined. For example, we noticed that the mass of (704) Interamnia is 4.751 km³ s² in Farnocchia et al. (2013b) whereas the mass is 2.464 km³ s² in Folkner et al. (2009). Moreover, we compared our results for the 46 selected NEAs (see Sect. 3.4) by including the 25 biggest asteroids from Folkner et al. (2009) as well as the 10 biggest asteroids, and we noticed only minor differences. Consequently, and also for time computation, we choose the ten biggest asteroids corresponding to a mass greater than 10^-11M_⊙.

2.2. Yarkovsky acceleration modelling

We consider the Yarkovsky effect a comet-like model as presented in Marsden et al. (1973), and we use the effect for non-gravitational perturbation. The effect is considered as a transverse force: $\begin{array}{ccc}{{a}}_{\mathit{Y}}\mathrm{=}{\mathit{A}}_{\mathrm{2}}\mathit{g}\mathrm{\left(}\mathit{r}\mathrm{\right)}{t}& & \end{array}$(1)where t is the normalised transverse vector (r × v) × r with r and v the heliocentric position and the velocity of the asteroid, respectively.

The expression g(r) is an r-function developed in Marsden et al. (1973) and A₂ is a parameter associated with this force. As in Farnocchia et al. (2013b), we used a simple expression with g(r) = (r₀/r)² where r₀ = 1 au is a normalising parameter.

The Gauss equations give the variation of the semi-major axis related to the transverse acceleration $\begin{array}{ccc}\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathit{ȧ}\mathrm{=}\frac{\mathrm{2}\mathit{a}\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}{\mathit{A}}_{\mathrm{2}}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{n}{\mathit{r}}^{\mathrm{3}}}& & \end{array}$where a is the semi-major axis, e the eccentricity, n the mean motion, and r is the heliocentric distance of the asteroid.

Consequently, by averaging the previous relation on an orbital period, the mean variation of ȧ is $\begin{array}{ccc}\mathrm{⟨}\mathit{ȧ}\mathrm{⟩}& \mathrm{=}& \frac{\mathrm{2}{\mathit{A}}_{\mathrm{2}}{\mathit{r}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{n}{\mathit{a}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}\mathrm{·}\end{array}$Finally, the Yarkovsky acceleration is written as a transverse force only depending on the orbital elements and the drift in semi-major axis (now written as ȧ): $\begin{array}{ccc}{{a}}_{\mathit{Y}}\mathrm{=}\frac{\mathit{n}}{\mathrm{2}}\frac{{\mathit{a}}^{\mathrm{2}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\mathrm{)}}{{\mathit{r}}^{\mathrm{2}}}\mathrm{(}\mathit{ȧ}\mathrm{)}{t}\mathit{.}& & \end{array}$(2)This expression is equivalent to those given in Chesley et al. (2008) and Mouret & Mignard (2011), and is consistent with models in Farnocchia et al. (2013b) where they estimated the parameter A₂. We can estimate either A₂ parameter or the mean drift in semi-major axis ȧ. Typical values for the drift are 10^-3 to 10^-4 au/My (Vokrouhlický et al. 2000), corresponding to a range of 10^-13 to 10^-15 au/d² for A₂. Thereafter we use 10^-4 au/My as a unit for the drift in semi-major axis ȧ and its uncertainty σ_ȧ, and 10^-15 au/d² for A₂ parameter and its uncertainty.

2.3. Observations and treatement

The asteroid observations come from the Minor Planet Center database³ as of 12 May 2014. Observations are also available on the AstDyS database⁴. Observations suffer from several sources of systematic errors that can be crucial for asteroid orbit determination (Carpino et al. 2003). Zonal errors of stellar catalogues are one of the systematic errors for observations. Chesley et al. (2010) proposed a method to remove the bias for several stellar catalogues. Bias corrections are given for each observation in the AstDyS database. Finally, the outliers are eliminated thanks to the procedure described in Carpino et al. (2003).

2.4. Fitting and weighting processes

The fitting process consists in the determination of seven parameters (the six components of the position and velocity at a given date and the drift in semi-major axis ȧ or the A₂ parameter defined in Sect. 2.2) that minimises the residuals (observed positions minus computed positions). This determination uses a Levenberg-Marquardt algorithm by iterative corrections of each parameter. For each iteration, the corrections to apply are determined thanks to the equations of variation and the least square method (LSM; for more information, see for example Desmars et al. 2009). In the LSM, a weighting matrix V_obs is required and usually considered as a diagonal matrix where the diagonal components are ${\mathit{ϵ}}_{\mathit{i}}^{\mathrm{2}}\mathrm{=}\mathrm{1}\mathit{/}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}$ where ${\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}$ is the estimated variance of the observation i.

In orbit determination, the main difficulty is to give an appropriate weight to each observation. Carpino et al. (2003) and Chesley et al. (2010) discussed weighting schemes and showed the orbit determination improvement by weighting observations, according to the observatory and the stellar catalogue used for the reduction.

The AstDyS database provides the weight for each observation. As a general rule, we use the value given in AstDyS, except for specific cases. As in Farnocchia et al. (2013b), we consider that observations before 1950 have a weight of 10 arcsec, and observations from 1950 to 1990 a weight of 3 arcsec. These observations are mainly photographic plates affected by low sensitivity and reduced with less accurate stellar catalogue. The other case is when observations have been carefully reduced again, as in Vokrouhlický et al. (2008), with four precovery observations of (152563) 1992BF in January 1953.

3. Detection of the drift

In this section, we try to detect a drift in the semi-major axis of NEAs. We first study the consistency of the drift and its uncertainty determined by LSM by comparing with statistical methods for three representative NEAs (Sect. 3.1). We also look after the robustness of the drift we determined (Sect. 3.2) and we present the drift in semi-major axis for NEAs (Sects. 3.3 and 3.4) and a comparison with the previous studies (Sect. 3.5). Finally, we discuss the constrains in some physical parameters of the asteroids given by the drift (Sect. 3.6).

3.1. Estimation of uncertainty

The fitting process gives the orbital elements (more specifically the state vector) of the asteroid and the drift in the semi-major axis. The LSM also provides the covariance matrix Λ and the standard deviation of each parameter is given by the root square of diagonal elements. Thus, for the drift in semi-major axis, the standard deviation is given by $\sqrt{{\mathrm{\Lambda }}_{\mathrm{77}}}$. This can be considered as a measure of the uncertainty of the drift value σ_ȧ. According to LSM theory, the uncertainty of the drift only depends on the date of observations and the weights given to these observations. In particular, it does not depend on the observation coordinates whereas the drift value depends on the date, the weights, and the coordinates of the observations. This aspect is important because we can study the evolution of the drift uncertainty related to the precision (or the weight) of the observations without the knowledge of the coordinates themselves and thus we can determine the drift uncertainty with simulated observations (see Sects. 4 and 5).

To verify if this measure is a good estimation of the uncertainty, we compare σ_ȧ with accuracies obtained by two statistical methods. In this context, we apply a Monte Carlo process on the observations and a bootstrap re-sampling as described in Desmars et al. (2009). The Monte Carlo method on the observations (MCO) comes from a technique developed by Virtanen et al. (2001) and consists of adding a random noise on each observation and then fitting the dynamical model to the new set of observations, giving a new drift value. The bootstrap re-sampling (BR) was developed by Efron (1979) and consists in creating a new set of observations by sampling the original set of observations with a draw with replacement. The new set of observations has the same number of elements as the original set. Consequently, some observations can appear several times whereas others do not appear at all. The dynamical model is then fitted to the new set, giving a new drift value.

For both methods, the process is repeated many times and a measure of the uncertainty is given by the standard deviation of all the drift values. We consider three significative NEAs:

• 1.

  (2340) Hathor for which we determined ȧ = (−14.312 ± 3.297) × 10^-4 au/My, represents NEAs whose S/N is greater than 3;

• 2.

  (138852) 2000WN10 with ȧ = (15.025 ± 7.159) × 10^-4 au/My, represents NEAs with an average S/N;

• 3.

  (153792) 2001VH75 with ȧ = (12.812 ± 165.373) × 10^-4 au/My, represents NEAs with very bad detection of the drift.

For both methods (MCO and BR), we generate 1000 samples. Table 2 presents the value of the drift and the standard deviation obtained with the fitting process (LSM) and the two different methods (MCO and BR).

The estimations of the uncertainty of the drift are equivalent for LSM and MCO. With the bootstrap method, we have slight differences in particular for NEAs with a bad detection of the drift. Actually, LSM and MCO make use of equivalent assumptions (independence of observations, Gaussian errors of observations) whereas the BR makes use of minimal assumptions (independence of observations only). For NEAs with a good S/N, which is one of our criterions of selection (see Sect. 3.4,) the three methods give similar results. The standard deviation σ_ȧ determined with LSM is a good estimator of the drift in semi-major axis uncertainty and can be used as a measure of the uncertainty.

3.2. Solution stability

The stability of the solution can be tested with a random holdout method. It consists of randomly extracting a percentage of the observations, fitting the model with these observations, and repeating the process several times. The changes in the fitted parameter gives information on the stability. In this section, we use the same representative NEAs as in Sect. 3.1 and we randomly remove 10% of the initial set of observations. The process is repeated 1000 times. Table 3 presents the mean (μ) and the standard deviation (σ) of the drift in semi-major axis ȧ and its uncertainty σ_ȧ, obtained with random holdout method for (2340) Hathor, (138852) 2000WN10, and (153792) 2001VH75.

The means for ȧ and σ_ȧ have to be compared with values given in Table 2 with LSM. By removing 10% of the initial set, the values obtained for the drift and its uncertainty do not vary too much. It shows that the drift in semi-major axis ȧ we detected and its uncertainty σ_ȧ are robust and stable, even for inaccurate determination of ȧ.

3.3. Drift determination for numbered NEAs

In this section, we try to determine a drift in semi-major axis for the 1629 numbered NEAs as they have been observed over several oppositions. Figure 1 represents the length of the observed arc, the uncertainty of the drift in semi-major axis, and the number of observations for the 1629 numbered NEAs.

Fig. 1 Observed arc, uncertainty of the drift in semi-major axis (in au/My), and absolute magnitude H for the 1629 numbered NEAs. NEAs with radar measurements are represented with crosses and NEAs without radar measurements with bullets.

For most NEAs, the uncertainty of the drift we detected is greater than 10^-3 au/My, leading to a bad S/N and inconsistent detections. As the order of magnitude for the drift is about 10^-3−10^-4 au/My (Vokrouhlický et al. 2000), an uncertainty of the same order is at least required for a consistent determination of the drift.

Figure 1 shows that a long period of observations (at least 20 years) is necessary to have a small uncertainty of the drift. Nevertheless, with radar measurements, the period of observations can be smaller because they lead to important constraints in the orbit determination (see for example Desmars et al. 2013) and in particular for the detection of Yarkovsky acceleration (see the case of (6489) Golevka in Chesley et al. 2003). As the length of orbital arc is longer for bright NEAs, we do not notice a clear correlation between H magnitude and drift uncertainty.

3.4. Drift determination for selected NEAs

We try to detect a drift in semi-major axis for the 10953 NEAs (1629 numbered and 9324 unnumbered) as of 12 May 2014. To select reliable detections of drift in semi-major axis, we adopt two filters.

The first filter is the S/N defined as S/N = | ȧ | /σ_ȧ. Farnocchia et al. (2013b) used S/N ≥ 3 as a threshold to define reliable detections but they also considered lower values for some specific cases. This threshold appears restrictive because a detection with S/N ~ 2 also leads to consistent drifts. In contrast, Nugent et al. (2012b) used S/N ≥ 1 ,which seems to be a nonrestrictive threshold. If we assume a Gaussian distribution for the drift value and if S/N = 1, the probability of detecting a drift value with an opposite sign from the real value (for example a positive drift where the real drift is negative or, conversely, a negative drift where the real drift is positive) is about 15.9%, whereas if S/N = 2, this probability is 2.3%. In that context, we adopt S/N ≥ 2 as an appropriate filter.

The second filter is to avoid spurious detections: for example, the detection of a strong drift for large NEAs. We estimated a maximum value for the drift in the semi-major axis according to the absolute magnitude H of the NEA, $\begin{array}{ccc}\mathit{ȧ}\max \mathrm{=}{\left|{\left(\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}\right)}_{\mathit{d}}\right|}_{\max }\mathrm{=}{\mathrm{10}}^{\mathrm{0.2}\mathit{H}\mathrm{-}\mathrm{6.85}}\mathrm{au}\mathrm{/}\mathrm{My}\mathit{.}& & \end{array}$(3)The maximum value assumes average values of semi-major axis, albedo, and bulk density for NEA according to Chesley et al. (2002). The full details of the estimation are given in Appendix A. This parameter remains an approximation and can be underestimated by a factor of 3 for NEA with low density and by a factor of 5 for NEAs with high eccentricity (Vokrouhlický et al. 2000). The ratio ξ = | ȧ | /ȧ_max can help to filter spurious detection. The parameter ξ is expected to be smaller than 1 but by considering the approximate value of the maximum drift, we select NEAs for which ξ ≤ 3.

Finally, Table 4 presents the selected NEAs, for which S/N ≥ 2 and ξ ≤ 3, with the drift in semi-major axis ȧ, its uncertainty σ_ȧ, the A₂ parameter and its uncertainty σ_A2, the S/N, the absolute magnitude H, the diameter D (from EARN database⁵), the semi-major axis (in au), the eccentricity, the ξ ratio, and the time span of observations. An additional asterisk indicates that radar measurements are available for this NEA. We also add three NEAs with σ_ȧ ≤ 10^-4 au/My because as discussed later, this sort of a precision indicates a drift close to zero, which provides an initial constraint.

Most of ξ values are close to 1, indicating that a drift detection is currently possible when the drift rate is close to its maximum value. Moreover, ξ close to 2 or 3 are associated with a high eccentricity.

3.5. Comparison with previous results

As already mentioned, previous studies (Farnocchia et al. 2013b; Nugent et al. 2012b) also detected a drift in semi-major axis for several NEAs. We compare our values with those of previous studies (Fig. 2).

Fig. 2 Comparison of the drift in semi-major axis we obtained and that of previous studies (Farnocchia et al. 2013b; Nugent et al. 2012b).

Our results are in good agreement with the results previously published. The difference can be mainly explained by additional observations. For example, Nugent et al. (2012b) used observations until 2012, Farnocchia et al. (2013b) until 2013, and we used observations until May 2014. A different threshold for rejection may also explain some differences. For example, the drift for 1999MN is sensitive to the rejection threshold because by using a less restrictive threshold, we have ȧ = (43.2 ± 10.7) × 10^-4 au/My , which is closer to the drift obtained by Farnocchia et al. (2013b). This may indicate a problem with astrometry for this object that should be remeasured.

Finally, considering the 1σ uncertainty, the greater discrepancy with previous results is for (101955) Bennu, with a difference of about 4σ. An attentive control of the astrometry and radar measurements do not allow us to explain this discrepancy.

3.6. Discussion

The drift in semi-major axis can be related to physical parameters. For example, as the diurnal component is dominant, Vokrouhlický et al. (2008) provides a relation between the drift, the obliquity of spin axis γ, the bulk density ρ_b, the diameter D , and the diurnal thermal parameter Θ, i.e. $\begin{array}{ccc}\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}\mathrm{\propto }\frac{\cos \mathit{\gamma }}{{\mathit{\rho }}_{\mathrm{b}}\mathit{D}}\frac{\mathrm{\Theta }}{\mathrm{1}\mathrm{+}\mathrm{\Theta }\mathrm{+}\mathrm{0.5}{\mathrm{\Theta }}^{\mathrm{2}}}\mathrm{·}& & \end{array}$(4)On the one hand, the knowledge of physical parameters can provide expected drift rates. For example, Nugent et al. (2012a) used diameters and albedo from WISE data to predict drift rates due to the Yarkovsky effect for 540 NEAs. On the other hand, the knowledge of da/dt can constrain physical parameters D,ρ,orΘ. In particular, if one of these parameters is already known, it would be helpful for the determination of other unknown parameters.

We also note the relation between the diameter of the NEA and its drift value. Diameter, absolute magnitude, and albedo of NEAs are available on the EARN database⁶ and when they are not available, we consider the diameter (in km) given by the relation with the absolute magnitude H and the albedo p_v and a parameter K = 1329 km (Pravec & Harris 2007), i.e. $\begin{array}{ccc}\mathit{D}\sqrt{{\mathit{p}}_{\mathit{v}}}\mathrm{=}\mathit{K}\mathrm{\times }{\mathrm{10}}^{\mathrm{-}\mathit{H}\mathit{/}\mathrm{5}}\mathit{.}& & \end{array}$(5)As most of the NEAs we selected have a diameter smaller than 0.6 km, we can obtain statistics for the distribution of the diameter and the mean value of the drift (in absolute value). For NEAs bigger than 0.6 km, we do not have enough objects to obtain consistent statistics. Figure 3 represents this distribution for small NEAs (D ≤ 0.6 km). The number of components of each bar of the histogram is also indicated. Even if the drift in semi-major axis does not depend only on the diameter (see Eq. (4)), we note the expected 1/D dependence of | ȧ |.

Fig. 3 Distribution of diameter and mean value of ȧ for NEAs smaller than 0.6 km. Number of component is indicated above each bar of the histogram and the standard deviation of the mean is represented by error bar.

Another interesting result can be found in the mechanism of NEA transportation described in Bottke et al. (2002) because one of the consequences of this mechanism is an excess of retrograde rotators in the NEA population. According to La Spina et al. (2004), 67% of the NEAs should be retrograde.

Actually, in Eq. (4), the sign of ȧ allows us to distinguish retrograde NEAs (negative drift) and prograde NEAs (positive drift). For the 43 selected NEAs from Table 4, we found 77% of retrogade rotators, which is close to the distribution given in La Spina et al. (2004).

4. The precovery context

As a long period of observations is necessary to detect an accurate drift in semi-major axis, data mining and precovery observations can appear very important in that context. The case of (152563) 1992BF studied in Vokrouhlický et al. (2008) is a significative example because four precoveries were available and reduced with a modern catalogue to refine astrometry and then orbit determination.

Nowadays, photographic plates used until the 1990s can be digitised and reduced with modern techniques and modern stellar catalogues (Robert et al. 2011) as the current observations on CCD frames. In this context, old observations can have a precision equivalent to current observations.

This section aims to investigate the case of NEAs with existing precoveries for which a new reduction of these observations would help to better determine the drift in semi-major axis. Then, by studying two specific NEAs, we highlight favourable periods that allow a better uncertainty of the drift, showing that data mining has to be considered a great opportunity in the detection of drift.

4.1. Existing precoveries

In this section, we identify NEAs for which precovery observations have been isolated but not used for orbit determination because of their inaccuracy. These observations allow us to extend the orbital arc for several tens of years and, in the case of a new reduction, they would allow for a better determination of the drift in semi-major axis, or at least a smaller uncertainty. Depending on the quality of the precovery observations, we consider two different precisions for the precoveries with a new reduction: 300 mas⁷ corresponding to accurate observations, and 1 arcsec corresponding to low-quality observations. We determine σ_ȧ considering both cases. In any case, the new reduction will decrease the uncertainty in the drift, but to select NEAs for which the new reduction would be helpful for the drift detection, we filter NEAs with the ratio ȧ_max/σ_ȧ, where σ_ȧ is obtained considering a 300 mas precision for precovery. In fact, in the case in which the drift is close to its maximum ȧ_max, this ratio gives an approximate S/N obtained with the new reduction. Table 5 presents 14 NEAs with a ratio ȧ_max/σ_ȧ ≥ 1 and indicates the year of the precovery, the length of additional arc (the time between the first observation used in orbit determination and the precovery), the number of precoveries, the current uncertainty in the drift in semi-major axis and the expected uncertainty of the drift in the case the precoveries can be reduced with an precision of 1 arcsec and 300 mas and the ratio ȧ_max/σ_ȧ.

Currently, these NEAs have a bad determination of the drift, but a new reduction of the precoveries would help to decrease the uncertainty. Even if the precision of precoveries with the new reduction is 1 arcsec (corresponding to a low-quality observation), we can expect a large improvement of the drift uncertainty especially as the length of orbital arc is increasing. Obviously, the uncertainty of the drift is much smaller as we consider 300 mas for the precision of precoveries.

4.2. Importance of data mining

Photographic plates used before the 1990s were used for various purposes and surveys. The field of view covered is usually several degrees and many objects can appear on the plates, in particular, NEAs leading to precovery observations. In this section, we investigate the impact of the date and the precision of the precovery on the drift uncertainty. To measure the importance of the precovery observations, we consider the case of (152563) 1992BF , except with only the 1992–2011 observations (i.e. without the four precoveries of 1953). The idea is to understand why the 1953 precoveries are important and what happen if the precovery date changes.

In this context, we consider one additional precovery observation made between 1950 and 1990 and we determine the uncertainty of the detected drift with this precovery observation. For this, we consider various values of precision for the precovery observation: 3 arcsec, which corresponds to a bad astrometry (inaccurate stellar catalogue, bad quality of photographic plates), 300 mas, which corresponds to mean precision of current observations, and 30 mas , which could correspond to the precision obtained with the future Gaia stellar catalogue (see Sect. 5).

Figure 4 shows the uncertainty of the drift in semi-major axis obtained in that context for (152563) 1992BF. The grey zone under the curve indicates period for which σ_ȧ is smaller than 10^-4 au/My, allowing for an accurate detection of the drift.

If the precovery observation has an precision of 3 arcsec, then most of the time σ_ȧ is close to 4.4 × 10^-4 au/My , which is the uncertainty obtained with only 1992–2011 observations. During short period of time (1953, 1960), the uncertainty becomes close to 2 × 10^-4 au/My. These periods correspond to close approaches with Earth, when the geocentric distance of the asteroid is smaller than 0.4 au. With a 300 mas precision for the precovery observation, we still have these favourable periods (1953, 1960, 1966) where σ_ȧ becomes smaller than 10^-4 au/My, but it becomes even smaller when the precovery is farther in the past. This is even more the case with a precision of 30 mas for the precovery observation, where σ_ȧ ≤ 10^-4 au/My , whatever the date of precovery observation before 1965 might be.

Close approaches with Earth appear as good opportunities to improve the uncertainty of the drift in semi-major axis. Indeed, NEAs move faster in the celestial sphere leading to stronger constraints on their positions and in orbit determinations. Besides, during close approaches, NEAs appear brighter than usual and can leave traces on photographic plates.

Fig. 4 Uncertainty of the drift in semi-major axis related to the date of the precovery observations for (152563) 1992BF, with a precision of 3 arcsec, 300 mas, and 30 mas from top to bottom, for the precovery observation. Unit is 10-4 au/My for σȧ in ordinate.

We repeat the study for (99942) Apophis, which is a NEA with a very close approach to Earth on April 2029, at about 38 000 km. Currently, we have σ_ȧ = 8.778 × 10^-4 au/My and very small S/N, which is obviously not significant, but provides an initial constraint.

Figure 5 shows that the uncertainty of the drift could be reduced by a precovery observation performed during favourable periods (1965, 1972, 1980, 1990), which corresponds to close approach with Earth when a geocentric distance is smaller than 0.1 au. On these occasions, Apophis could have been imaged on photographic plates. Actually, we found a Beijing Observatory plate that covers the position of Apophis on January 1981. Unfortunately, because of its apparent magnitude of about 18, Apophis could not be detected on the plate. Nevertheless, in the case of a positive precovery observation at this date and with a precision of 300 mas, it would be possible to have a drift uncertainty of about 0.8 × 10^-4 au/My.

Fig. 5 Uncertainty of the drift in semi-major axis related to the date of the precovery observations for (99942) Apophis, with a precision of 3 arcsec, 300 mas, and 30 mas from top to bottom, for the precovery observation. Unit is 10-4 au/My for σȧ in ordinate.

Using astrometric observations, there are two ways to decrease the drift uncertainty: by extending the orbital arc of the NEA (see Fig. 1) and during NEA close approaches with Earth. Figures 4 and 5 show that the drift uncertainty σ_ȧ is smaller when the geocentric distance is small. During these periods, NEA appears brighter and during past close approaches, bright NEA could have been imaged on photographic plates. This is the case for several objects, including (152563) 1992BF, for which the photographic plates come from the Palomar Mountain series that was used for the DSS survey. This is an example of an exploited photographic archive, but other observatories in the world may have unexploited archives and could bring new observations to improve orbit determination and drift detection.

5. The Gaia mission context

The space mission Gaia was launched successfully on December 2013. The aim is to have a 3D-map of our Galaxy. In the case of solar system objects, the satellite will enable us to map thousands of main belt asteroids and NEAs down to magnitude 20. Gaia will observe several times each NEA for five years of the mission. Mouret & Mignard (2011) identified 64 promising NEAs for which Gaia observations could help the detection of a drift in semi-major axis with an uncertainty better than 5 × 10^-4 au/My.

Beyond the Gaia observations themselves, one of the best expected improvements will be the Gaia stellar catalogue allowing for orbital improvement and detection of drift in asteroid semi-major axis. The Gaia catalogue (Mignard et al. 2007) will provide unbiased positions of a billion stars down to magnitude 20 and with a precision depending on magnitude (7 μas at ≤10 mag; 12–25 μas at 15 mag and 100–300 μas at 20 mag).

With such a precision, new processes of reduction will be necessary, such as propagation models to third order for star proper motions, differential aberration, atmospheric absorption, colour of stars, etc. Finally, by including these corrections and considering the limitation due to the seeing, an uncertainty of 10–20 mas can be expected on the position of any asteroid reduced with the Gaia stellar catalogue.

Moreover, the Gaia stellar catalogue will enable us to reduce old photographic plates. The reduced position will be affected by other errors or limitations, including the distortion or instrument deformations; the granularity of the photographic plate; the digitisation; the exposure time; and the proper motions of stars. Depending on the magnitude and the epoch, the error would be, for example, about 12 mas for 20 mag star and 1.5 mas for a 16 mag star at epoch 1950 (Mignard 2013).

Another problem could be raised from a lack of meta-data in time observations. The time of the observations could be given with a low precision (to the minute, for example), or an error clock in the time of observation could possibly lead to inaccurate or even worse unusable observations.

In this study, we do not consider clock error, but only errors due to digitisation, photographic plates, or instruments. Assuming a good plate quality, fine grain emulsion, and optimal plate measuring, Zacharias (2013) estimated an error budget of about 10–20 mas for an astrometric position extracted from digitised plate. For these reasons, we assume that the astrometric precision of the new reduction with Gaia stellar catalogue will be 10 mas for observations since 1990 (for mainly CCD observations), 30 mas for 1950–1990 observations, 50 mas for 1920–1950 observations, and 100 mas for observations before 1920 (mainly photographic plates).

To measure the improvement of the drift detection with the Gaia stellar catalogue, we choose several assumptions as to new reduction methods that make sue of the Gaia catalogue. It seems obvious that all of the observations could not be reduced again, thus we assume that only a few would be reduced. As mentioned in Desmars et al. (2013), it appears that in such a case, the best solution for orbit determination is to reduce the first and last observations rather than a random part of them. Actually, the first and last observations, by virtually increasing the length of orbital arc, bring more important constraints on the orbital motion of the asteroid. In that context, we consider five cases (also considered in Desmars et al. 2013):

• Case 1: we examine current observations with their current uncertainty;

• Case 2: we assume that the first and the last observations can be reduced with the Gaia stellar catalogue;

• Case 3: we assume that the first five and last five observations can be reduced with the Gaia stellar catalogue;

• Case 4: we assume that the first ten and the last ten observations can be reduced with the Gaia stellar catalogue;

• Case 5: we assume that all observations can be reduced with the Gaia stellar catalogue.

Table 6 presents the distribution of drift uncertainty for the 1629 numbered NEAs according to various cases of new reduction with Gaia catalogue.

By considering different thresholds for the uncertainty, we note that the Gaia catalogue will be helpful to get accurate detection of drift in semi-major axis. Even in the minimal case where only two observations can be reduced with the Gaia catalogue (the first and last observations), 31 NEAs will have σ_ȧ ≤ 10^-4 au/My, whereas only eight NEAs currently reach this level of uncertainty. The more we reduce observations with Gaia catalogue, the more the quality of the drift detection will improve. In the extreme case where all observations could be reduced with Gaia catalogue, 688 (more than 40% of numbered NEAs) will have σ_ȧ ≤ 10^-4 au/My.

Considering the maximum value of the drift, the ratio ȧ_max/σ_ȧ also helps to measure the impact of the Gaia catalogue in the detection of drift. Table 7 presents the distribution of ratio for the different cases of reduction. As mentioned in previous section, the ratio can provide an approximation of the expected S/N in the case of the drift is close to its maximal value. So ȧ_max/σ_ȧ ≥ 3 could lead to reliable detection of the drift. Currently, only 15 objects reach this ratio, but with the minimal case with only two observations reduced with Gaia catalogue, 70 NEAs will have ȧ_max/σ_ȧ ≥ 3. Reducing 161 and 223 NEAs with 10 and 20 observations, respectively, using the Gaia catalogue, will fulfill the same condition. Finally, in the extreme case where all observations could be reduced again, 938 NEAs (more than 57%) will have a ratio greater than 3.

The Gaia catalogue will be very useful for orbit determination and for accurate detection of drift in semi-major axis, even with a few number of new reduction with Gaia catalogue. Of course, when the Gaia catalogue will be available in 2022 (or even with intermediate releases during the mission), new observations will be done and the drift uncertainty will be even better for most of NEAs.

6. Conclusion

The Yarkovsky effect produces a drift in the semi-major axis of NEAs. With a precise process of orbit determination, we have detected consistent drift for 46 NEAs. The drift value we measured are in agreement with values obtained by statistical methods. Moreover, the robustness has been tested for three representative NEAs. The drift in semi-major axis can be related to other physical parameters, such as the diameter, the obliquity, or the thermal parameter. The knowledge of the drift allows us to constrain these physical characteristics. For example, the prograde or retrograde rotation of the NEA can be determined by the sign of ȧ. The ratio of 78% of retrogade rotators is consistent with the expected ratio explaining the mechanism of NEA transportation.

We have also investigated the importance of precovery observations in the accurate detection of drift. In particular, we have shown that the drift uncertainty can be deeply improved by reducing only one precovery with modern techniques and modern stellar catalogue. The improvement is more important during close approach with Earth, which usually corresponds to a good period of visibility for the NEA (smaller apparent magnitude). During previous close approaches, some NEAs may have been imaged on photographic plates, and thus could provide accurate observations in order to refine orbit determination or, even more importantly, to help to detect a drift. We have already identified 14 NEAs with unused precoveries because of bad astrometry that can help to refine the detection of the drift.

Finally, we have highlighted and quantified the impact of the Gaia catalogue in the drift determination. With this catalogue, future and also past observations could be reduced with very high precision (a few tens of mas). If we reduce a few current observations with the Gaia catalogue, we can deeply increase the uncertainty in the drift detection. For example, with only two observations reduced again with the Gaia precision, 31 NEAs will have σ_ȧ ≤ 10^-4 au/My whereas only 8 NEAs currently reach this uncertainty. The Gaia catalogue opens up new perspectives in the drift detection.

The reduction of old photographic plates of solar system objects in the context of the Gaia mission is the topic of the NAROO⁸ project (Arlot 2013). This project investigates natural satellites and asteroids (trans-Neptunian objects, NEA, or Main Belt asteroids) in order to refine the dynamics of these objects. The detection of the drift in semi-major axis of NEAs is an important part of this project.

Acknowledgments

This work is supported by CNPq grant 161605/2012-5.

References

Appendix A: Maximum drift

To highlight unreliable detection, we introduce an estimation of a maximal drift in semi-major axis related to the absolute magnitude H. Vokrouhlický (1998) then Vokrouhlický et al. (2000) gave an relation between the drift in semi-major axis due to diurnal component of the Yarkovsky effect and several physical parameters, valid for circular orbit. We find $\begin{array}{ccc}{\left(\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}\right)}_{\mathit{d}}\mathrm{=}\mathrm{-}\frac{\mathrm{8}\mathit{\alpha }}{\mathrm{9}\mathit{n}}\mathrm{\Phi }\mathrm{\left(}\mathit{a}\mathrm{\right)}\frac{\mathit{G}\sin \mathit{\delta }}{\mathrm{1}\mathrm{+}\mathit{\lambda }}\cos \mathit{\gamma }& & \end{array}$(A.1)where α is the absortivity (complement of Bond albedo A), n is the mean motion, γ is the obliquity, and Φ(a) is the standard radiation function (inversely proportional to the diameter D, the bulk density ρ, and the semi-major axis a). The expression $\frac{\mathit{G}\sin \mathit{\delta }}{\mathrm{1}\mathrm{+}\mathit{\lambda }}$ is fully detailled in Vokrouhlický et al. (2000) and can be approximated by $\begin{array}{ccc}\frac{\mathit{G}\sin \mathit{\delta }}{\mathrm{1}\mathrm{+}\mathit{\lambda }}\mathrm{=}\frac{\mathrm{-}\mathrm{0.5}\mathrm{\Theta }}{\mathrm{1}\mathrm{+}\mathrm{\Theta }\mathrm{+}\mathrm{0.5}{\mathrm{\Theta }}^{\mathrm{2}}}& & \end{array}$(A.2)where Θ is the diurnal thermal parameter.

By using μ = n²a³ and numerical values for Φ function, we have $\begin{array}{ccc}{\left(\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}\right)}_{\mathit{d}}& \mathrm{=}& \mathrm{6.412}\mathrm{\times }{\mathrm{10}}^{\mathrm{3}}\frac{\mathit{\alpha }}{\mathit{D\rho }\sqrt{\mathit{a}}}\frac{\mathrm{0.5}\mathrm{\Theta }}{\mathrm{1}\mathrm{+}\mathrm{\Theta }\mathrm{+}\mathrm{0.5}{\mathrm{\Theta }}^{\mathrm{2}}}\cos \mathit{\gamma ,}\end{array}$(A.3)with the convenient units, a in au and $\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}$ in au/My. The expression $\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}$ is maximum when γ = 0 or 180 degrees and $\frac{\mathrm{0.5}\mathrm{\Theta }}{\mathrm{1}\mathrm{+}\mathrm{\Theta }\mathrm{+}\mathrm{0.5}{\mathrm{\Theta }}^{\mathrm{2}}}$ reaches its maximum value 0.207 for $\mathrm{\Theta }\mathrm{=}\sqrt{\mathrm{2}}$.

Finally, $\begin{array}{ccc}{\left|{\left(\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}\right)}_{\mathit{d}}\right|}_{\max }\mathrm{=}\mathrm{1.327}\mathrm{\times }{\mathrm{10}}^{\mathrm{3}}\frac{\mathit{\alpha }}{\mathit{D\rho }\sqrt{\mathit{a}}}\mathit{.}& & \end{array}$(A.4)To obtain a simple relation, we can consider a = 1 au, a mean density ρ = 2600 kg m^-3 et a mean albedo p_v = 0.154 for NEAs (Chesley et al. 2002). The Bond albedo will be A = 0.060 using the relation from Bowell et al. (1989)A = (0.290 + 0.684G)p_v (with G = 0.15) and α = 0.939.

The diameter can be approximated by the relation (Pravec & Harris 2007) $\begin{array}{ccc}\mathit{D}\sqrt{{\mathit{p}}_{\mathit{v}}}\mathrm{=}\mathrm{1329}\mathrm{km}\mathrm{\times }{\mathrm{10}}^{\mathrm{-}\mathit{H}\mathit{/}\mathrm{5}}& & \end{array}$with H the absolute magnitude.

Finally, the maximum drift can be approximated with the absolute magnitude H as $\begin{array}{ccc}{\left|{\left(\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}\right)}_{\mathit{d}}\right|}_{\max }& \mathrm{=}& \end{array}$(A.5)

All Tables

All Figures

	Fig. 1 Observed arc, uncertainty of the drift in semi-major axis (in au/My), and absolute magnitude H for the 1629 numbered NEAs. NEAs with radar measurements are represented with crosses and NEAs without radar measurements with bullets.
In the text

	Fig. 2 Comparison of the drift in semi-major axis we obtained and that of previous studies (Farnocchia et al. 2013b; Nugent et al. 2012b).
In the text

	Fig. 3 Distribution of diameter and mean value of ȧ for NEAs smaller than 0.6 km. Number of component is indicated above each bar of the histogram and the standard deviation of the mean is represented by error bar.
In the text

	Fig. 4 Uncertainty of the drift in semi-major axis related to the date of the precovery observations for (152563) 1992BF, with a precision of 3 arcsec, 300 mas, and 30 mas from top to bottom, for the precovery observation. Unit is 10-4 au/My for σȧ in ordinate.
In the text

	Fig. 5 Uncertainty of the drift in semi-major axis related to the date of the precovery observations for (99942) Apophis, with a precision of 3 arcsec, 300 mas, and 30 mas from top to bottom, for the precovery observation. Unit is 10-4 au/My for σȧ in ordinate.
In the text