© ESO 2018

1 Introduction

Several complex phenomena cause asteroid orbital evolution to be a difficult problem. By definition, near-Earth asteroids (NEAs) experience close approaches with terrestrial planets, which are the main source of chaos in their orbital evolution. Small perturbations, such as non-gravitational perturbations, can significantly affect a NEA trajectory because of this chaoticity.

The Yarkovsky effect is due to the recoil force undergone by a rotating body as a consequence of its anisotropic thermal emission (Vokrouhlický et al. 2000, 2015a). The main manifestation of the Yarkovsky effect is a secular semimajor axis drift da∕dt, which leads to a mean anomaly run-off that grows quadratically with time. Typical values of this perturbation for sub-kilometre NEAs are da∕dt ≃ 10⁻⁴–10⁻³ au My⁻¹. Because of its small size, the Yarkovsky effect can only be detected for asteroids with a well-constrained orbit.

The Yarkovsky effect is a non-gravitational perturbation that is generally split into a seasonal and a diurnal component. The seasonal component arises from the temperature variations that the heliocentric asteroid experiences as a consequence of its orbital motion. An explicit computation of the corresponding acceleration is not easy, even for a spherical body, and becomes very difficult for complex shaped bodies. On the other hand, the diurnal component is due to the lag between the absorption of the radiation coming from the Sun and the corresponding re-emission in the thermal wavelength. The surface of a rotating body illuminated by the Sun is warmed by solar radiation during the day and cools at night. It is worth noticing that the diurnal and seasonal components have different consequences on the secular semimajor axis drift. In particular, the diurnal effect produces a positive drift for prograde rotators and a negative drift for retrograde rotators, whereas the seasonal secular drift is always negative. The magnitude of the diurnal effect is generally larger than that of the seasonal effect (Vokrouhlický et al. 2000).

The Yarkovsky effect is key to understanding several aspects of asteroids dynamics:

1. Non-gravitational forces can be relevant for a reliable impact risk assessment when in the presence of deep planetary encounters or when having a long time horizon for the potential impact search (Farnocchia et al. 2015b). As a matter of fact, both of these factors call for a greater consideration of the sources of orbit propagation uncertainty, such as non-gravitational perturbations. Currently, there are four known cases that required the inclusion of the Yarkovsky effect in term of hazard assessment: (101955) Bennu (Milani et al. 2009; Chesley et al. 2014), (99942) Apophis (Chesley 2006; Giorgini et al. 2008; Farnocchia et al. 2013a; Vokrouhlický et al. 2015b), (29075) 1950 DA (Giorgini et al. 2002; Farnocchia & Chesley 2014), and (410777) 2009 FD (Spoto et al. 2014).

2. The semimajor axis drift produced by the Yarkovsky effect has sculpted the main belt for millions of years (Vokrouhlický et al. 2006). The Yarkovsky effect is crucial to understanding the ageing process of asteroid families and the transport mechanism from the main belt to the inner solar system (Vokrouhlický et al. 2000). The Yarkovsky effect has still not been measured in the main belt, thus Spoto et al. (2015) used a calibration based on asteroid (101955) Bennu to compute the ages of more than 50 families in the main belt. The improvement in the detection of the Yarkovsky drift for NEAs represents a new step forward in creating a reliable chronology of the asteroid belt.

3. Yarkovsky detections provide constraints on asteroid physical properties. Two remarkable efforts in estimating the bulk density from the Yarkovsky drift have already been carried out for two potentially hazardous asteroids: (101955) Bennu (Chesley et al. 2014) and (29075) 1950 DA (Rozitis et al. 2014). Furthermore, the dependence of the Yarkosvky effect on the obliquity can be useful to model the spin axis obliquity distribution of NEAs (Tardioli et al. 2017).

There have been several efforts to model and determine the Yarkovsky effect on the NEA population. The first detection of the Yarkovsky effect was predicted for the asteroid (6489) Golevka in Vokrouhlický et al. (2000) and achieved in 2003 thanks to radar ranging of this object (Chesley et al. 2003). Later, the Yarkovsky effect played a fundamental role in the attribution of four 1953 precovery observations to the asteroid (152563) 1992 BF (Vokrouhlický et al. 2008). Moreover, Chesley et al. (2014) detected the Yarkovsky effect acting on (101955) Bennu from astrometric observations and high-quality radar measurements over three apparitions. Currently, asteroid Bennu has the best determined value for the Yarkovsky acceleration, which also led to an estimation of its bulk density (Chesley et al. 2014). In general, Nugent et al. (2012) provided a list of 13 Yarkovsky detections with S∕N > 3 and later work increased this number to 21 (Farnocchia et al. 2013b). The most recent census is from Chesley et al. (2016), which identified 42NEAs with valid Yarkovsky detections. Both Farnocchia et al. (2013b) and Chesley et al. (2016) flag spurious cases based on whether the detected drift is compatible with the physical properties of the corresponding object and the Yarkovsky mechanism. Since the number of significant Yarkovsky detections in the NEA population is steadily increasing, we decided to update the list.

2 Method

2.1 Force model

Usually, there is not enough information on the physical model of an asteroid to directly compute the Yarkovsky acceleration through a thermophysical model. Instead, evidence of the Yarkovsky-related drift may be detectable from the observational dataset via orbit determination. Indeed, a gravity-only model may not provide a satisfactory fit to the available data. A Yarkovsky detection is more likely when a very accurate astrometric dataset is available, especially in case of radar measurements at multiple apparitions or when the observational arc is long, thus allowing the orbital drift to become detectable. In such cases, a force model that also includes the Yarkovsky acceleration could result in a better fit to the observations.

We modelled the Yarkovsky perturbation with a formulation that depends on a single dynamical parameter, to be determined in the orbital fit together with the orbital elements. Since the secular perturbation caused by the Yarkovsky effect is a semimajor axis drift, we used a transverse acceleration $$a_t=A_{2}g(r)\widehat{t},$$(1)

as in Marsden et al. (1973) and Farnocchia et al. (2013b). In Eq. (1), A₂ is a free parameter and g(r) is a suitable function of the heliocentric distance of the asteroid. In particular we assumed a power law $$g(r)={\left(\frac{r_0}{r}\right)}^d,$$

where r₀ = 1 au is used as normalization factor. The exponent d depends on the asteroid and is related to the thermophysical properties of the asteroid. Farnocchia et al. (2013b) showed that the value of d is always between 0.5 and 3.5. They used d =2 for all asteroids because no thermophysical data are available. In our analysis we adopted the same values for d, apart from (101955) Bennu for which the value d = 2.25 is assumed (Chesley et al. 2014).

Typical values of the Yarkovsky acceleration for a sub-kilometre NEA are 10⁻¹⁵–10⁻¹³ au d⁻². As a consequence, to reliably estimate the Yarkovsky effect, the right-hand side of the equations of motion has to include all the accelerations down to the same order of magnitude. Our force model includes the gravitational accelerations of the Sun, the eight planets, and the Moon based on the Jet Propulsion Lab (JPL) planetary ephemerides DE431 (Folkner et al. 2014). To ensure a more complete force model, we also included the contributions coming from 16 massive main belt bodies and Pluto. All the gravitational masses we used are listed in Table 1. Since we compared our results with those obtained by JPL, we point out that the JPL team used the 16 most massive main belt asteroids as estimated by Folkner et al. (2014), which produces a slight difference, both in the list and masses.

The relativity model includes the relativistic contribution of the Sun, the planets and the Moon. In particular, we used the Einstein–Infeld–Hoffman equations, namely, the equations of the approximate dynamics ofa system of point-like masses due to their mutual gravitational interactions, in a first order post-Newtonian expansion, as described in Einstein et al. (1938), Will (1993), and Moyer (2003).

2.2 Statistical treatment of the astrometry

The statistical treatment of the astrometry is key to a reliable orbit determination. The differential corrections procedure provides the nominal orbit of the asteroid and its uncertainty (Milani & Gronchi 2010, Chap. 5), which strongly depend upon the observationsaccuracy and error modelling. For the computations carried out for this paper, we used the debiasing and weighting scheme provided in Farnocchia et al. (2015a). The JPL team uses the more recent Vereš et al. (2017) weighting scheme.

The astrometric data usually can contain outliers that can affect the solution of the orbit determination. To remove erroneous observations from the fit, we applied the outlier rejection procedure described in Carpino et al. (2003).

In addition to our default data treatment, we applied ad hoc modifications for the following cases:

(152563) 1992 BF. The four 1953 precovery observations of this NEA have been carefully re-measured in Vokrouhlický et al. (2008). We adopt the given positions and standard deviations, the latter being 0.5 arcsec in right ascension and 1 arcsec in declination.

2009 BD. This object is one of the smallest NEAs currently known (Mommert et al. 2014b) and thus solar radiation pressure (SRP) affects its orbit. A direct detection of the area-to-mass ratio is contained in Micheli et al. (2012), which provided high-qualityastrometry from Mauna Kea and replaced all the observations from the Tzec Maun Observatory (H10) with a single position. For these observations we set data weights based on the uncertainties provided by Micheli et al. (2012) and we included both the Yarkovsky effect and SRP in the orbital fit (see Sect. 6).

2011 MD. As 2009 BD, this asteroid is very small and is among those for which we determined both the Yarkovsky effect and SRP. 2011 MD has been observed during the 2011 very close approach with the Earth. Despite the short arc of three months, a very large number of optical observations of 2011 MD were collected, precisely 1555. Following the strategy presented in Mommert et al. (2014a), we relaxed the weights for the observations collected during the close approach¹ and we added the Spitzer detection (on February 11, 2014), which extends the observation arc by almost three years.

2015 TC₂₅. Asteroid 2015 TC₂₅ was discovered by the Catalina Sky Survey in October 2015, just two days before an Earth flyby at 0.3 lunar distances. It is one of the smallest asteroids ever discovered, about 2 m in diameter (Reddy et al. 2016), and the 2017 astrometry permits us to achieve an estimate of SRP. We are aware that for 2015 TC₂₅ the JPL team carried out a specific study (Farnocchia et al. 2017), which adopted ad hoc weights based on observer-provided uncertainty estimates. To handle this case, we used the same data treatment as JPL.

We note that it is desirable to keep the number of “manual” interventions on the observational data as small as possible. Indeed we are trying to figure out how to automatize the determination of the set of NEAs with significant and reliable Yarkovsky effect. Anyway, in some cases a manual intervention is needed to properly handle observational issues, for example too many close observations taken during a very close approach and affected by timing errors or remeasurement of old observations.

2.3 Starting sample of NEAs

As first sample of asteroids, we started selecting those objects in NEODyS² that have a formal uncertainty on the semimajor axis σ(a) < 3 × 10⁻⁹ au. The choice of the threshold for σ(a) comes from an order of magnitude estimate: for an asteroid with diameter 1 km the Yarkovsky drift is about 3 × 10⁻¹⁰ au y⁻¹, thus it causes a variation of 3 × 10⁻⁹ au in ten years. The value of σ(a) has to be computed at the mean epoch of the observations, since it is the best choice for the orbital fit quality. Moreover, this uncertainty threshold corresponds to a gravity-only fit. After the Yarkovsky coefficient is estimated, the uncertainty of the semimajor axis sharply increases because of the strong correlation between A₂ and the semimajor axis.

The list of asteroids satisfying this criterion contained 519 objects (as of February 2018). As a second step, we extracted from the JPL database a set of 89 asteroids having A₂ determined³. Among these, only 16 were not contained in our first list, thus we added them. Furthermore, we considered all the reliable detections from Farnocchia et al. (2013b) and it turned out that only 3 asteroids did not belong to any of the previous lists, thus we added these to our sample as well.

Summarizing, we started with a sample of 519 + 16 + 3 = 538 objects. For each one of these objects we performed an orbital fit for the initial conditions together with the Yarkovsky parameter A₂, without any a prioriconstraint. For a few of these objects we also estimated SRP. As a result of the fit, we derived the signal-to-noise ratio ${\text{S/N}}_{A_2}$ of the A₂ parameter, obtaining 101 detections with $S/N_{A_2}\ge 3$ and 437 with $S/N_{A_2}<3$, most of which showing a negligible S/N.

2.4 Yarkovsky expected value

By means oforbit determination, we determined a transverse acceleration directly from the astrometry. However, to claim that the measured acceleration is caused by the Yarkovsky effect we needed to make sure that its magnitude is compatible with the physical properties of the object and the Yarkovsky mechanism. Therefore, we provide an expected value of the Yarkovsky-related orbital drift.

In Farnocchia et al. (2013b), an expected value for A₂ is computed by exploiting the diameter of the asteroid and scaling from the corresponding value of (101955) Bennu, the best determined and reliable Yarkovsky detection. In this paper, we make use of the Yarkovsky calibration as in Spoto et al. (2015) as follows: $${\left(\frac{\text{d}a}{\text{d}t}\right)}_{\text{exp}}={\left(\frac{\text{d}a}{\text{d}t}\right)}_{\mathcal{B}}\times \frac{\sqrt{a_{\mathcal{B}}}\left(1-e_{\mathcal{B}}^2\right)}{\sqrt{a}\left(1-e^2\right)}\frac{D_{\mathcal{B}}}{D}\frac{{\rho }_{\mathcal{B}}}{\rho }\frac{\cos \varphi }{\cos {\varphi }_{\mathcal{B}}}\frac{1-A}{1-A_{\mathcal{B}}},$$(2)

where D is the diameter of the asteroid, ρ is the density, ϕ is the obliquity (angle between the spin axis and the normal to the orbit plane), and A is the Bond albedo. We computed the latter from the geometric albedo p_v, using $A=\frac{1}{3}{p}_v$ (Muinonen et al. 2010). The symbols with a “$\mathcal{B}$” refer to asteroid (101955) Bennu, and the assumed values are listed in Table 2 with their references.

For the diameter D we used the known physical value when available. When the shape of the asteroid is not simple enough to be approximated by an ellipsoidal model, we used the dynamically equivalent equal volume ellipsoid dimensions to compute the equivalent diameter. In particular, this effort has been carried out for three asteroids, namely (4179) Toutatis (Hudson et al. 2003), (162421) 2000 ET₇₀ (Naidu et al. 2013), and (275677) 2000 RS₁₁ (Brauer et al. 2015). Otherwise, when no physical information are available, we estimated the diameter from the absolute magnitude H following the relation (Pravec & Harris 2007) $$D=1329\text{ km}\times {10}^{-H/5}\times \frac{1}{\sqrt{p_v}},$$

where the geometric albedo p_v is assumed to be p_v = 0.154 if unknown.

As shown in Eq. (2), the density is required to estimate the strength of the Yarkovsky effect for asteroids with small diameters. Carry (2012) reported a large number of asteroid densities that we used as starting point. However, in general, density estimates are more reliable and accurate for massive bodies and there is a trend of a decreasing density with diameter because of the increasing macroporosity⁴ resulting from the cascade of collisions suffered by the body (Carry 2012; Scheeres et al. 2015). We thus extrapolated the density of small asteroids from the density of large asteroids belonging to the same taxonomic class by increasing their macroporosity to that of Bennu (${\mathcal{P}}_{\mathcal{B}}=(40\pm 10)\%$, from Chesley et al. 2014). Such macroporosity is typical for (sub-)kilometre-sized asteroids, as illustrated by (25143) Itokawa, visited by the JAXA Hayabusa mission (Fujiwara et al. 2006). This is a modified version of the approach given in Spoto et al. (2015), still using Bennu for the scaling since it has the best estimated Yarkovsky acceleration and a comprehensive physical characterization⁵. Thus, the scaled density is given by $${\rho }_s=(1-{\mathcal{P}}_{\mathcal{B}})\rho ,$$(3)

where the density scaling factor is $1-{\mathcal{P}}_{\mathcal{B}}=0.60$ and ρ is the known density of the large asteroid. Equation (3) follows from the above assumptions and from the definition of macroporosity. We selected the large asteroids (4) Vesta, (10) Hygiea, (15) Eunomia, and (216) Kleopatra as representative of the taxonomic classes V, C, S, Xe, respectively. The density of the representative asteroids and their scaled values are listed in Table 3.

We used three sources of asteroid physical information: the database of physical properties of NEAs provided by E.A.R.N., the JPL Small-Body Database, and the data provided by the WISE mission, such as diameters and albedos (Mainzer et al. 2011). It is important to point out that we have no physical information concerning the large majority of the objects discussed in this paper. For instance, for the 44% of our detections with $S/N_{A_2}>2.5$ we have no physical data, for 62% we have no measured albedo values, and less than half of our detections can be assigned to a taxonomic class.

2.5 Filtering criterion

We used theYarkovsky-related expected value as a filtering criterion to understand whether the estimated orbital drift d a∕d t is physically consistent with the Yarkovsky effect. If the estimated da∕dt is significantly larger than the maximum absolute expected value (assuming cos ϕ= ±1), the result is inconsistent with the Yarkovsky mechanism. We computed the indicator parameter $$\mathcal{S}=\left|\frac{\text{d}a/\text{d}t}{{(\text{d}a/\text{d}t)}_{\text{exp}}}\right|.$$(4)

Since most times there is very little to no physical information, we needed some margin on the upper threshold for $\mathcal{S}$, which therefore should be larger than 1. We filtered out the candidate detections with $\mathcal{S}>2$. The current maximum value allowed for $\mathcal{S}$ is empirical, but it could be refined. In particular, this upper threshold can be lowered when better data are available. Improving the computation of the expected value – thus decreasing the uncertainty of the indicator parameter $\mathcal{S}$ – requires at least a reliable taxonomic type (for the scaling needed for the density) and better diameters. Values of $\mathcal{S}$ greater than the maximum threshold indicate questionable results. These spurious detections should be investigated to find possible causes and solutions. In general, either the $\mathcal{S}$ value is too high to be compatible with an acceptable detection or it is barely above the maximum threshold, in such a way that additional information would clarify the situation and allow us to decide whether the detection is accepted or refused. For further details, see Sect. 4.

We point out that values $\mathcal{S}\ll 1$ are permitted. Indeed, Eq. (2) employed asteroid size and bulk density, thus $\mathcal{S}\ll 1$ means that the orbital drift is significantly lower than the maximum expected value. Several phenomena can lower the Yarkovsky effect: obliquity ϕ ≃ 90°, very large or very small thermal inertia, larger density than expected, or small rotation angular velocity. For instance, asteroid (85774) 1998 UT₁₈ has a rotation period of about 34 h, and indeed the indicator $\mathcal{S}$ is low (≃0.3, cf. Table A.1). Detections of this kind are significant detections of a weak Yarkovsky drift. A second class of weak Yarkovsky drifts can be defined. These are non-detections, that is $S/N_{A_2}<3$, but the asteroid has physical properties that would permit a significant detection if the Yarkovsky effect were maximized. Chesley et al. (2016) referred to these detections as weak detections. Despite the low S/N, the result of the A₂ estimation can provide useful constraints on the physical properties of the asteroid.

By combining the value of the $S/N_{A_2}$ coming from the orbital fit with the value of the filtering parameter $\mathcal{S}$, we divided our detections in three categories as follows:

• the detections satisfying both $S/N_{A_2}\ge 3$ and $\mathcal{S}\le 2$ are considered accepted;

• the category called marginal significance includes the asteroids for which $2.5<S/N_{A_2}<3$ and $\mathcal{S}\le 2$, plus (410777) 2009 FD and (99942) Apophis, both remarkable for their impact monitoring;

• the detections with $S/N_{A_2}>3$ and $\mathcal{S}>2$ are rejected because they have a value that is too high for the indicator parameter $\mathcal{S}$, suggesting that the detected A₂ signal isunrealistic or not explicable with the Yarkovsky effect (see 2003 RM in Chesley et al. (2016), or (4015) Wilson-Harrington in Sect. 4).

The results only include detections, that is we do not list the asteroids for which we did not find a significant Yarkovsky signal from the observational dataset (85% of the initial sample). Figure 1 provides an overall view of our classification. In particular, we consider the plane $(S/N_{A_2},\mathcal{S})$ and we indicate the detections of each class (but the rejected) with a different colour as follows: the accepted detections are indicated with a green dot; the marginal significance detections are represented with a blue dot, except for (410777) 2009 FD and (99942) Apophis, which are indicated with a blue asterisk (special cases); the rejected detections are indicated with a red cross.

Fig. 1 Graphical representation of the partition of the detections set into classes on the plane $(S/N_{A_2},\mathcal{S})$. We plot theaccepted detections (green dots), the marginally significance detections (blue dots) plus the special cases (410777) 2009 FD and (99942) Apophis (blue asterisk), and the rejected detections (red crosses).

3 Accepted and significant results

As explained in Sect. 2.3, we started from the initial sample of NEAs and we performed an orbital fit including the Yarkovsky parameter A₂. Then we applied the filtering procedure and obtained 87 detections, which are listed in Table A.1 and Table A.2. For each asteroid, we report the value of the absolute magnitude H, the A₂ parameter along with its uncertainty, the S/N of A₂, the value of the semimajor axis drift da∕dt, the indicator parameter $\mathcal{S}$, and the available physical data such as the geometric albedo p_v, the diameter D, density ρ, and taxonomic class. As explained in Sect. 2.4, when no information on the diameter is available, we infer a diameter from the absolute magnitude. We denote these cases with an asterisk (*) in the diameter column. When the albedo is not directly measured but the taxonomic class is known, we assign the albedo according to Binzel et al. (2002) and denote the albedo value with a dagger (†). When the albedo is not known we assume p_v = 0.154, and indicate thisvalue with a “d”. We note that 0.154 is a mean value, which has a low probability of being accurate because of the bimodality of the albedo distribution of NEAs. Most asteroids are either significantly brighter or significantly darker. Thus, when a diameter D is derived from the absolute magnitude and this default albedo, its relative uncertainty is larger, and in turn the value of $\mathcal{S}$ (containing the factor 1∕D) is uncertain.

Some of the asteroids included in Tables A.1 and A.2 deserve dedicated comments.

(1566) Icarus. It is known that the 1968 observations of (1566) Icarus are affected by large timing errors. A possible solution to this problem is to include timing errors in the observations uncertainty possibly even removing the systematic timing errors. A possible alternative is to properly treat the correlation between the right ascension and declination. This operation will be made easier after the adoption of the new Astrometric Data Exchange Standard (ADES⁶). For now, it is possible to adapt the weighting scheme to underweight the observations during the 1968 close approach. The JPL team took this approach, but the Pisa team did not. By comparing the result of the two groups and also that of Greenberg et al. (2017), we can claim that the detection of the Yarkovsky effect is confirmed, even if there is a significant difference between the standard deviations (see Sect. 7), which is explained by the different weighting scheme.

(3908) Nyx. Asteroid (3908) Nyx is classified as V type, but it has many properties that are inconsistent with (4) Vesta. Thus, the density scaling is not performed using the density of Vesta as for the other V-type asteroids. Asteroid (5381) Sekhmet is a V-type with a diameter that is comparable to that of Nyx. The density of Sekhmet is (1.30 ± 0.65) g cm⁻³ (Carry 2012), compatible with the estimate in Farnocchia et al. (2014) and we also assume this value for Nyx.

(4179) Toutatis This asteroid shows a significant Yarkovsky detection with $S/N_{A_2}>3$ and $\mathcal{S}<2$ and is therefore included in the list of accepted detections. However, we point out that this detection is subject to substantial uncertainties beyond the formal uncertainties resulting from the fit to the optical and radar astrometry. The large aphelion of the orbit of Toutatis and the small magnitude of the Yarkovsky perturbation make the A₂ estimate sensitive to the set of main belt perturbers included in the force model and the uncertainty in their masses.

4 Rejected results

In this section we consider the significant detections that we rated as spurious, i.e. for which we obtained a Yarkovsky detection greater than reasonably expected from the Yarkovsky effect. Despite the S/N being less than 3, we also add(4015) Wilson-Harrington to this category as a special case, as explained below.

There are several reasons for refusing a detection. First, dynamical model problems can occur in a few cases, such as (4015) Wilson-Harrington. Second, the results are sometimes strongly dependent on few observations, typically old isolated observationsthat are separated by a long time interval from the bulk of the dataset. In these cases we usually reject the detection unless the precovery has been carefully remeasured, as for (152563) 1992 BF. Third, there are solutions with Yarkovsky that are affected by observational data of questionable reliability. The list of all the rejected detections contains 14 cases; see Table A.3. Below we provide dedicated comments for each rejected detection.

(4015) Wilson-Harrington. This asteroid was initially discovered in 1949 as a comet at the Palomar Sky Survey. It was named 107P/Wilson-Harrington, but then it was lost. Thirty years later the asteroid 1979 VA was discovered and, after the 1988 apparition, it was numbered as (4015) 1979 VA. On August 13, 1992, the IAU circular 5585 (Bowell et al. 1992) reported that the asteroid (4015) 1979 VA and the comet 107P/Wilson-Harrington were indeed the same object. Furthermore, no cometary activity was noted during the well-observed 1979–80 apparition, confirming that it is actually an extinct comet. The detection of the Yarkovsky parameter is indeed significant, but the $\mathcal{S}$ value indicates a value of the non-gravitational acceleration that is larger than that expected from the Yarkovsky effect. Since the observed arc contains the time span of cometary activity, the most likely interpretation is that the large detected transverse acceleration is caused by the out-gassing rather than the Yarkovsky effect. Furthermore, a dynamical model assuming a constant value for A₂, as this is the model we employed, is not sufficiently representative of the real orbital dynamics, given that the cometary activity has ceased. Thus in this case we consider that a non-gravitational effect has been detected, but not Yarkovsky.

In some cases, spurious detections are due to poor optical astrometry, often affecting isolated old observations. In this case we rejected the detection and a remeasurement of these old observations would be desirable to clarify the situation: the Yarkovsky signal could significantly either increase or disappear. This is the case for the asteroids listed below.

(260141) 2004 QT₂₄. The detected signal is strongly dependent on four observations in 1993 and 1998 from Siding Spring Observatory DSS.

(350751) 2002 AW. This asteroid has two isolated observations in 1991 from Palomar Mountain-DSS.

(39565) 1992 SL. This asteroid has one isolated observation in 1950 from Palomar Mountain.

(4486) Mithra. This asteroid has a signal that is strongly dependent on a single isolated observation in 1974 from Crimea-Nauchnij (MPC code 095).

(474158) 1999 FA. This object has one isolated observation in 1978, from Siding Spring Observatory. In agreement with Farnocchia et al. (2013b), we consider that the 1978 observation would need to be remeasured before accepting the Yarkovsky detection for (474158) 1999 FA.

(162421) 2000 ET₇₀. This asteroid has two isolated observations in 1977 from European Southern Observatory, La Silla DSS.

There are detections that have to be rated as spurious because $\mathcal{S}$ indicates a Yarkovsky drift that is way larger than expected even though the S/N of the A₂ parameter is greater than 3. This holds for 2010 KP₁₀, (308635) 2005 YU₅₅, (139359) 2001 ME₁, (142561) 2002 TX₆₈, and (192563) 1998 WZ₆. To confirm the reliability of our filtering criterion, we carefully checked each of these spurious detections. They show problematic astrometry, which resulted in a incorrect determination of the Yarkovsky effect.

A separate comment holds for (175706) 1996 FG₃, since it is a binary asteroid (Scheirich et al. 2015). The signal found for the Yarkovsky detection is likely due to the astrometric data treatment, as confirmed by the fact that it disappears when the weighting scheme proposed in Vereš et al. (2017) is applied. Once we have a Yarkovsky detection for this object it will be possible to compare with the Yarkovsky theory for binary asteroids as described by Vokrouhlický et al. (2005).

A remarkable case is (433) Eros. This asteroid shows a significant value for the Yarkovsky effect, but a very high value for the indicator parameter. Moreover, the obliquity of Eros is known to be ≃ 89° (Yeomans et al. 2000), therefore we would expect a value for $\mathcal{S}$ much less than 1. Thus, this detection is spurious, likely caused by historical data dating back to 1893 for which it is challenging to come up with a reliable statistical treatment.

5 Marginal significance

We now consider the marginal significance class, containing the detections for which $2.5<S/N_{A_2}<3$ and $\mathcal{S}\le 2$. These detections are physically meaningful since they satisfy the filtering on $\mathcal{S}$, but the signal to noise for the A₂ parameter as determined from the observations is not enough for a reliable detection. In addition, as mentioned before, we include two special casesin this category, namely (410777) 2009 FD and (99942) Apophis. These two objects show acceptable values of the indicator $\mathcal{S}$ but the S/N of the A₂ parameter isvery low (cf. Table A.4). Nevertheless, we decided to keep these objects because the Yarkovsky drift plays a fundamental role for its impact predictions (see the introduction). In this way we grouped 24 detections in this class, which are listed in Table A.4.

(99942) Apophis. Similar to Toutatis, Apophis also has a complex rotation, as shown in Pravec et al. (2014); Vokrouhlický et al. (2015b). However, the Yarkovsky effect is not significantly weakened by the tumbling state. Vokrouhlický et al. (2015b) used the available rotation state, shape, size, and thermophysical model of Apophis to predict the Yarkovsky semimajor axis drift. The drift obtained by fitting the astrometric data is compatible with the model prediction. We obtained d a∕dt = (−24.50 ± 13.58) × 10⁻⁴ au Myr⁻¹ for the fitted value, which is completely consistent with Vokrouhlický et al. (2015b). There is no question that the Yarkovsky effect has to be taken into account for Apophis to predict future motion, especially for impact hazard assessment (Chesley 2006; Giorgini et al. 2008; Farnocchia et al. 2013a).

(410777) 2009 FD. The Yarkovsky effect found is below the significance level, and nevertheless it has to be taken into account for long-term impact monitoring purposes (Spoto et al. 2014).

Maintaining a list of marginal significance detections is useful because they are candidates for future detections as observational data improves and increases.

6 Direct radiation pressure detection

Solar radiation pressure is a more complicated perturbation to detect. So far, SRP has only been detected for very small objects (H > 27) that experienced Earth encounters. Thus we started selecting the smallest asteroids of the initial sample (more precisely, those with H > 24) since the effect becomes larger for smaller size objects and we tried to detect SRP along with the Yarkovsky parameter.

The acceleration caused by SRP is radial and can be modelled with a single parameter A₁, $$a_r=A_{1}g(r)\widehat{r}.$$

In this equation, A₁ is a free parameter, and g(r) = 1∕r², where r is the heliocentric distance in astronomical units. Physically, the value of A₁ depends mostly on the area-to-mass ratio $\mathcal{A}/M$. The relationbetween these parameters is the following: $$A_1=\frac{{\Phi }_{\odot }}{c}\times C_R\times \mathcal{A}/M,$$

where c is the speed of light, Φ_⊙ is the solar radiation energy flux at 1 au whose value is Φ_⊙ ≃ 1.361 kW m⁻², and C_R is a coefficient (on the order of 1) depending upon the shape and optical properties of the surface.

The starting sample of asteroids for which we attempted an 8-dimensional fit contained ten objects. We found four accepted detections, i.e. $S/N_{A_1}\ge 3$, which are listed in Table A.5. For three asteroids of this category, namely, 2011 MD, 2012 LA, and 2015 TC₂₅, the Yarkovsky detection is not significant and thus the $\mathcal{S}$ value, although above the threshold in one case, does not provide any information. Concerning the area-to-mass ratio, we would compare the value of $\mathcal{A}/M$ with an expected value, as we do with the secular semimajor axis drift da∕dt, but this is not possible for now since the diameter is very uncertain and the other physical properties are currently unknown. Other fitted values of the area-to-mass ratio have already been determined for 2009 BD (Micheli et al. 2012), 2012 LA (Micheli et al. 2013), and 2011 MD (Micheli et al. 2014), although without including the Yarkovsky effect in the dynamical model. Asteroid 2006 RH₁₂₀ is listed separately from the others since we consider it spurious, as we explain in what follows.

2006 RH₁₂₀. This strange detection has already been discussed in Chesley et al. (2016). Our results are very compatible with those of that paper and we agree with the motivations provided to reject this detection. The most likely explanation for the high transverse acceleration can be the presence of some non-conservative force, for example mass shedding, outgassing, or micrometeorite flux, which can become as relevant as the Yarkovsky effect for objects of this size. The area-to-mass ratio, which results in a significant detection, is not compatible with the hypothesis that 2006 RH₁₂₀ is an artificial object.

7 Comparison with JPL results

As we already mentioned, the JPL database is regularly updated with the asteroids for which the orbital fit shows evidence of the Yarkovsky effect. The SRP is also regularly updated as well when appropriate. The results produced by two independent software areexpected to be different, but compatible. In order to compare them we compute the relative errors $${\epsilon }_r(A_2)=\frac{|A_2-A_2^{\text{JPL}}|}{{\sigma }_{A_2}}\text{and}{\epsilon }_r^{\text{JPL}}(A_2)=\frac{|A_2-A_2^{\text{JPL}}|}{{\sigma }_{A_2,\text{JPL}}},$$

where the superscript “JPL” refers to the JPL solution. To quantify the difference between the results presented in this paper and the JPL values, we use the quantity $${\chi }_{A_2}=\frac{|A_2-A_2^{\text{JPL}}|}{\sqrt{{\sigma }_{A_2}^2+{\sigma }_{A_2,\text{JPL}}^2}},$$

from (Milani & Gronchi 2010, Sect. 7.2). We consider compatible two solutions for which ${\chi }_{A_2}\le 1$.

Starting from the list of our accepted detections, we compared the results every time an asteroid is included in the JPL database of Yarkovsky effect detections. The results of the comparison are contained in Tables A.6 and A.7. Just for five asteroids in this list both the relative errors are greater than 1, even though never above 2.5. Using the metric given by ${\chi }_{A_2}$, we identify just 3 asteroids (marked with a star in Tables A.6 and A.7) whose detections are not fully compatible with respect to our criteria. Anyway, a ${\chi }_{A_2}$ moderately above 1 for 3 cases out of 92 being compared shows a strong agreement between our results and the results from JPL.

We note that this result is not a null test, that is the expected value of the difference in the estimated values of A₂ is not zero. This because the two computations have used two different astrometric error models, that is Farnocchia et al. (2015a) at NEODyS and Vereš et al. (2017) at JPL. The comparative results described in the last three columns indicate an exceptionally good agreement. This agreement may be interpreted as a validation of the procedures used both at NEODyS and at JPL, both to compute the Yarkovsky effect constants and to select the cases in which the results are reliable.

The comparison was also performed for the shorter list of objects for which we have both A₂ and A₁, that is both Yarkovsky effect and direct radiation pressure were included in the dynamical model. Table A.8 contains theS/Ns for both parameters in both solutions and all the metrics for the discrepancies. Apart from the results for 2006 RH₁₂₀, which are rated as spurious, the accepted results are fully consistent.

8 Impactmonitoring with non-gravitational parameters

A force model including non-gravitational forces is sometimes needed to make reliable impact predictions, especially if we want to extend the hazard analysis time span to longer intervals with respect to one century (the default time span adopted by the current impact monitoring systems). More precisely, the non-gravitational model also plays a fundamental role in the Line Of Variations (LOV) computation and propagation (Milani et al. 2000, 2005a,b). If an asteroid with a very well-constrained orbit experiences a very deep close approach, the post-encounter situation is equal to that of a poorly determined orbit; the difference is that the large uncertainty of the asteroid state is due to the divergence of nearby orbits caused by the encounter and not to the poor constraints of the initial conditions. In this case, the initial confidence region is very small, thus the use of the linear approximation of the LOV is allowed. In case such an encounter occurs the linear LOV direction is derived by analysing that encounter and mapping back the corresponding LOV trace on the target plane (TP; Valsecchi et al. 2003) to the space of initial conditions. This method has been used in Spoto et al. (2014) to properly assess the impact risk of (410777) 2009 FD, exploiting its 2185 scattering encounter with the Earth. The same formalism can be used even when we are not in the presence of a scattering encounter, but the close encounter isso deep that the LOV turns out to be very stretched at the next encounter, as in the cases analysed below.

So far, just four asteroids required such special treatment for a proper impact risk assessment, namely (101955) Bennu, (99942) Apophis, (29075) 1950 DA, and (410777) 2009 FD, but this list is expected to grow as a consequence of the work presented inthis paper. Below we show two examples of asteroids for which we found virtual impactors using a non-gravitational model and that have no possible impacts with a purely gravitational model. We are aware that such a work could be performed on many asteroids with accepted Yarkovsky detections, but this is beyond the scope of this paper.

2001 BB₁₆. Currently, this asteroid has a low Minimum Orbit Intersection Distance (MOID) value, ≃ 0.0043 au, but no chance of impacting the Earth in the next century. 2001 BB₁₆ has a deep close approach with the Earth in 2082, which causes an increase of the stretching of two orders of magnitude with respect to the next 2086 encounter, whereas the stretching value remains essentially constant until the 2082 close approach. We used this close approach to derive the LOV direction and we performed the impact monitoring through 2200 employing a non-gravitational model including the Yarkovsky effect. The results are shown in Table 4. In particular we found two VIs at the very end of the 22nd century, which we would not find with a gravity-only model.

2011 MD. This is a very small asteroid, about 6 m in diameter, as determined in Mommert et al. (2014a). In this case as well, the MOID value is very low, ≃ 0.00036 au and it has no virtual impactor in the next century. In 2049, this asteroid will experience two very close approaches with the Earth, causing an increase of two orders of magnitude in the stretching between these encounters and the following approach in 2067. We used thefirst 2049 close approach (the deepest of the two) to compute the LOV direction in the space of initial conditions. We thus performed the impact monitoring using a dynamical model including both the Yarkovsky effect and SRP. The results are shown in Table 5. When we only include SRP, the orbit uncertainty shrinks and thus the number of VIs is much lower than before (see Table 6). Both Tables 5 and 6 list the virtual impactors with IP ≥ 10⁻⁷, since this threshold is the completeness limit used for the LOV sampling (Del Vigna et al. 2018).

It is worth noting that this asteroid is so small that it would not reach the Earth in case of a real impact because it would be burnt in theatmosphere. This case is studied to show that, in some cases, a non-gravitational model is needed to make reliable impact predictions and also that different models of non-gravitational perturbations can give very different results.

9 Conclusions and future work

In this paper, we significantly increased the knowledge of non-gravitational perturbations on NEAs based on actual measurements rather than on modelling. The number of significant and reliable Yarkovsky detections in the NEA catalogue is expected to grow continuously. In fact, the data volume of future surveys, the increased astrometric accuracy for optical observations, more accurate star catalogue debiasing techniques, and expanded efforts in radar astrometry provide ever better constraints to measure this small effect. We identified 87 NEAs with significant and reliable Yarkovsky detections, thus doubling the list provided in Chesley et al. (2016). For few exceedingly small asteroids, we attempted to directly detect SRP together with the Yarkovsky-related acceleration. For such cases, SRP is needed to obtain a more satisfactory orbital fit.

There are several research centres handling the computation of asteroid orbits as an industrial production. These centres are recomputing either all the orbits of more than 500 000 numbered asteroids every time a change in the error model occurs or a large portion of these orbits just to take into account new observations and new asteroid discoveries⁷. There are important scientific goals such as asteroid families and impact monitoring that can only be achieved by maintaining and constantly updating such large lists of orbits.

We dedicated a significant effort to clarifying a number of marginal and/or spurious cases, not only to recover a few dubious cases but also to refine the methodology and therefore be ready for the future increase of significant detections. Indeed, the problem to be faced in the near future is not another increase by a factor two, rather an increase by orders of magnitude. The second Gaiadata release (April 2018) will contain about 1.7 billion of sources that are brighter than magnitude 21 and ≃ 14 000 asteroids with astrometry reaching the sub-milliarcsec accuracy in an optimal range of magnitude G ≃ 12−17 (Gaia Collaboration 2018a,b). The stellar catalogue produced by Gaia will represent the starting point for a new debiasing and weighting scheme. Moreover, the combination of Gaia asteroid observations with the already available observations will produce a sharp increase in the number of objects for which the Yarkovsky effect will be detectable. Thus, the challenge in papers like this is not to establish a new record list of Yarkovsky and/or radiation pressure detections, but rather to develop an automated calculation of orbits with estimated non-gravitational parameters.

The computations of orbits with non-gravitational effects is still very far from an automated process. To avoid spurious detections, we used the most recent error models for the observations and a filtering criterion based on an estimate of the Yarkovsky effect based upon a physical model of the asteroid. Unfortunately, both of these tools are still incomplete. The error models suffer from the continued unavailability of metadata, such as the signal-to-noise of individual observations, with the result that observations with different quality are bundled together and the statistical analysis of the residuals does not yet allow a correct derivation of uncertainty of the measurement error. The physical models of asteroids, which are needed to estimate the expected Yarkovsky effect, are very rough approximations when the main physical data are not available, as is the case for the majority of the asteroids in our tables. Moreover, such small perturbations can be sensitive to old isolated and possibly bad astrometric positions.

In conclusion, we made a step in the right direction by developing and testing the use of different error models and by using the difference in the results as an estimate of the sensitivity of the results on the error model. We developed and tested the use of a filter for spurious cases, which is based on an estimate of the expected Yarkovsky effect, which is roughly the same as the Yarkovsky calibration used to compute the age of asteroid families (Milani et al. 2014; Spoto et al. 2015). Both tools improved our capability to obtain a list of reliable Yarkovsky detections and a much shorter list of radiation pressure detections for natural bodies.

Acknowledgements

We thank the referee Dr. David Vokrouhlicky for his useful comments that have improved the quality of the paper. A. D.V. and L.F. acknowledge support by the company SpaceDyS. D.F. conducted this research at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. This research was conducted under European Space Agency contract No. 4000113555/15/DMRP “P2-NEO-II Improved NEO Data Processing Capabilities”.

Appendix A Additional tables

References

All Tables

All Figures

	Fig. 1 Graphical representation of the partition of the detections set into classes on the plane $(S/N_{A_2},\mathcal{S})$. We plot theaccepted detections (green dots), the marginally significance detections (blue dots) plus the special cases (410777) 2009 FD and (99942) Apophis (blue asterisk), and the rejected detections (red crosses).
In the text