An automated procedure for the detection of the Yarkovsky effect and results from the ESA NEO Coordination Centre

Context: The measurement of the Yarkovsky effect on near-Earth asteroids (NEAs) is common practice in orbit determination today, and the number of detections will increase with the developments of new and more accurate telescopic surveys. However, the process of finding new detections and identifying spurious ones is not yet automated, and it often relies on personal judgment. Aims: We aim to introduce a more automated procedure that can search for NEA candidates to measure the Yarkovsky effect, and that can identify spurious detections. Methods: The expected semi-major axis drift on an NEA caused by the Yarkovsky effect was computed with a Monte Carlo method on a statistical model of the physical parameters of the asteroid that relies on the most recent NEA population models and data. The expected drift was used to select candidates in which the Yarkovsky effect might be detected, according to the current knowledge of their orbit and the length of their observational arc. Then, a nongravitational acceleration along the transverse direction was estimated through orbit determination for each candidate. If the detected acceleration was statistically significant, we performed a statistical test to determine whether it was compatible with the Yarkovsky effect model. Finally, we determined the dependence on an isolated tracklet. Results: Among the known NEAs, our procedure automatically found 348 detections of the Yarkovsky effect that were accepted. The results are overall compatible with the predicted trend with the the inverse of the diameter, and the procedure appears to be efficient in identifying and rejecting spurious detections. This algorithm is now adopted by the ESA NEO Coordination Centre to periodically update the catalogue of NEAs with a measurable Yarkovsky effect, and the results are automatically posted on the web portal.


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. 2015), 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. 2015).Additionally, measurements of the Yarkovsky effect can be used to extrapolate the ⋆ Full tables B.1, B.2 and B.3 are only available in electronic form at the CDS via anonymous ftp to cdsarc.u-strasbg.fr(130.79.128.5) or via http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/.physical properties of the asteroids (see e.g.Chesley et al. 2014;Rozitis et al. 2014;Fenucci et al. 2021Fenucci et al. , 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 1 (NEOCC), the NEODyS 2 service, and the NASA JPL Solar System Dynamics 3 (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 sev-A&A proofs: manuscript no.47820corr eral 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. 2022) 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 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 Sec. 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 Sec. 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 Sec. 4.

Orbit determination
The dynamical model we used to determine the orbit of a 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 t of the form 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 where n is the mean motion, and p = a(1 − e 2 ) 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 where r = r/r is the radial unit vector, and r is the heliocentric position of the asteroid.The coefficient A 1 is computed as in Montenbruck & Gill (2000), where φ = 1.361 kW m −2 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 2 (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.4 under ESA contracts.

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 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 In Eq. ( 6), r = r/r is the heliocentric unit position vector, and ŝ is the unit vector of the asteroid spin axis.In addition, where S = πR 2 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 γ 1 , γ 2 are expressed as In Eq. ( 8), R ′ = R/l d is the scaled radius, where l d = K/(ρCω) 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 1 , k 2 , and k 3 are positive analytic functions of the scaled radius (Vokrouhlický 1998(Vokrouhlický , 1999)).In addition, Θ d = √ ρKCω/(εσT 3 ⋆ ), where σ is the Stefan-Boltzmann constant, and T ⋆ is the subsolar temperature, which is given by 4εσT 4 ⋆ = αF.
The seasonal component is expressed as where n is the unit vector normal to the orbital plane.The coefficients γ1 , γ2 have the same expressions as Eq. ( 8), but evaluated with Θ s = √ ρKCn/(εσT 3 ⋆ ), and with a scaled radius R ′ s = R/l s , where l s = K/(ρ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.

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 −1 K −1 , which is appropriate for NEAs (Delbo' et al. 2015).However, changing this value in the range 500 − 1200 J kg −1 K −1 , 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(Fenucci et al. , 2023)).Because we are interested in estimating the maximum semimajor 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 service5 (Berthier et al. 2023) through the Python rocks API6 .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 dependency 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 A&A proofs: manuscript no.47820corr with the diameter as 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 meters 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 meters, asteroids may rotate very fast without being disrupted, and therefore, we assumed a trend that scales as D −2 , which is a better fit to the data than the D −1 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. (2021Fenucci et al. ( , 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 We assumed Γ to be Gaussian distributed with a mean value of 250 J m −2 K −1 s −1/2 (Delbo' et al. 2007) and a standard deviation of 100 J m −2 K −1 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.10) Notes.Reference [1] above refers to Fenucci et al. (2021).

Treatment of poor astrometry
The fact that a high signal-to-noise ratio (S/N) detection of A 2 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 3) they were separated by at least five years 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 2 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 2 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 tracklets 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.

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 −5 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 drift7 (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 Sec.2.2.1, using 500 000 combinations of input parameters chosen from the parameter distributions described in Sec.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 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 is fulfilled, the asteroid is considered a good candidate for Yarkovsky effect detection, hence the value of A 2 is fitted through a 7-dimensional orbit determination (7-D OD).If the signal-to-noise ratio S/N = |A 2 /σ A 2 | 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 2 obtained from orbit determination could be actually caused by the Yarkovsky effect, we determine whether A 2 is statistically compatible with the physical properties of the given object or not.To this purpose, we convert A 2 and its uncertainty σ A 2 to semi-major axis drift values (da/dt) fit and σ (da/dt) fit respectively, by using Eq. ( 2).Then, we check if the value (da/dt) fit obtained from orbit determination fulfills the condition where k = 1.645 gives the 5-th percentile of the measured Yarkovsky drift.If inequality of Eq. ( 14) is fulfilled, it means that semi-major axis drift values down to |(da/dt) fit | − kσ (da/dt) fit 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 Sec.2.3, and keep only those detections passing the test.

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 Center8 (MPC).Optical observations are downloaded from the MPC, while radar observations are taken from the JPL Solar System Dynamics website9 .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   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 2 , for each positive detections.In this process, we assumed A 2 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 R/P = 2 +1 −0.7 (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 years, 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 R/P = 2 +0.8 −0.7 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 15-th and the 85-th per-centile of the diameter distribution, using the model described in Sec.2.2.2.Analytical models of the Yarkovsky effect predict that da/dt scales with the diameter as D −1 .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 −5 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.

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 4 , and 2015 TC 25 (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 1 , an NEA that is currently in the NEOCC Risk List 10 .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 2 t −1 (see Table 3 for the complete orbital solution), providing a new SRP detection that was not available in the literature 11 .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 120 (Seligman et al. 2023) and (523599) 2003 RM (Farnocchia et al. 2023).In addition to the radial component A 1 and the tangential component A 2 , it is possible to also fit a nongravitational out-of-plane force component A 3 for these two NEAs, showing a dynamical behavior similar to that of comets.In these two cases, the value of A 2 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. 10https://neo.ssa.esa.int/risk-list 11We note that this detection was however already present in the JPL SBDB since May 2021.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 years, and A 2 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 Sec. 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 2 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 2 obtained without including the isolated tracklet in the orbital fit, following the procedure described in Sec.2.3, and 30 of them still present a signal in A 2 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 124 and 2014 RO 17 , for which accurate measurements were performed by MM; 2015 LA 2 , 2018 KR, and 2016 UW 80 , 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 can not be maintained manually.Other 11 objects presented intraarc isolated tracklets, but the A 2 detection was affected only for 2005 TG 50 and 2011 XC, that were therefore considered as spurious.In addition, the S/N drops down also for 2009 OS 5 when observations from 2014 from Mauna Kea Observatory (IAU code 568) are excluded in the fit, however we considered them as reliable 12 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).

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. 2021Fenucci et al. , 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 13 and available through the HTTPS APIs 14 .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.

Comparison with the JPL SBDB
We downloaded the list of determinations of A 2 from the JPL Small Body Database 15 (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 1 , and 6 objects have a determination of both A 1 and A 3 .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 2 determinations, together with their 1-σ uncertainty.We fitted the data with a linear function of the form α|A 2 | + β with the ODR algorithm, that gave an estimation of the parameters of α = 1.00 ± 0.0001 and β = (−1.90± 0.97) × 10 −16 , 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 where the superscript JPL denotes the value provided by the JPL SBDB.For the vast majority of the asteroids in common, the relative errors ε 1 , ε 2 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 3 , 2003 AF 23 , (1685) Toro.However, also in these cases the values of ε 1 and ε 2 are only slightly larger than 3.We also performed a test by checking that 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 23 and (65717) 1993 BX 3 , 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 ε 1 and ε 2 , and of χ for all the NEAs with χ > 1.The full table in extended precision, used to generate Fig. 6, is available through the CDS.
Asteroids with a determination of A 2 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 2 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 2 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 2 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 2 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.

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 −4 au My −1 obtained before the close approach, to 22.90 × 10 −4 au My −1 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 Sec.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 2 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 2 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 120 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.

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 2 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.

Fig. 1 .
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.

Fig. 2 .
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.

Figure 3 Fig. 3 .
Figure3shows 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 9 , (99942) Apophis, (480883) 2001 YE 4 , (524522) 2002 VE 68 , (2340) Hathor, and 2002 BF 25 , which are not shown in Fig. 3.

Fig. 4 .
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 15-th and 85-th percentile of the diameter distribution (see Sec. 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 A 2 , 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.

Fig. 5 .
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).

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

Notes.
Columns are: Ast.= Designation of the asteroid; H = Absolute magnitude; RMS = RMS of normalized residuals for 7-D OD; RMS 6 = RMS of normalized residuals for 6-D OD; A 2 = Transversal nongravitational component obtained with 7-D OD; σ(A 2 ) = Error in A 2 ; da/dt = Semi-major axis drift obtained by converting A 2 ; σ(da/dt) = Error in the semi-major axis drift; Y M = 95-th percentile of the maximum expected Yarkovsky drift; S/N = Signal-to-noise ratio of the detection; Res.= Acceptance flag, can be only accepted (Acc.) or rejected (Rej.);#OO = Number of Optical Observations; #ReO = Number of rejected Optical Observation in the 7-D OD; #ReO 6 = Number of rejected Optical Observation in the 6-D OD; #RaO = Number of Radar Observations; #ReRa = Number of rejected Radar Observations in the 7-D OD; #ReRa 6 = Number of rejected Radar Observations in the 7-D OD; D 15 = 15-th percentile of the diameter distribution; D 50 = 50-th percentile of the diameter distribution; D 85 = 85-th percentile of the diameter distribution; P = Rotation period; Tax = Taxonomic class; ∆T = Length of the observational arc.The full list of results in extended precision is available at the CDS.The A 2 and σ(A 2 ) values are given in 10 −14 au d −2 , while the values of the semi-major axis drifts da/dt, σ(da/dt) and Y M are given in 10 −4 au My −1 .The diameter percentiles D 15 , D 50 , and D 85 are given in meters, the rotation period P in hours, and the arc length ∆T in years.

Notes.=
Columns are: Ast = Designation of the asteroid; A 2 = Transversal nongravitational component estimated by NEOCC; σ A 2 = Error in A 2 estimated by NEOCC; S/N = Signal-to-noise ratio of the NEOCC detection; A JPL 2 Transversal nongravitational component estimated by JPL; σ JPL A 2 = Error in A 2 estimated by JPL; S/N JPL = Signal-to-noise ratio of the JPL detection; JPL Model = Dynamical model used by JPL (1 = A 2 ; 2 = A 1 and A 2 ; 3 = A 1 , A 2 and A 3 ); ε 1 = Relative error w.r.t.NEOCC; ε 2 = Relative error w.r.t.JPL; χ = Chi value.The full list in extended precision is available at the CDS.The A 2 and σ(A 2 ) values are given in 10 −14 au d −2 .

Notes.
Columns are: Ast = Designation of the asteroid; da/dt = Semi-major axis drift computed by the NEOCC; σ(da/dt) = Uncertainty in the semi-major axis drift computed by the NEOCC; S/N = Signal-to-noise ratio of semi-major axis drift detection of NEOCC measurement; Y M = 95-th percentile of the expected maximum Yarkovsky semi-major axis drift; da/dt JPL = Semi-major axis drift computed by JPL; σ(da/dt JPL ) = Uncertainty in the semi-major axis drift computed by JPL; S/N JPL = Signal-to-noise ratio of semi-major axis drift detection of JPL measurement; JPL Model = Dynamical model used by JPL (1 = A 2 ; 2 = A 1 and A 2 ; 3 = A 1 , A 2 and A 3 ).The full list in extended precision is available at the CDS.The values of da/dt, σ(da/dt) and Y M are given in 10 −4 au My −1 .

Table 1 .
Gravitational parameters GM of the 16 most massive main-belt asteroids included in the dynamical model, and that of Pluto.
Brozović et al. (2015)l parameters of the 16 main-belt asteroids are extracted from the header of the JPL ephemerides DE441 (https://ssd.jpl.nasa.gov/ftp/eph/planets/ascii/de441/).The gravitational parameter of Pluto is taken fromBrozović et al. (2015).The values are reported in the units of the

Table 2 .
Physical parameters for the Yarkovsky effect model and the corresponding value or model.

Table B .
1. Results of the automated Yarkovsky effect determination procedure.

Table B .
3. Asteroids with measured Yarkovsky effect for the JPL SBDB, not identified or rejected by our procedure.