© The Authors 2024

1 Introduction

The Yarkovsky effect is a nongravitational force caused by thermal emission. It affects the dynamics of asteroids whose diameter is smaller than about 40 km (see e.g. Vokrouhlický 1998; Bottke et al. 2006; Vokrouhlický et al. 2015a), and it mainly manifests as a drift in the semi-major axis of the object. It is known to cause several phenomena in the dynamics of small bodies in the Solar System, such as the spreading of asteroid families (see e.g. Spoto et al. 2015), and the transport of objects from the main belt to the near-Earth region (see e.g. Granvik et al. 2017). From the point of view of planetary defence, measuring nongravitational effects acting on near-Earth asteroids (NEAs) is important for a reliable assessment of the impact risk, especially for the long-time horizon search of potential impacts (Chesley et al. 2014; Farnocchia & Chesley 2014; Spoto et al. 2014; Vokrouhlický et al. 2015b).

Additionally, measurements of the Yarkovsky effect can be used to extrapolate the physical properties of the asteroids (see e.g. Chesley et al. 2014; Rozitis et al. 2014; Fenucci et al. 2021, 2023), which are key parameters to know for the design of lander and sample-return spacecraft missions (Murdoch et al. 2021), or to predict the efficiency of a kinetic impactor-deflection mission (Cheng et al. 2023; Daly et al. 2023; Thomas et al. 2023).

Since the first measurement of the Yarkovsky effect on asteroid (6489) Golevka (Chesley et al. 2003), many detections on NEAs have been found and announced in a series of related works (Nugent et al. 2012; Farnocchia et al. 2013; Del Vigna et al. 2018; Greenberg et al. 2020). Today, the information about the detections is regularly provided by the three main centers for asteroid orbit computation: the ESA NEO Coordination Centre¹ (NEOCC), the NEODyS² service, and the NASA JPL Solar System Dynamics³ (SSD) group. However, the process of detecting new measurements is not fully automated, and the identification of spurious detections still relies on manual checks that are often based on personal judgment. There are several reasons that make a detection spurious, such as dynamical modeling issues or poor-quality old astrometry. Some detections may also appear to be incompatible with the Yarkovsky effect, but they are in fact real and are caused by other poorly understood phenomena, as in the case of asteroid (523599) 2003 RM (Farnocchia et al. 2023). Therefore, in some cases, it is questionable whether a Yarkovsky effect has really been detected or not. On the other hand, the new and improved observational technologies and the high-quality data provided by the ESA Gaia mission (Tanga et al. 2023) mean that the number of detected Yarkovsky effects is only expected to increase in the future. This will become especially true after the beginning of the operational phase of new-generation ground-based surveys, such as the ESA Flyeye telescope (Conversi et al. 2021) and the Vera Rubin Observatory (Ivezić et al. 2019), and of the space-based infrared telescopes NEO Surveyor by NASA (Mainzer et al. 2021) and NEOMIR by ESA (Conversi et al. 2023). These new facilities are expected to increase the number of known asteroids by roughly an order of magnitude (Jones et al. 2018). This situation calls for the development of new and automated methods for the detection of the Yarkovsky effect and for the identification of spurious measurements.

In contrast to previous efforts, we propose a systematic procedure here to search for candidates for measurements of the Yarkovsky effect on NEAs, and for the identification of spurious detections. We use a statistical model of the physical parameters of the asteroids that relies on the most recent NEA data sets and population models (Granvik et al. 2018; Morbidelli et al. 2020; Berthier et al. 2023). We first estimate the expected Yarkovsky effect with a Monte Carlo method. Then, a list of NEA candidates for the detection of the Yarkovsky effect is determined according to the expected semi-major axis drift, the observational arc, and the uncertainty of the current orbit. For each object of the list, a nongravitational force along the transverse direction is estimated through orbit determination, and only detections with a good signal-to-noise ratio that are statistically compatible with the prediction from the physical model are kept. Then, we determine the dependence of the detections on isolated tracklets, and only those that passed the check are accepted. Among known NEAs, the algorithm identified 348 detections. The procedure described here is now adopted by the ESA NEO Coordination Centre, and data are automatically uploaded on the web portal.

The paper is organized as follows. In Sect. 2, we describe the model for the Yarkovsky effect and the model of the physical parameters of the NEAs that we used to predict the semi-major axis drift, the method we used to determine the orbit, and the procedure we used to automatically detect the Yarkovsky effect. In Sect. 3, we provide the results we obtained by running the procedure on the known NEA population, we compare them with the results reported by the JPL SSD group, and we discuss the effectiveness and limitations of our algorithm. Finally, we summarize our conclusions in Sect. 4.

2 Methods

2.1 Orbit determination

The dynamical model we used to determine the orbit of an NEA includes the gravitational perturbations of the Sun, the eight planets, the Moon, the 16 most massive main-belt asteroids, and Pluto. The masses and positions of all these bodies were all retrieved from the JPL ephemeris DE441 (Park et al. 2021). The masses of the 16 main-belt asteroids included in the dynamical model and the mass of Pluto are given in Table 1. Relativistic effects of the Sun, the planets, and the Moon were also added to the dynamical model as a first-order Newtonian expansion (Will 1993). The Yarkovsky effect is expressed as an acceleration along the transverse direction $\widehat{t}$ of the form $$a_t=A_2{\left(\frac{1\text{au}}{r}\right)}^2\widehat{t},$$(1)

where r is the distance from the Sun (Farnocchia et al. 2013). The corresponding semi-major axis drift is obtained through the Gauss planetary equations by $$\frac{\text{d}a}{ \text{d}t}=\frac{2A_2\left(1-{\text{e}}^2\right)}{n{p}^2},$$(2)

where n is the mean motion, and p = a(1 − e²) is the semilatus rectum. Additionally, it is possible to optionally add solar radiation pressure (SRP) to the dynamical model. The SRP is modeled as a radial force a_r of the form $$a_r=A_1{\left(\frac{1\text{\hspace{0.17em}au}}{r}\right)}^2\widehat{r},$$(3)

where $\widehat{r}=r/r$ is the radial unit vector, and r is the heliocentric position of the asteroid. The coefficient A₁ is computed as in Montenbruck & Gill (2000), $$A_1=(1-\text{Θ})\times \frac{\varphi }{c}\times \frac{A}{m},$$(4)

where ϕ = 1.361 kW m⁻² is the solar radiation energy flux at 1 au, c is the speed of light, A/m is the area-to-mass parameter, and Θ is the occultation function, which determines whether the object is occulted by a major planet.

The parameter A₂ (and possibly the A/m parameter, when SRP is included in the model) is determined together with the orbital elements by fitting the model to the observations through a least-squares procedure (Milani & Gronchi 2009). Observational outliers were identified and discarded with an automatic rejection algorithm developed by Carpino et al. (2003). In all the computations, we used the weighting scheme based on observatory statistics developed by Vereš et al. (2017), and the debiasing scheme based on star catalogs developed by Farnocchia et al. (2015). The orbital elements estimated from the fit refer to the weighted mean of the observation times, as indicated in Milani & Gronchi (2009). To perform the orbital fit, we used the ESA Aegis orbit determination and impact monitoring software (Faggioli et al. 2023), developed and maintained by SpaceDyS s.r.l.⁴ under ESA contracts.

2.2 Expected Yarkovsky drift

2.2.1 Semi-analytical formula for the semi-major axis drift

The instantaneous semi-major axis drift produced by the Yarkovsky effect can be expressed by analytical formulas (Vokrouhlický et al. 2017), assuming a spherical body and a linearization of the boundary conditions. It is given by $$\frac{\text{d}a}{ \text{d}t}=\frac{2}{n^{2}a}{f}_{\text{Y}}\cdot v,$$(5)

where a is the semi-major axis of the asteroid orbit, v is the heliocentric orbital velocity, and f_Y is the instantaneous value of the Yarkovsky acceleration. The term f_Y is the result of the sum of the diurnal component f_Y,d and the seasonal component f_Y,s. The diurnal component is expressed as $$f_{\text{Y},\text{d}}=\kappa \left[(\widehat{r}\cdot \widehat{s})\widehat{s}+{\gamma }_1(\widehat{r}\times \widehat{s})+{\gamma }_2\widehat{s}\times (\widehat{r}\times \widehat{s})\right].$$(6)

In Eq. (6), $\widehat{r}=r/r$ is the heliocentric unit position vector, and ŝ is the unit vector of the asteroid spin axis. In addition, $$\kappa =\frac{4\alpha }{9}\frac{SF}{mc},$$(7)

where S = πR² is the cross section of the asteroid, R is the radius, F is the solar radiation flux at a heliocentric distance r, α is the absorption coefficient, m is the asteroid mass, and c is the speed of light. The coefficients γ₁, γ₂ are expressed as $$\begin{array}{l}{\gamma }_1=-\frac{k_1\left({R^{\prime }}_{\text{d}}\right){\text{Θ}}_{\text{d}}}{1+2k_2\left(R^{\prime }\right){\text{Θ}}_{\text{d}}+k_3\left({R^{\prime }}_{\text{d}}\right){\text{Θ}}_{\text{d}}^2},\\ {\gamma }_2=-\frac{1+k_2\left({R^{\prime }}_{\text{d}}\right){\text{Θ}}_{\text{d}}}{1+2k_2\left(R^{\prime }\right){\text{Θ}}_{\text{d}}+k_3\left({R^{\prime }}_{\text{d}}\right){\text{Θ}}_{\text{d}}^2}.\end{array}$$(8)

In Eq. (8), R′ = R/l_d is the scaled radius, where $l_{\text{d}}=\sqrt{K/(\rho C\omega )}$ is the penetration depth of the diurnal thermal wave, obtained from the thermal conductivity K, the heat capacity C, the density ρ, and the rotation frequency ω. The functions k₁, k₂, and k₃ are positive analytic functions of the scaled radius (Vokrouhlický 1998, 1999). In addition, ${\text{Θ}}_{\text{d}}=\sqrt{\rho KC\omega }/\left(\epsilon \sigma T_{\star }^3\right)$, where σ is the Stefan-Boltzmann constant, and T_★ is the subsolar temperature, which is given by $4\epsilon \sigma T_{\star }^4=\alpha F$.

The seasonal component is expressed as $$f_{\text{Y},\text{s}}=\kappa \left[{\overline{\gamma }}_1(\widehat{r}\cdot \widehat{s})+{\overline{\gamma }}_2(\widehat{n}\times \widehat{r})\cdot \widehat{s}\right]\widehat{s},$$(9)

where $\widehat{n}$ is the unit vector normal to the orbital plane. The coefficients ${\overline{\gamma }}_1,{\overline{\gamma }}_2$ have the same expressions as Eq. (8), but evaluated with ${\overline{\text{Θ}}}_{\text{s}}=\sqrt{\rho KCn}/\left(\epsilon \sigma T_{\star }^3\right)$, and with a scaled radius ${R^{\prime }}_S=R/l_S$, where $l_{\text{s}}=\sqrt{K/(\rho Cn)}$ is the penetration depth of the seasonal thermal wave. The average Yarkovsky drift da/dt is finally obtained by numerically averaging the instantaneous Yarkovsky drift of Eq. (5) over an orbital period.

2.2.2 Modeling the orbital and physical parameters

The orbital parameters in the analytical formulas of Eq. (5) are the semi-major axis a and the eccentricity e. The uncertainties on these parameters are usually small when the orbit is well determined, and they typically do not affect the results given by Eq. (5). We therefore fixed them to their nominal value.

The other parameters in Eq. (5) are the diameter D, the density ρ, the obliquity γ, the rotation period P, the heat capacity C, the thermal conductivity K, the emissivity ɛ, and the absorption coefficient α. The emissivity was fixed to ɛ = 0.984, which corresponds to the average value measured on meteorites (Ostrowski & Bryson 2019). The heat capacity was fixed to a value of C = 800 J kg⁻¹ K⁻¹, which is appropriate for NEAs (Delbo’ et al. 2015). However, changing this value in the range 500–1200 J kg⁻¹ K⁻¹ , which is thought to be valid for most of materials present on asteroids, produces very small changes in the predicted semi-major axis drift (Fenucci et al. 2021, 2023). Because we are interested in estimating the maximum semi-major axis drift that can be reached by a NEA, we assumed that the obliquity γ can be either 0° or 180°. For the same reason, we also fixed α = 1.

To obtain the values of the other physical parameters, we used the SsODNet service⁵ (Berthier et al. 2023) through the Python rocks API⁶. The diameter D, the rotation period P, and the taxonomic class were retrieved through ssoCard, which provides a best estimate by averaging all the available measurements of a given parameter. When the values of D or P were provided, we assumed that they are Gaussian distributed with a mean value equal to the nominal value and a standard deviation equal to the error. On the other hand, the taxonomic complex gives only indirect information about the density. Table A.1 reports all the taxonomic complexes found for NEAs in the SsODNet database. The grouping of the different taxonomic types (Tholen 1984; Bus & Binzel 2002; DeMeo et al. 2009) and complexes can be found in Berthier et al. (2023). The density was assumed to be distributed according to a log-normal distribution, with values of µ and σ provided by Table A.1. These values were obtained by a statistics on the density values extracted from the SsODNet database (see Appendix A).

Whenever the diameter D, the taxonomic complex, or the rotation period P were lacking in the SsODNet database, we modeled them according to the properties and available models of the whole NEA population. The distribution of the measured rotation periods shows a dependence on the diameter (Pravec & Harris 2000), where objects rotating in less than about 2.2 h are found only at sizes with diameters smaller than roughly 150– 200 m. Here we assumed that the rotation period varied with the diameter as $$P(D)=\{\begin{array}{ll}5.6 \text{h}\hfill & \text{ if }D\ge 200 \text{m},\hfill \\ 5.6 \text{h}\times {\left(\frac{D}{0.2 \text{km}}\right)}^2\hfill & \text{ if }D<200 \text{m}.\hfill \end{array}$$(10)

Figure 1 shows the trend of the rotation period of Eq. (10) we adopted, superimposed on the distribution of rotation period measurements of NEAs extracted from the Asteroid Lightcurve Database (LCDB; Warner et al. 2009). Above the threshold of 200 m of diameter, we used a constant rotation period of 5.6 h, corresponding to the median value of the measured rotation periods of asteroids larger than the threshold. At sizes smaller than 200 m, asteroids may rotate very fast without being disrupted, and therefore, we assumed a trend that scales as D⁻², which is a better fit to the data than the D⁻¹ trend that was previously suggested by Farinella et al. (1998). The biases in the rotation period measurements for asteroids with diameters smaller than 200 m are unclear, and it is difficult to extrapolate a general trend for the whole population. However, we found this solution to be the best that can be done based upon the currently available observations.

For the diameter D and the density ρ, we used the model by Fenucci et al. (2021, 2023), which combines the orbital distribution model for NEAs by Granvik et al. (2018) and the albedo distribution model for NEAs by Morbidelli et al. (2020). The model first produces a distribution of the visual albedo p_V that is later converted into a joined distribution (D,ρ) of the diameter and density. For this purpose, we also assumed that the absolute magnitude H is Gaussian distributed, with a standard deviation given by the root mean square (rms) error obtained by converting the visual magnitude measurements into absolute magnitude. When D is given in the SsODNet database, but not the taxonomic complex (or when the taxonomic complex is given in the SsODNet database, but not the diameter D), we used the marginal distribution of ρ (the marginal distribution of D) obtained by the same model.

The thermal conductivity K is related to the thermal inertia Γ by the relation $$\text{Γ}=\sqrt{\rho KC}.$$(11)

We assumed Γ to be Gaussian distributed with a mean value of 250 J m⁻² K⁻¹ s^−1/2 (Delbo’ et al. 2007) and a standard deviation of 100 J m⁻² K⁻¹ s^−1/2, and the value of K was then computed by inverting Eq. (11).

Table 2 summarizes the physical parameters for the model of the Yarkovsky effect in Eq. (5) and the corresponding value or model we used to compute the predicted drift of the semi-major axis.

Fig. 1 Distribution of the rotation period measurements of NEAs extracted from LCDB, together with the rotational disruption limit of 2.2 h and the dependence of P(D) on the diameter.

2.3 Treatment of poor astrometry

The fact that a high signal-to-noise ratio (S/N) detection of A₂ is found from orbit determination does not imply that it is due to a true dynamical effect. In fact, previous works (Farnocchia et al. 2013; Del Vigna et al. 2018) found that old and poorly measured astrometry or isolated tracklets that significantly increase the observational arc (see also Hu et al. 2023) may lead to erroneous detections of the Yarkovsky effect. To treat old astrometry, we adopted data weights at 10 arcsec for measurements taken before 1890 and at 5 arcsec for those taken before 1950, so that very old observations were less significant in the orbit determination process.

To ensure that a detection with an S/N higher than 3 is reliable and does not depend on short isolated tracklets, we adopted the following strategy. We defined an isolated tracklet as a set of observations such that (1) they were all taken from the same observatory, (2) the length of the arc was shorter than 15 days, and 93) they were separated by at least 5 yr from the previous and the following observation sets. When an isolated tracklet was found, we removed it from the dataset of observations, and we repeated the process until a dataset without isolated tracklets was obtained. Then, we performed the 7-D OD with the dataset that did not contain isolated tracklets. When the S/N in A₂ was lower than 3, the detection strongly depends on the old isolated observations and may therefore not be reliable or possibly spurious. As a conservative choice, the detection was therefore discarded. On the contrary, if a signal with a large S/N in A₂ is found also with the smaller set of observations and with a possibly shorter observational arc, the signal is most likely true and due to dynamical effects rather than astrometric issues. In this case the detection is considered valid.

Note that detections are more sensitive to isolated track-lets appearing at the beginning of the observational arc, as they extend the arc length, however they could be also affected by intra-arc isolated tracklets. In addition to this, some isolated tracklets may be known to be reliable. For example, some old observations may have been carefully remeasured from an experienced observer who examined the detection and confirmed its authenticity and the correctness of the astrometry (see Vokrouhlický et al. 2008, for an example), or they may be produced by modern surveys using more advanced technologies. To handle these cases, we introduce a list of exceptions to the rule described above, where the isolated tracklet is always included in the orbital fit. Regarding modern surveys, we assume that only tracklets from the Pan-STARRS survey (IAU codes F51 and F52) are always included in the orbital fit, as a conservative choice.

2.4 Automated Yarkovsky effect detection

To select candidates for Yarkovsky effect detection, we first discard all the NEAs with uncertainty in semi-major axis σ_a larger than 10⁻⁵ au, since their orbit is not known with a sufficient accuracy to attempt for the detection. Here σ_a refers to the uncertainty of the orbit obtained from 6-dimensional orbit determination (6-D OD), estimated at the weighted mean observational epoch. For the asteroids passing the first criterion, we compute the distribution of the predicted maximum Yarkovsky drift⁷ (da/dt)_max that each of them can be subjected to. This is done with a Monte Carlo method by computing the expected Yarkovsky drift with the model of Sect. 2.2.1, using 500 000 combinations of input parameters chosen from the parameter distributions described in Sect. 2.2.2. Then, we determine the 95-th percentile of the distribution of |(da/dt)_max|, i.e. the value Y_M such that $$P\left(\left|{(\text{d}a/\text{d}t)}_{\max }\right|<Y_M\right)=0.95,$$(12)

where P(⋅) denotes the probability. Values of semi-major axis drifts larger than Y_M in absolute values are attained only in very unlikely combinations of the input parameters. By denoting with ∆T the total length of the observational arc of the asteroid, if the condition $${\sigma }_a<\text{Δ}T\times Y_M$$(13)

is fulfilled, the asteroid is considered a good candidate for Yarkovsky effect detection, hence the value of A₂ is fitted through a 7-dimensional orbit determination (7-D OD). If the signal-to-noise ratio $S/N=\left|A_2/{\sigma }_{A_2}\right|$ is smaller than 3, then the detection is considered not statistically significant, and therefore discarded.

To understand whether a detection is spurious or if the value of A₂ obtained from orbit determination could be actually caused by the Yarkovsky effect, we determine whether A₂ is statistically compatible with the physical properties of the given object or not. To this purpose, we convert A₂ and its uncertainty ${\sigma }_{A_2}$ to semi-major axis drift values (da/dt)_fit and ${\sigma }_{{(\text{d}a/\text{d}t)}_{\text{fit}}}$ respectively, by using Eq. (2). Then, we check if the value (da/dt)_fit obtained from orbit determination fulfills the condition $$\left|{(\text{d}a/\text{d}t)}_{\text{fit}}\right|-k{\sigma }_{{(\text{d}a/\text{d}t)}_{\text{fit}}}<Y_M,$$(14)

where k = 1.645 gives the 5th percentile of the measured Yarkovsky drift. If inequality of Eq. (14) is fulfilled, it means that semi-major axis drift values down to $\left|{(\text{d}a/\text{d}t)}_{\text{fit}}\right|-k{\sigma }_{{(\text{d}a/\text{d}t)}_{\text{fit}}}\mid$ are compatible with the predictions obtained with the Yarkovsky model. On the other hand, if inequality of Eq. (14) is not valid, the value detected from astrometry is not compatible with the Yarkovsky physical model, the detection is considered spurious and therefore discarded. Note that the criterion of Eq. (14) is adopted because the distribution of the expected maximum Yarkovsky drift is not Gaussian in general. Finally, for those detections fulfilling Eq. (14) we check the dependence on isolated tracklets by using the method described in Sect. 2.3, and keep only those detections passing the test.

Fig. 2 Distribution of the rms of the normalized residuals for the NEAs with a positive Yarkovsy effect detection. The blue (red) histogram represents the rms of the normalized residuals obtained with (without) the Yarkovsky effect included in the dynamical model. The dashed vertical lines correspond to the median of the two distributions, where the numerical value is reported in the legend.

3 Results and discussion

3.1 Accepted detections

We ran the above procedure on the NEOCC list of known NEAs on 7 August 2023, soon after the last monthly orbital update issued by the Minor Planet Center⁸ (MPC). Optical observations are downloaded from the MPC, while radar observations are taken from the JPL Solar System Dynamics website⁹. Out of 32 396 NEAs, 5 308 were selected as candidates for a possible Yarkovsky effect detection. After the orbit determination process, we found 461 NEAs with S/N ≥ 3, and 373 of them fulfilled also the condition of Eq. (14). Among them, 69 objects presented isolated tracklets, and 25 were found to be spurious detections. In conclusion, we found 348 positive Yarkovsky effect detections. A short list of detections with high S/N is shown in Table B.1, while the complete list of accepted detections in extended precision can be found at the CDS. Figure 2 shows the distribution of the rms error of the normalized observational residuals of the accepted detections, obtained with 7-D OD (blue histogram). As a comparison, the red histogram shows the distribution of the rms of normalized residuals of the same NEAs obtained with the 6-D OD. The median value decreases from 0.616 for the 6-D OD, to 0.548 for the 7-D OD. The distribution obtained with 7-D OD is skewed towards smaller values if compared to that obtained with 6-D OD, meaning that including the Yarkovsky effect in the orbital fit of these NEAs generally improves the quality of the fit of the accepted detections.

Figure 3 shows the distribution of the S/N of the accepted Yarkovsky effect detections. About 110 objects have S/N between 3 and 4, and the presence of radar observations generally improves the S/N, even though they are not necessarily needed to obtain good detections. The asteroids with S/N > 40 are (101955) Bennu, 1998 SD₉, (99942) Apophis, (480883) 2001 YE₄, (524522) 2002 VE₆₈, (2340) Hathor, and 2002 BF₂₅, which are not shown in Fig. 3.

The Yarkovsky effect determinations can be used to extrapolate the ratio of retrograde-to-prograde rotators (R/P), because negative (positive) semi-major axis drift values are associated with retrograde (prograde) rotators. To this purpose, we randomly generate 10 000 values of A₂, for each positive detections. In this process, we assumed A₂ to be a Gaussian random variable, with mean value equal to the nominal value and standard deviation equal to the 1-σ uncertainty of the detection. With this method, we obtained R/P = 2.976. Spin pole determinations of NEAs show that ${\text{R/P=2}}_{-0.7}^{+1}$ (La Spina et al. 2004), in good agreement from what NEA distribution models predict. Although the ratio R/P we found is still consistent with spin pole determinations, the discrepancy with the distribution models can possibly be explained by the biases in the orbital distribution of NEAs with a measured Yarkovsky effect (Farnocchia et al. 2013). Another possible cause contributing to this discrepancy might be that the Yarkovsky detections are driven by cases with extreme obliquity, which may have a larger R/P ratio than the general NEA population. However, it is worth noting that the percentage of retrograde rotators found here is about 75%, which is lower than the 81% found in Farnocchia et al. (2013), and closer to the 69% value expected in the NEA population. This is an indication that the observation capabilities improved in the last 10 yr, and the new generation surveys may help in making the sample of NEAs with measured Yarkovsky effect representative of the whole debiased NEA population. Note that a value of ${\text{R/P=2}}_{-0.7}^{+0.8}$ consistent with observations was previously found in Tardioli et al. (2017), by fitting the NEA obliquity distribution to the Yarkovsky effect determinations.

Figure 4 shows the measured semi-major axis drifts |(da/dt)_fit| as a function of the diameter D. Error bars in the diameter were obtained by computing the 15th and the 85-th percentile of the diameter distribution, using the model described in Sect. 2.2.2. Analytical models of the Yarkovsky effect predict that da/dt scales with the diameter as D⁻¹. To verify that spurious detections are correctly discarded by our procedure, we fit the data about accepted detections with a function of the form α × D^β using the Orthogonal Distance Regression (ODR; Boggs et al. 1987) method. The best fit values computed by the ODR algorithm are α = (31.65 ± 2.42) × 10⁻⁵ and β = −0.97 ± 0.03. Note that the exponent β is statistically compatible with the value of −1 predicted by Yarkovsky analytical models within 1-σ. The power law function computed with the nominal parameters is also shown in Fig. 4. Detections with S/N ≥ 3 that are considered spurious are also shown in Fig. 4 with red points, and most of them are far off the best-fit power law computed from the accepted detections. This already shows that our method is able to appropriately find spurious detections. The dispersion from the best fit power law of Fig. 4 is explained by the differences in the physical properties of asteroids, that are able to slightly modify the magnitude of the semi-major axis drift. To verify this we retrieved taxonomic data from the SsODNet database about the accepted Yarkovsky effect detections, finding 116 NEAs with such information. A histogram showing the number of NEAs with measured Yarkovsky effect for each taxonomic complex is shown in Fig. 5. As expected, the vast majority of the detections belong to the S-complex, however almost all the complexes (see Table A.1) are represented in the Yarkovsky measurement sample. These facts indicate that our procedure is able to discriminate between positive and negative Yarkovsky effect detections from a statistical point of view.

Fig. 3 Distribution of the S/N of the accepted Yarkovsky effect detections, obtained with equi-spaced bins of width equal to 1. NEAs with only optical observations are represented in the blue histogram, while the red histogram refers to NEAs observed also with radar.

Fig. 4 Distribution of the accepted Yarkovsky effect detections in the plane (D, |(da/dt)fit|). Green points are NEAs for which an estimate of the diameter is available, while black points are those for which the physical model is used. The error bars in the diameter are represented by using the 15th and 85th percentile of the diameter distribution (see Sect. 2.2.2), which roughly correspond to the 1-σ uncertainty in the case the diameter is estimated. The error bars in the semi-major axis drift represent the 1-σ uncertainty obtained by the conversion of A2, that has been estimated with orbit determination. The dashed blue line represents the fit with a function of the form α× Dβ. Red points are detections with S/N ≥ 3 that are considered spurious, hence rejected by our procedure.

3.2 Notable rejections

To further show the effectiveness of the algorithm introduced here, we highlight some notable cases among the identified spurious detections. In addition to the Yarkovsky effect, other nongravitational forces such as SRP can influence the dynamics of very small bodies. So far, SRP has been detected through orbit determination only on a short list of NEAs: 2009 BD (Micheli et al. 2012), 2011 MD (Micheli et al. 2014), 2012 LA (Micheli et al. 2013), 2012 TC₄, and 2015 TC₂₅ (Del Vigna et al. 2018). All these objects do not appear in our list of accepted Yarkovsky effect detections, because the semi-major axis drift obtained by the 7-D OD is not compatible with what expected from the physical model. This issue is actually fixed when SRP is added to the dynamical model, and the area-to-mass ratio A/m of the object is added to the list of parameters to be estimated in the orbital fit.

A similar result has been found for 2021 GM₁, an NEA that is currently in the NEOCC Risk List¹⁰. The Yarkovsky effect detection of this object is rejected by our procedure, but a nongravitational effect is needed to better fit the observations. Adding SRP to the dynamical model, we fitted an area-to-mass ratio of 0.2564 ± 0.00257 kg² t⁻¹ (see Table 3 for the complete orbital solution), providing a new SRP detection that was not available in the literature¹¹. This shows that the procedure may also be helpful in finding asteroid candidates for detection of other nongravitational effects.

Other interesting cases are 2006 RH₁₂₀ (Seligman et al. 2023) and (523599) 2003 RM (Farnocchia et al. 2023). In addition to the radial component A₁ and the tangential component A₂, it is possible to also fit a nongravitational out-of-plane force component A₃ for these two NEAs, showing a dynamical behavior similar to that of comets. In these two cases, the value of A₂ obtained by a 7-D OD only is not compatible with the Yarkovsky physical model (see also Table B.3), while it is likely due to outgassing (Chesley et al. 2016; Taylor et al. 2024). Still, these detections are recognized as spurious by our procedure.

Another remarkable detection that was flagged as spurious is (433) Eros. This NEA was selected as a candidate for Yarkovsky effect detection because of its long observational arc of more than 100 yr, and A₂ was fitted with a S/N of about 5. However, the detection was not statistically compatible with the prediction made by the Yarkovsky physical model, and therefore the detection is spurious. This detection needs to be necessarily discarded, because bad quality astrometric measurements performed before the 1950 cause an artificial nongravitational acceleration that is actually not expected. In fact, (433) Eros is known to have an obliquity of almost 90° (Yeomans et al. 2000), that minimizes the Yarvkosky effect, but the procedure is able to discard the detection even by assuming an obliquity of 0° or 180° (see Sect. 2.2.2). All these considerations provide evidences of the effectiveness of the algorithm introduced here, in both the selection of possible candidates for the Yarkovsky effect determination, and in the identification of spurious detections.

The algorithm found 69 A₂ detections with S/N larger than 3 and fulfilling Eq. (14) that were affected by the presence of isolated tracklets in their dataset. Among them, 58 presented isolated tracklets at the beginning of the observational arc, and the complete list of these asteroids is reported in Table B.4. The table reports also the S/N of A₂ obtained without including the isolated tracklet in the orbital fit, following the procedure described in Sect. 2.3, and 30 of them still present a signal in A₂ even with the shorter observational arc. These objects were therefore accepted as positive Yarkovsky effect detections. Regarding the list of exceptions for which we know that astrometry is reliable, up to now it includes: 2012 BB₁₂₄ and 2014 RO₁₇, for which accurate measurements were performed by MM; 2015 LA₂, 2018 KR, and 2016 UW₈₀, that present isolated tracklets from F51; (152563) 1992 BF for which it is known that the tracklet has been accurately remeasured (Del Vigna et al. 2018). While this list may increase in the future, it is not expected to grow to a point that cannot be maintained manually. Other 11 objects presented intra-arc isolated tracklets, but the A₂ detection was affected only for 2005 TG₅₀ and 2011 XC, that were therefore considered as spurious. In addition, the S/N drops down also for 2009 OS₅ when observations from 2014 from Mauna Kea Observatory (IAU code 568) are excluded in the fit, however we considered them as reli-able¹² and we therefore accepted this detection. Finally, among the detections that were not accepted, we note that the algorithm flagged (469219) Kamo’oalewa as spurious, in agreement with Hu et al. (2023).

Fig. 5 Taxonomy distribution of the accepted Yarkovsky effect detections. Data about the taxonomic complex were retrieved from the SsODNet database Berthier et al. (2023).

3.3 Publication on the NEOCC portal

The availability of a large number of reliable Yarkovsky effect detections opens possible future developments in other fields of asteroid research, such as physical properties determination (Rozitis & Green 2014; Fenucci et al. 2021, 2023), long term dynamical studies (Fenucci & Novaković 2021), or population studies on NEAs (Tardioli et al. 2017). To this end, it is important to continuously provide up-to-date data to the scientific community. The procedure described here has been now adopted by the ESA NEOCC to periodically update the Yarkovsky effect detections on NEAs, that are automatically posted on the NEOCC web portal¹³ and available through the HTTPS APIs¹⁴. The catalogues of observations and of known NEAs are continuously growing, and the algorithm is scheduled to run after every monthly orbital update issued by the MPC, so that old determinations are updated and possible new detections are added to the database.

3.4 Comparison with the JPL SBDB

We downloaded the list of determinations of A₂ from the JPL Small Body Database¹⁵ (SBDB) maintained by the JPL SSD group on 8 August 2023. A total of 333 measurements, 308 of which with S/N ≥ 3, were found. Differences in the number of detections with a good S/N are possibly due to different debiasing and weighting schemes of observations, different observational datasets and rejections, and timing errors handling. These differences become particularly significant when the S/N is close to the threshold of 3. Among the JPL detections, 13 objects have also a determination of the component A₁, and 6 objects have a determination of both A₁ and A₃. We did not make the comparison for these objects, because the results are obtained with different dynamical models. Of the remaining objects, 270 were found in our list of accepted Yarkovsky effect measurements. Figure 6 shows the scatter plot of the common A₂ determinations, together with their 1-σ uncertainty. We fitted the data with a linear function of the form α|A₂| + β with the ODR algorithm, that gave an estimation of the parameters of α = 1.00 ± 0.0001 and β = (−1.90 ± 0.97) × 10⁻¹⁶, meaning that our determinations and those of the JPL SBDB are statistically compatible, and no clear outlier can be recognized. This is certainly due to the fact that the orbit determination algorithms and the datasets used by the JPL SSD group are similar to those used here for these cases. To make a quantitative comparison of the results for each asteroid, we computed the relative errors $${\epsilon }_1=\frac{\left|A_2-A_2^{\text{JPL}}\right|}{{\sigma }_{A_2}},{\epsilon }_2=\frac{\left|A_2-A_2^{\text{JPL}}\right|}{{\sigma }_{A_2^{\text{JPL}}}},$$(15)

where the superscript JPL denotes the value provided by the JPL SBDB. For the vast majority of the asteroids in common, the relative errors ɛ₁ , ɛ₂ were both smaller than 3, indicating that the results provided by us and by the JPL SBDB are all statistically comparable, even case by case. The only cases in which the relative errors are larger than 3 are (65717) 1993 BX₃, 2003 AF₂₃, (1685) Toro. However, also in these cases the values of ɛ₁ and ɛ₂ are only slightly larger than 3. We also performed a test by checking that $$\chi =\frac{\left|A_2-A_2^{\text{JPL}}\right|}{\sqrt{{\left({\sigma }_{A_2}\right)}^2+{\left({\sigma }_{A_2^{\text{JPL}}}\right)}^2}}<2,$$(16)

meaning that the hypothesis that the two distributions are not the same is rejected at 95% confidence level. For all the common asteroids with a Yarkovsky effect measurement we found χ < 2, with the only exceptions of 2003 AF₂₃ and (65717) 1993 BX₃, that have χ = 2.23 and χ = 2.18, respectively. However, the two measurements are still considered the same if the test is verified with a more relaxed confidence level of 90%. In addition, we found that the vast majority of the measurements have χ < 1, implying an excellent agreement between our results are those reported in the JPL SBDB. Table B.2 reports the numerical values of the relative errors ɛ₁ and ɛ₂, and of χ for all the NEAs with χ > 1. The full table in extended precision, used to generate Fig. 6, is available at the CDS.

Asteroids with a determination of A₂ that are reported in the JPL SBDB, but are either discarded or not found by our procedure, are 54. Among them, 44 have S/N ≥ 3 and their list is shown in Table B.3, while the complete list in extended precision is available through the CDS. Notable cases are those of (65803) Didymos, (3200) Phaethon, and (25143) Itokawa.

For asteroid (65803) Didymos the JPL SBDB reports an orbital solution with S/N = 16 in the A₂ parameter, while the value we obtained here is only of 2.57. The reasons for this discrepancy are probably due to the availability of some occultation observations took during the September 2022 campaign (Dunham et al. 2023), the availability of data from the Double Asteroid Redirection Test (DART, Cheng et al. 2018) mission, and the inclusion of Gaia astrometry in the orbit computation. These high accuracy observations are not reported by the MPC, and therefore they were not considered in our data set.

For (3200) Phaethon the JPL SBDB reports an orbital fit with a S/N in A₂ of 8.3. The Yarkovsky effect measurement for this object has been reported by Hanuš et al. (2018), where the authors included also very precise Gaia astrometric observations taken between 2014 and 2016 in the dataset. After the publication of that work, occultation observations taken on 25 October 2019 permitted to further improve the orbit of (3200) Phaethon (Dunham et al. 2021). However, both the Gaia astrometry and the occultations are currently not reported by the MPC, and we believe this being the reason we are not able to fit the A₂ parameter with such high accuracy.

The last case is that of (25143) Itokawa, for which the JPL SBDB gives an orbital solution with S/N of 3.2 in the A₂ parameter. In this case, it is not easy to identify what may cause this discrepancy in the results, and the only noticeable difference that we could find is in the number of observations used in the orbital fit. For other objects, we also obtain a S/N larger than 3 as for the JPL SBDB, but they are discarded because the drift is too large to be explained by the Yarkovsky effect, and therefore Eq. (14) does not hold.

Fig. 6 Distribution of the common A2 measurements found in this work and at the JPL SBDB. The blue dashed line is the diagonal.

3.5 Limitations of the procedure

A possible limitation of the procedure may be due to the fact that the Yarkovsky effect modeling takes into account the orbital elements of the NEA at a certain epoch to compute the expected drift da/dt. However, this may have consequences if an object has a close encounter with a planet that is able to change its orbital elements. While this is certainly true, the changes induced in the 95-th percentile Y_M are negligible because the distribution of the expected semi-major axis drift is mostly controlled by the uncertainties in the physical properties. To show this, we took into account the case of (99942) Apophis, for which we have an estimate of the diameter, a measurement of the rotation period, and we know that the spectral type is S. Therefore, the physical model of (99942) Apophis has certainly more constraints than the generic NEA. To show that changes in the orbital elements have a little influence on the results, we computed Y_M with the orbital elements computed at the mean epoch, and with those obtained by propagating the nominal orbit to an epoch right after the close approach of 13 April 2029, that changes the semi-major axis by about 0.18 au. The value of Y_M went from 22.45 × 10⁻⁴ au My⁻¹ obtained before the close approach, to 22.90 × 10⁻⁴ au My⁻¹ computed after the close approach. For the generic NEA we expect the change to be even smaller, because the distribution of the expected maximum semi-major axis drift is dominated by the uncertainties in the physical properties, rather than by the uncertainties in the orbital elements.

From the orbit determination point of view, an effect that is not taken into account in the model described in Sect. 2.1 is that of timing errors in the astrometric observations. It has been shown that the astrometric errors are affected by time calibration errors made by observers. These calibration errors can be either random or systematic, and they are especially important when an object has a fast apparent motion (Farnocchia et al. 2022). In turn, this may result in an erroneous estimation of the A₂ parameter, thus affecting the Yarkovsky effect detection. In this respect, the Aegis software used at the NEOCC is already able to treat timing errors by correcting the weights of the observations using the correlation with the right ascension and declination, with the method introduced in Farnocchia et al. (2022). However, the information about the timing errors is communicated only in the new Astrometry Data Exchange Standard (ADES, Chesley et al. 2017), which is not currently used to automatically compute asteroid orbits at the NEOCC (Faggioli et al. 2023). The switch of the Aegis system to automatically work with observations in ADES format is planned for the future, and the procedure described here will be able to automatically take into account also timing errors in the orbit determination.

Finally, the detection of a nongravitational force A₂ along the transverse direction does not necessarily mean that the Yarkovsky effect is responsible for it. In fact, other dynamical effects, such as outgassing or mass-shedding, may produce a transversal acceleration (Farnocchia et al. 2023). The usage of a statistical model for the physical properties of NEAs seems to help identifying these outliers, as confirmed by the cases of (523599) 2003 RM and 2006 RH₁₂₀ that were flagged as spurious. Still, there may be cases in which the acceleration produced by such phenomena is small enough to result compatible with the acceleration produced by the Yarkovsky effect. In these cases, the detection of such acceleration is real, and it is not possible to make a clear distinction between the Yarkovsky effect and activity events, unless there are proofs of such phenomena obtained from detailed studies of the single object.

4 Conclusions

In this paper, we introduced a procedure for the automated Yarkovsky effect detection on NEAs and for the identification of spurious detections. The procedure makes use of a statistical model to estimate the physical properties of an NEA, based on the most recent data and models about the NEA population. Physical properties are then used to estimate the expected maximum Yarkovsky semi-major axis drift. Candidates for detection are selected according to the current precision in the knowledge of the orbit, the observational arc, and the expected semi-major axis drift. Then, a 7-dimensional orbit determination of the candidates, that includes the estimation of a transversal nongravitational component, is performed. If the S/N of the A₂ detection is larger than 3, a statistical test to check the compatibility with the Yarkovsky effect physical model is carried out. Finally, the dependence of the detection on isolated tracklets is analyzed.

Among the known NEAs, we identified 348 positive Yarkovsky effect detections. This procedure is now adopted by the ESA NEO Coordination Centre, and data are automatically uploaded on the web portal and updated at every monthly orbital update issued by the MPC. A comparison with the results obtained by the JPL SSD is also shown, and the measurements are generally in a very good agreement, being based on similar orbit determination methods and datasets.

Acknowledgements

We thank Davide Farnocchia for the constructive comments on an early version of the manuscript that permitted us to improve the quality of the work.

Appendix A Densities of asteroid complexes

To determine values of the densities of asteroid complexes, we extracted data from the SsODNet database using the best value provided by the ssoCard. We selected only asteroids for which the taxonomic complex and the density were both determined, and we found a total of 386 objects. To determine the density distribution of a complex, we used the following procedure: 1) the density distribution of each asteroid is assumed to be log-normal; 2) a density sample for all the asteroids belonging to a given complex was randomly generated, obtaining a sample for the whole complex; 3) the sample for the whole complex was fitted with a log-normal distribution.

Table A.1 shows the list of asteroid complexes found in the SsODNet, together with the total number of density determinations, and the number of determinations with signal-to-noise ratio (S/N) larger than 3, 2, and 1.5. The sixth and seventh columns contain the value adopted in this work for each complex. The values for the C−, Ch−, M−, Q−, and S −complexes were fitted by using only the data with S/N > 3. The values for the L−, X−, and Z−complexes were obtained from the measurements with S/N > 1.5. The B−complex showed two outliers with relatively high density (i.e. (2) Pallas and (426) Hippo), that we removed from the sample used for the density fit. The E−complex also has a similar issue, because asteroids (44) Nysa and (3671) Dionysus show a density that is considered too small for these asteroid types, and therefore we removed them from the final fit. For the A−, P−, R−, K−, and U−complexes it was not possible to reasonably fit the sample with a log-normal distribution, therefore we used the population model mentioned in Sect. 2.2.2. Finally, the densities determined for asteroids in the V−complex show a large variety of values, therefore we adopted the value reported by Carry (2012). The log-normal fits obtained for the with the method described above for the B−, C−, Ch−, E−, L−, M−, Q−, S−, X−, and Z−complexes are all shown in Fig. A.1

Fig. A.1 Fit of the density distribution for the B−, C−, Ch−, E−, L−, M−, Q−, S−, X−, and Z−complexes. The histogram is obtained by generating log-normal distributions of asteroid density determinations of the corresponding complex, while the red curve is the best-fit log-normal distribution. Numerical values of the parameters µ and σ of the best-fit distribution are shown in the title of each panel.

Appendix B Additional tables

References

All Tables

All Figures

	Fig. 1 Distribution of the rotation period measurements of NEAs extracted from LCDB, together with the rotational disruption limit of 2.2 h and the dependence of P(D) on the diameter.
In the text

	Fig. 2 Distribution of the rms of the normalized residuals for the NEAs with a positive Yarkovsy effect detection. The blue (red) histogram represents the rms of the normalized residuals obtained with (without) the Yarkovsky effect included in the dynamical model. The dashed vertical lines correspond to the median of the two distributions, where the numerical value is reported in the legend.
In the text

	Fig. 3 Distribution of the S/N of the accepted Yarkovsky effect detections, obtained with equi-spaced bins of width equal to 1. NEAs with only optical observations are represented in the blue histogram, while the red histogram refers to NEAs observed also with radar.
In the text

Fig. 4 Distribution of the accepted Yarkovsky effect detections in the plane (D, |(da/dt)fit|). Green points are NEAs for which an estimate of the diameter is available, while black points are those for which the physical model is used. The error bars in the diameter are represented by using the 15th and 85th percentile of the diameter distribution (see Sect. 2.2.2), which roughly correspond to the 1-σ uncertainty in the case the diameter is estimated. The error bars in the semi-major axis drift represent the 1-σ uncertainty obtained by the conversion of A2, that has been estimated with orbit determination. The dashed blue line represents the fit with a function of the form α× Dβ. Red points are detections with S/N ≥ 3 that are considered spurious, hence rejected by our procedure.

In the text

	Fig. 5 Taxonomy distribution of the accepted Yarkovsky effect detections. Data about the taxonomic complex were retrieved from the SsODNet database Berthier et al. (2023).
In the text

	Fig. 6 Distribution of the common A2 measurements found in this work and at the JPL SBDB. The blue dashed line is the diagonal.
In the text

Fig. A.1 Fit of the density distribution for the B−, C−, Ch−, E−, L−, M−, Q−, S−, X−, and Z−complexes. The histogram is obtained by generating log-normal distributions of asteroid density determinations of the corresponding complex, while the red curve is the best-fit log-normal distribution. Numerical values of the parameters µ and σ of the best-fit distribution are shown in the title of each panel.

In the text