© The Authors 2024

1 Introduction

Photometric observations of asteroids can be used to determine their rotation state (in most cases, a unique rotation period and spin-axis direction) and shape. In contrast to the shape, which is most often only a convex approximation due to the limited information content of the disk-integrated light curves, the spin state can be determined very precisely. This precision increases with increasing number of available observations. The ability to accurately constrain the direction of the spin axis principally stems from a number of independent viewing geometries, defined by the observer (usually Earth), the asteroid, and the Sun, that are represented by the data. In the case of a few exceptional near-Earth asteroids, a sufficient number of observations may be achieved in one or two apparitions (Monteiro et al. 2020; Kwiatkowski et al. 2021). More often, only data accumulated over years or decades eventually help to determine the spin-axis direction to within a few degrees of accuracy. The precision of the rotation-period determination mainly depends on the time span covered by observations. With light curves that are observed over several decades and a typical rotation period of some hours, the rotation-period precision can reach about 0.01 s.

The simplest lowest-energy rotation model, characterized by a fixed direction of the spin axis in the inertial space and by a constant rotation frequency, is adequate to fit the observations of most asteroids. This model is about as sophisticated as the Keplerian description of the asteroid heliocentric motion on a fixed ellipse, however. As data spanning sufficiently long intervals of time become available and help determine the spin state very accurately, the basic model may need generalization by effects of even tiny torques, especially if they accumulate over time into a strong perturbation. The relevant torques originate (i) in the solar gravitational field and (ii) in the solar radiation.

The first torque (i) is a well-known effect that has been empirically known and was later theoretically analyzed by astronomers for centuries in the case of the rotation of Earth (the difference is only that two centers, the Sun and the Moon, act together to produce the resulting lunisolar precession of the Earth axis). The explanation of the effect was first given by Isaac Newton (e.g., Chandrasekhar 1995, Sect. 23) and was later mathematically mastered by Jean d’Alembert (D’Alembert 1749). While it is obvious for the Earth, it took centuries before the solar-induced gravitational precession was discussed in the context of asteroid rotation. The pioneering works in this respect came from the Uppsala school during the 1990s (e.g. Skoglöv et al. 1996; Skoglöv 1997, 1998; Skoglöv & Erikson 2002). Still, in spite of significant improvements in mathematical modeling, the effect was elusive to direct detection. Only when very precise measurements from the Dawn spacecraft visiting (4) Vesta became available was the asteroidal precession observationally determined (e.g., Konopliv et al. 2014). However, the accuracy was not impressive, to the point that some consider the Vesta case inconclusive (see, e.g., Archinal et al. 2018). The reason is that the effect is rather small for a roundish shape like that of Vesta, only ≃0.28° cy⁻¹, and Dawn stayed near Vesta for only about a year. The detection of asteroidal spin precession using much less precise Earth-bound observations can only be successful when the effect at a properly chosen target is much stronger than was found for Vesta. Ideally, this would require a slowly rotating, highly irregular body on a near-Earth orbit for which the effect may increase to a fraction of a degree in a year. Still, a good-quality data set spanning a long period of time is needed. The first hint that it might be possible to detect the asteroidal spin precession has recently been given by Ďurech et al. (2022b) for (1620) Geographos and (1685) Toro. However, the first convincing detection awaited the current paper. For the first time, we provide the detection of the spin-axis precession for (2100) Ra-Shalom. This object is nearly an optimum target for this goal: it resides on an Aten orbit, it rotates slowly, it has an irregular shape and obliquity away from the unfavorable values of 0°, 90° and 180^°, and the available numerous photometric data span a long interval of 44 yr. The details are given in Sect. 3.2.

The second torque (ii) is a novel process without relevance for planets, satellites, or even large asteroids. It has to do with radiation pressure that is imparted on small, irregularly shaped asteroids, and it consists of two components: a smaller torque due to directly reflected sunlight in the optical band, and a more important component due to the recoil of the thermally emitted radiation by the asteroid itself. The analysis of this phenomenon has undergone an impressive revival in the last two decades after the pioneering work of Rubincam (2000). Rubincam also coined the acronym YORP, which stands for Yarkovsky-O’Keefe-Radzievski-Paddack effect. Many of the important applications of YORP in planetary science can be found in the review by Vokrouhlický et al. (2015). While the gravitational torque only affects the direction of the spin axis, the YORP effect triggers secular perturbation of both the spin-axis direction and the spin frequency. Both are tiny effects, so that it is not surprising that the direct detection of YORP proved to be a tricky task. Vokrouhlický et al. (2004) pointed out that the only foreseeable way to detect the YORP effect is the secular change in the rotation frequency, which fortuitously produces a quadratic perturbation in the rotation phase. This is directly observable with photometric data. Taylor et al. (2007); Lowry et al. (2007); Kaasalainen et al. (2007) indeed followed this method and obtained the first YORP detections for the small near-Earth asteroids (54509) YORP and (1862) Apollo. Since then, the observations of a handful of other asteroids allowed us to detect the YORP torque directly, but the sample still remains very limited. Most importantly, the available results are still insufficient to solve the quantitative challenges posed by the theory of the YORP effect (see Vokrouhlický et al. 2015, for a review and detailed discussion). This also motivated our work: We aimed to increase the sample of asteroids with YORP detections, to determine their physical properties, and to compare the observed YORP with the theory.

Building partly on our previous research and also reporting on new targets, we present evidence for the YORP effect in four near-Earth asteroids: (i) An improved accuracy of the YORP effect detection for the two bodies (1862) Apollo and (161989) Cacus, and (ii) a determination of the YORP effect for the first time for another two bodies, (2100) Ra-Shalom and (138852) 2000 WN10. Several aspects related to these detections are notable. The YORP effect detection in (1862) Apollo now has the highest signal-to-noise ratio, with a very good prospect of further improvements. In the case of (138852) 2000 WN10, the YORP effect is detected for a prograde-rotating body for the first time. Finally, (2100) Ra-Shalom represents a difficult analysis in which the spin-axis precession is detected for the first time together with the YORP effect on a body with the lowest rotation frequency so far. Additionally, we analyzed available data for asteroid (85989) 1999 JD6, for which the YORP effect has previously been reported by Tian et al. (2022). However, the addition of new observations to the data set of this object shows that the YORP effect has not yet been detected at a statistically significant level. Nevertheless, we find it interesting that while a no-YORP solution is statistically admissible at a level of about one sigma for (85989) 1999 JD6, most solutions require the rotational frequency to be decelerated. If this is confirmed in the future, this would be a breakthrough case in the context of the YORP effect theory because all YORP detections so far have shown an acceleration of the rotation frequency (see Sect. 4).

The paper is organized as follows. In Sect. 2 we describe the available observation data sets, the reconstructed models, and the uncertainty in the rotation parameters obtained using our light-curve inversion method. We compare these values with the theoretical prediction in Sect. 3, and we discuss the broader context of our results in Sect. 4.

2 New photometric observations and model reconstruction

This section describes new photometric observations and model reconstruction for each asteroid in our sample. The final model is clearly based not only on our new data, but also uses previously published observations whenever available. The observatories contributing to our campaign are listed in Table A.1. Observations carried out with the Danish Telescope at La Silla (DK154) and the Ondřejov Observatory telescope (D65) were calibrated, the other photometry was relative without all-sky calibration. Observing times were converted from the reported UTC into TDB.

The new observations are summarized in the tables in the appendix. For simplicity, we only list individual nights, not individual light curves when more of them (due to observations in different filters or calibration issues) were observed in a single night with the same telescope.

To reconstruct the physical model, we used the light-curve inversion method of Kaasalainen et al. (2001) and Kaasalainen & Torppa (2001), in which the angular frequency ω evolves linearly in time as $\omega ={\omega }_0+\dot{\omega }t$, where the rate υ ≡ $\dot{\omega }$ = dω/dt and the initial value ω₀ are free parameters of the optimization (a comment on the validity of this approach can be found in Sect. 4). In the following text, we call υ a YORP parameter, although the connection between the change in the rotation rate detected from light curves and the YORP effect is discussed below in Sect. 3. Instead of ω₀, we report the sidereal rotation period P₀ = 2π/ω₀. Because the rotation period evolves with time, it is necessary to also report the epoch JD₀, for which P₀ is given. In our previous works (e.g., Ďurech et al. 2008, 2018, 2022b), we set JD₀ to the epoch of the first observation. However, this causes a correlation between υ and P₀, a shorter initial rotation period (at the beginning of the observing data set) and a smaller YORP produce about the same evolution in the rotation phase angle as a longer initial period that evolves faster due to a larger YORP parameter u. The rotation phase angle φ evolves over time t from some initial value φ₀ as $$\phi (t)={\phi }_0+\frac{2\pi }{P_0}\left(t-{\text{JD}}_0\right)+\frac{1}{2}\upsilon {\left(t-{\text{JD}}_0\right)}^2.$$(1)

With JD₀ at the beginning of the observing data set, the second term is always positive, so that by decreasing P₀, the phase φ increases for all observations, and this can partly be balanced by decreasing υ, from which a positive correlation between P₀ and υ arises. However, when JD₀ is somewhere in the middle of the observing interval, the second term is negative for half of the observations, and a small change in P₀ leads to an increase in φ for half of the observations and a decrease for the remaining half, which cannot be compensated for by changing υ. To avoid a correlation between υ and P₀, we now set JD₀ somewhere close to the center of the interval covered by observations.

To estimate the uncertainties in the spin and YORP parameters, we created 3000 bootstrap data sets and repeated the inversion. The bootstrap sample was a random selection of the same number of light curves; in each selected light curve, we bootstrapped its points. The parameter uncertainties reported below are the standard deviations of the parameter distributions.

We also included the precession effect of the spin axis due to solar torque for asteroids (2100) Ra-Shalom and (85989) 1999 JD6, for which this effect is strong enough to be taken into account (see Ďurech et al. 2022b, for details).

2.1 (1862) Apollo

Apollo was one of the two asteroids, together with (54509) YORP, for which the YORP effect was detected for the first time (Kaasalainen et al. 2007; Lowry et al. 2007; Taylor et al. 2007). Ďurech et al. (2008) published an update of the original detection and derived a YORP acceleration of (5.5 ± 1.2) × 10⁻⁸ rad d⁻² from 27 yr of observations from 1980–2007.

We observed Apollo during five more apparitions in 2014, 2017, 2019, 2021, and 2023. We also used other published observations from 2014 (Warner 2014) and 2021 (Warner & Stephens 2022) that are available in the ALCDEF database¹ (Warner et al. 2009). The full data set now consists of 87 light curves covering the time 1980–2023. The new light curves are listed in Table A.2. The observations taken with DK154 and D65 were absolutely calibrated in the Johnson-Cousins VR photometric system. We determined a color index (V – R) = 0.428 ± 0.010 mag, which is compatible with its Q-type spectral classification.

The new model we derived has the following spin parameters: an ecliptic longitude λ = (61.1 ± 4.6)°, an ecliptic latitude β = (−71.2 ± 2.2)°, a sidereal rotation period P₀ = (3.0654226 ± 0.0000002) h (given at JD₀ 2452276.0 epoch), and a secular change in the rotation rate υ = (4.94 ± 0.09) × 10⁻⁸ rad d⁻². The uncertainties of these spin parameters were estimated from the distribution of the bootstrap solutions shown in Fig. 1. The shape model is shown in Fig. 2. The difference in the rotation phase accumulated over 43 yr between a constant-period model with P₀ and our best-fit model is ~340°, almost one full rotation. In Fig. 3 we show the match between the data and model for the six new light curves as an example.

2.2 (2100) Ra-Shalom

Previous photometry data of Ra-Shalom were analyzed by Ďurech et al. (2012a) and later by Ďurech et al. (2018). In both cases, the data set was found to be too short to reveal the YORP effect, and the constant period model was found adequate to match the available observations. The best Ďurech et al. (2018) could find was that the 3σ uncertainty interval of the YORP parameter, namely −1.0 × 10⁻⁸ < υ < 1.5 × 10⁻⁸ rad d⁻², was slightly asymmetric with respect to the origin. This finding indicated a possible positive υ, which we aimed to detect here by adding new observations. The new photometric observations from 2019 and 2022 are listed in Table A.3. Part of these observations was conducted with the 28″ telescope at the WISE Observatory (Brosch et al. 2015). We also included data published by Warner & Stephens (2020a, 2023) in our analysis.

We also used light-curve observations from a citizen science project spearheaded by the SETI Institute, which connects telescope owners of the eVscope and eQuinox models developed by the France-based company, Unistellar. These digital, automated, and compact telescopes, boasting an 11.3 cm mirror, were tailored for people with very little background in astronomy. Consequently, the telescope is seamlessly operated using a smartphone through a dedicated mobile app.

An observing campaign focusing on Ra-Shalom commenced in August 2022. A myriad of citizen astronomers collated and submitted their data, covering a nearly two-month span. We selected a subset of data with adequate photometric quality, which is detailed in Table A.4. Nonetheless, these light curves were not incorporated into the final modeling because their photometric quality was inferior to that of our primary data set. This discrepancy can be attributed to the small aperture of the telescope and to the fact that the data typically encapsulated limited segments of the complete rotation phase (approximately one to two hours relative to the almost 20-h rotation period of Ra-Shalom). We still used this as a chance to corroborate the data reliability from the citizen project by juxtaposing it against our model predictions. This endeavor underscores the potential of the Unistellar burgeoning telescope network (comprising nearly 2000 potential contributors of various scientific modalities backed by Unistellar and its partner, the SETI Institute) for the photometry of pertinent targets, especially near-Earth asteroids (NEAs). For Ra-Shalom, a comprehensive data set from experienced professional and amateur astronomers was at our disposal, rendering the Unistellar observations superfluous for the model reconstruction. This burgeoning citizen astronomer network harbors significant potential for the future, however, particularly in capturing photometric data for any sufficiently luminous NEA (refer to the instance of (7335) 1989 JA, Lambert et al. 2023).

The first photometric observations of Ra-Shalom were taken in 1978. The interval covered by available observations therefore spans 44 yr, with 190 light curves in total. Being a near-Earth asteroid with the semimajor axis a = 0.832 au, an eccentricity e = 0.437, and a rotation period of about 19.82 h, the solar gravitational torque causes a potentially non-negligible precession of the Ra-Shalom rotation axis. In order to justify this extension of the model, along the lines developed in Ďurech et al. (2022b), we first theoretically estimated the corresponding perturbation of the rotation pole direction. The model without YORP has a pole direction (286.4°, −56.4°), a rotation period P = 19.820040 h, and the dynamical ellipticity of the shape model is ∆ = 0.232 (assuming a uniform density distribution). The formal uncertainty in the pole direction is only a few degrees. The precession constant a (see Eq. (1) in Ďurech et al. 2022b, for its definition) may be theoretically estimated with these data. We obtained α = 2430″ yr⁻¹, which corresponds to a shift of 27° on the precession cone, or 16° in λ and −4° in β over 44 yr. As a result, the pole precession is large enough and must be considered when interpreting the available light curves.

Following the approach of Ďurech et al. (2022b), we thus added the precession constant α, along with the YORP strength υ, as a second free parameter in the model. We scanned the values of υ and α on a grid to determine whether the fit level depends on these two parameters. Although the parameter α has to be positive, we formally tested values between −5000 and 5000″ yr⁻¹. The YORP parameter υ was let free to range the interval −5 and 5 × 10⁻⁹ rad d⁻². The result is shown in Fig. 4. The best solution with the lowest root mean square (RMS) residuals is for α = 3000″yr⁻¹ and υ = 2.9 × 10⁻⁹rad d⁻². We determined the uncertainty contours in the same way as in Vokrouhlický et al. (2011) or Polishook (2014), namely by defining an appropriate level of χ² with respect to the minimum χ².The 3σ interval defined in this way still covers the υ = 0 value, but the 1σ interval is about ±2 × 10⁻⁹rad d⁻², which is lower than the value of υ itself. Moreover, these χ² based error intervals are larger than the uncertainties estimated by bootstrapping (see Appendix A in Ďurech et al. 2022b). The bootstrap 1σ error of the YORP parameter for α = 3000″ yr⁻¹ is δυ = 8 × 10⁻¹⁰ rad d⁻². Thus, our conservative conclusion is that the YORP is detected at the 2σ level and the precession due to solar gravitation torque at the 1σ level.

The best model has a pole λ₀ = (257 ± 4)°, β₀ = (−52 ± 2)°, P₀ = (19.820072 ± 0.000008) h (both for JD₀ 2443763.0), and ∆ = 0.259. With this updated ∆ value, the theoretical precession constant α is now 2720″ yr⁻¹. This agrees excellently with the formally best value obtained from observations, given the uncertainty and simplification of the convex shape model. The final shape model is shown in Fig. 5, and its synthetic light curves are shown in Fig. 6.

Fig. 1 Bootstrap distribution of the pole direction in the ecliptic longitude λ and latitude β (top), the YORP parameter υ, and the rotation period P0 (bottom) for asteroid (1862) Apollo (the latter corresponding to epoch JD0 2452276.0). The color-coding corresponds to the density of points. The black cross marks the best-fit solution based on the full original (not bootstrapped) data set.

Fig. 2 Convex shape model of asteroid (1862) Apollo shown from the equatorial level (left and center, 90° apart) and pole-on (right).

Fig. 3 Example light curves of (1862) Apollo. The epoch of observation is shown by the label at the top of each panel. The blue points are the observed data, and the red curve is our best model including YORP. The geometry of the observation is described by the aspect angle θ, the solar aspect angle θ0, and the solar phase angle α.

Fig. 4 Color map of the reduced χ2 values of the light-curve fit for asteroid Ra-Shalom for different values of the precession rate α and the YORP strength υ. The red point marks the lowest χ2 value, the inner yellow contour is the estimated 1σ uncertainty limit, and the outer white contour is a 3σ boundary.

Fig. 5 Convex shape model of asteroid (2100) Ra-Shalom shown from the equatorial level (left and center, 90° apart) and pole-on (right).

2.3 (85989) 1999 JD6

Marshall et al. (2015) and later Marshall (2017) combined Goldstone and Arecibo radar observations obtained during the close approach in 2015 with optical light curves and reconstructed a physical model of this asteroid. The shape model of 1999 JD6 is a contact binary, with two connected lobes and dimensions 3.0 × 1.2 × 1.0km. The fortuitous geometry of the radar observations allowed Marshall (2017) to determine the pole direction very accurately, namely (λ,β) = (220.3°, −73.43°), with an uncertainty of 0.25°. The sidereal rotation period was 7.6643464 ± 0.0000056 h, and these authors provided an upper limit on the YORP value 1.6 × 10⁻⁶ deg d⁻², which corresponds to |υ| < 2.8 × 10⁻⁸rad d⁻².

Recently, Tian et al. (2022) analyzed the optical light curves spanning a 20-yr time interval (between 2000 and 2020). From this data set, they reported the YORP detection for this asteroid. In particular, their model has a rotation period 7.667749 ± 0.000009 h (for JD₀ 2451728.0) and a rotation pole direction λ = (232 ± 2)°, β (−59 ± 1)°.In addition, their analysis required a rotation acceleration υ = (2.4 ± 0.3) × 10⁻⁸ rad d⁻² (1σ error), which Tian et al. (2022) interpreted as the YORP effect.

To reconstruct our independent shape and spin model, we used a similar light-curve data set as Tian et al. (2022), but extended it by additional observations taken in May and June 2023. Altogether, our data extend over ten apparitions from 1999 to 2023 (see Table A.5). To our surprise, we obtained a different result than the previous two studies. Our best model is shown in Fig 7, has a pole direction (260^°, −60^°) and a period P = 7.664 354 h. When allowing for a nonzero υ, the best value is υ = −5.0 × 10⁻⁹ rad d⁻². The bootstrap estimates of the errors in the pole direction, period, and the YORP parameter are shown in Fig. 8. The 1σ uncertainties are 4.8° in λ, 4.0° in β, 1.5 × 10⁻⁶ h in P, and 4.1 × 10⁻⁹ rad d⁻² in υ. Therefore, even if the formally best fit and bootstrap samples suggest a negative YORP acceleration of about −5.0 × 10⁻⁹rad d⁻², its uncertainty is at about the same level, and zero YORP is still consistent with the current data set.

As to the difference with respect to the radar study, we note that our spin-axis direction does not agree with the precise value determined by Marshall (2017). Fixing the pole direction at their value still gives an acceptable fit, but the χ² value increases by about 5%. When we treat all data sets as relative light curves, the pole direction shifts to (252°, −67°), which is closer to the radar pole, but still outside the error intervals. The reason for this discrepancy is not known to us. We speculate that it might be related to our model being convex, while the real shape of 1999 JD6 is bilobed, as was clearly revealed by radar delay-Doppler images.

The convex shape model shown in Fig. 7 has a dynamical flattening of ∆ = 0.40. With a semimajor axis 0.883 au, an excentricity 0.633, and an inclination 17.06°, the theoretically estimated solar precession constant is α = 2130″ yr⁻¹. Over 24 yr of observations, this corresponds to a shift of almost 13^° on the precession cone, which is 7.3° in λ and −2.6° in β. These values are comparable to the pole uncertainties, so that the solar precession effect should also be included in the modeling because it has a non-negligible effect, and we wondered whether this omission might not contribute to the difference with respect to the results of Marshall (2017).

We proceeded similarly as in the case of asteroid Ra-Shalom: We probed values of α between −5000 and 5000″ yr⁻¹, and the YORP parameter value v was tested in between −3 and 3 × 10⁻⁸rad d⁻². The result is shown in Fig. 9. The formally best solutions with the lowest RMS residuals are obtained for α < 0, which is unrealistic. However, the estimated 3σ uncertainty interval is so wide that it includes all tested values of α. This means that our data set is not sensitive to α, and we currently cannot constrain it from observations. For any value of the precession constant, the constant period, namely υ = 0, is still consistent with the data (Fig. 10). However, as in the case of the bootstrap samples, the uncertainty interval is not symmetric around zero here either, and negative values of υ are preferred.

Fig. 6 Example light curves of (2100) Ra-Shalom. The blue points show observed data, and the red curve shows our best model with YORP and precession. The last two light curves are examples of Unistellar eVscope data that were not used for the modeling.

Fig. 7 Convex shape model of asteroid (85989) 1999 JD6 shown from the equatorial level (left and center, 90° apart) and pole-on (right).

Fig. 8 Bootstrap distribution of the pole direction in the ecliptic longitude λ and latitude β (top) and the YORP parameter υ and the rotation period P0 (bottom) for asteroid (85989) 1999 JD6. The color-coding corresponds to the density of points. The cross marks the solution based on the original (not bootstrapped) data set.

Fig. 9 Color map of the reduced χ2 values of the light-curve fit for asteroid (85989) 1999 JD6 for different values of the precession rate α and the YORP strength υ. The inner yellow contour is the estimated lσ uncertainty limit, and the outer white contour is a 3σ boundary.

2.4 (138852) 2000 WN10

The data set we used for the shape reconstruction consisted of light curves from 2008 and 2009 obtained by Skiff et al. (2012), from 2015 by Warner (2016) downloaded from the Asteroid Lightcurve Database² (Warner et al. 2009), and our new observations listed in Table A.6.

The light-curve inversion led to a strong YORP detection. The best spin solution has the following parameters: λ = (318 ± 4)°, β = (60 ± 8)°, P = 4.4636677 ± 0.0000006 h, (for JD₀ 2457023.5), and υ = (5.5 ± 0.7) × 10⁻⁸ rad d⁻². The corresponding best-fit shape model of (138852) 2000 WN10 is shown in Fig. 11. The bootstrap distribution of the poles and rotation parameters is shown in Fig 12.

In addition to the best-fit solution described above (Fig. 13), there was another local minimum in χ² for the sidereal rotation P = 4.465939 1 h. However, this solution was significantly worse, as it provided the same fit with nonzero YORP as the best period above without YORP. The pole direction was not unique, likely because of the limited geometry; there were also possible solutions with the pole within +10° around the direction β = 90°. Some of them are also shown in the upper plot in Fig. 12. However, for this spin-axis direction, a YORP value between 5 and 6 × 10⁻⁸ rad d⁻² also provides a significantly better fit than a constant period model.

2.5 (161989) Cacus

The previous model based on observations from 1978–2016 was published by Ďurech et al. (2018). Its spin parameters were λ = (254 ± 5)°, β = (−62 ± 2)°, and υ = (1.9 ± 0.3) × 10⁻⁸ rad d⁻². We updated the model by adding new light curves from 2022 that were observed with TRAPPIST and by Panfichi & Pajuelo (2023; see Table A.7).

The new model (Fig. 14) is essentially the same as the old one, with the same spin parameters, but smaller uncertainties: The pole direction (252 ± 3°, −63 ± 2°), period P = 3.7550527 (for JD₀ 2451544.5), and the YORP value υ = 1.86 × 10⁻⁸ rad d⁻².

However, in contrast to the old model by Ďurech et al. (2018), the YORP detection now does not depend on the isolated data from 1978. Even without the two light curves by Degewij et al. (1978) and Schuster et al. (1979), the best-fit YORP value is υ = 1.6 × 10⁻⁸ rad d⁻². The bootstrap results are shown in Fig. 15. The plot of the period versus YORP shows a core of points and diagonal tails. This pattern corresponds to the unusual time distribution of the photometric observations. Two light curves were taken in 1978, and the remaining part of the data set (another 33 light curves) consists of observations taken between 2003 and 2022. The outlying diagonal points correspond to bootstrap samples that by chance do not contain light curves from 1978. There were 381 such cases out of 3000. For them, the zeropoint lies even before the first observation, so that there is a strong correlation between the period and YORP. On the other hand, the bulk of the points in the center contains observations from 1978, so that the correlation is much weaker. The formal standard deviations for the whole bootstrap sample are 0.0000006 h and 1.3 × 10⁻⁹ rad d⁻². For a limited subsample including 1978 data, it is 0.0000003 h and 8.5 × 10⁻¹⁰ rad d⁻².

Fig. 10 Example light curves of (85989) 1999 JD6. The blue points show observed data, the red curve shows our best model without YORP, and the dashed curve shows the best model with the spin and YORP parameters fixed at values given by Tian et al. (2022; clearly not matching the available observations).

Fig. 11 Convex shape model of asteroid (138852) 2000 WN10 shown from the equatorial level (left and center, 90° apart) and pole-on (right).

3 Justification of the rotation rate change using a theoretical model

The secular change in the sidereal rotation rate, expressed by the parameter υ, has been empirically derived in Sect. 2, without a particular theoretical interpretation. This approach, namely including υ into the light-curve inversion method as a free, solved-for parameter, has been found to be very efficient in most previous studies. However, it ideally requires a subsequent step of connecting its value to a particular physical phenomenon. Spontaneous particle ejections (Scheeres et al. 2020) or micrometeoroid impacts (Wiegert 2015) may be examples of these candidate processes. We took a different standpoint here by adopting a link of the observationally determined υ values to the YORP effect (Vokrouhlický et al. 2015).

To set up this connection, we must evaluate υ_model using a numerical model. As always, the latter depends on many physical parameters p, thus υ_model(p). We may take the liberty of adjusting p to certain optimum values p_★ to match υ_model(p_★) ≃ υ, but the point is to ensure that p_★ has reasonable (expected) values. Unfortunately, the YORP is still slightly more complex: (i) Some of the parameters are well-constrained observationally or may be estimated using an educated guess (e.g., the rotation state, the large-scale convex shape, the size, the surface thermal inertia, and the bulk density, as an example of the latter group), but (ii) others are not easily accessible from ground-based observations (e.g., shape nonconvexities and small-scale surface irregularities). The predicted υ_model depends on (ii), however, which may represent a significant perturbation in specific cases (e.g., Statler 2009; Breiter et al. 2009). The specific physical phenomena that were considered to quantitatively test the YORP-dependence on these hidden parameters had to do with the lateral heat conduction in small-scale surface irregularities (e.g., Golubov & Krugly 2012; Ševeček et al. 2015; Golubov & Lipatova 2022), or with the anisotropy of the thermal emission (“thermal beaming”) and mutual irradiation of the surface facets (e.g., Rozitis & Green 2012, 2013a). We did not evaluate the contribution of (ii) here, and leave this part to future studies of specific targets. We only modeled (i). To do this we used our well-tested numerical approach that we presented in Čapek & Vokrouhlický (2004, 2005). In general terms, the model solves the 1D heat diffusion problem independently for each of the surface facets, with nonlinear Robin boundary conditions at the surface and deep interior below the facet. The time domain of the solution spans one revolution about the Sun, and the rotation period is slightly adjusted (typically within its uncertainty limits) to have an integer number of rotation cycles in one revolution period (in this way, the periodicity of the solution is imposed). The discretization steps in the space and time domains are chosen to satisfy the von Neumann criterion, and the iterations are repeated until a sub-degree tolerance for all surface facets is satisfied. After reaching convergence, we evaluated the total thermal torque and acceleration along the heliocentric revolution using the appropriate sum over all surface facets. More details are provided in the references mentioned above.

For the parameters p, we used the nominal (best-fit) rotation period, pole orientation, and shape model determined from the light-curve inversion in Sect. 2. Noting the strong dependence of the YORP effect on the asteroid shape, we chose an additional ten variant solutions of the spin state and shape, all providing a statistically acceptable fit to the light curves. We ran the model for all of them and report (i) the median and (ii) the range of the υ_model values. The other physical parameters, such as the size and bulk density, were either taken from previously published studies, constrained by the observed Yarkovsky effect, or assumed. Their values are individually discussed in the next sections. We note that υ_model does not depend on the surface thermal inertia in the simple model of the YORP effect we use here (see Čapek & Vokrouhlický 2004 for numerical evidence and Nesvorný & Vokrouhlický 2007; Breiter & Michalska 2008 for analytical proof). Therefore, we do not need to specify this parameter when we report the υ_model.

Even though our model caveat consists of the absence of the observationally hidden parameters, we further justify it by predicting the thermal acceleration in heliocentric motion (the Yarkovsky effect). All five NEAs analyzed in this paper have a reliably measured secular drift of the semimajor axis by now.

Similarly to the case of the steady increase in the sidereal rotation rate, we assume that the origin of this perturbation is due to the thermal accelerations in the heliocentric motion (the Yarkovsky effect). An important justification then arises from the common prediction of the YORP and Yarkovsky strength based on our model. Importantly, the Yarkovsky (orbital) part is somewhat less dependent on the details of the asteroid shape and surface roughness (our runs for the best-fit model and the sample of ten variant models usually provides a Yarkovsky effect prediction within 5–15% of the median value). As a result, it can provide a reasonable constraint on the bulk density. The value of the secular drift of the semimajor axis da/dt derived within the orbit determination process has been calculated in several independent studies published over the past decade or more (see, e.g., Chesley et al. 2008, 2016; Nugent et al. 2012; Farnocchia et al. 2013; Vokrouhlický et al. 2015; Tardioli et al. 2017; Del Vigna et al. 2018; Greenberg et al. 2020, listed here sequentially in time). As soon as the available astrometric information was accurate and abundant, the da/dt solution across these different sources converged to a single value. We used the most recently updated solution provided by the JPL Horizons website³. This is important because the statistically significant value of da/dt for (161989) Cacus only became available after the radar astrometry was taken in September 2022.

Fig. 12 Bootstrap distribution of the pole direction in the ecliptic longitude λ and latitude β (top), the YORP parameter υ, and the rotation period P0 (bottom) for asteroid (138852) 2000 WN10 (the latter at epoch JD0 2457022.5). The color-coding corresponds to the density of points, and the cross marks the solution based on the original (not bootstrapped) data set. Pole ecliptic latitudes β higher than 90° correspond to pole directions β − 90°, λ − 180°.

Fig. 13 Example light curves of (138852) 2000 WN10. The blue points show the observed data, the red curve shows our best model with YORP, and the dashed curve shows the best model without YORP.

Fig. 14 Convex shape model of asteroid (161989) Cacus shown from the equatorial level (left and center, 90° apart) and pole-on (right).

Fig. 15 Bootstrap distribution of the pole direction in the ecliptic longitude λ and latitude β (top), the YORP parameter υ, and the rotation period P0 (bottom) for asteroid (161989) Cacus (the P0 value is given at epoch JD0 2451544.5). The color-coding corresponds to the density of points, and the cross marks the solution based on the original (not bootstrapped) data set.

3.1 (1862) Apollo

While discovered already in April 1932 (see Wolf & Reinmuth 1932, becoming later a namesake of its orbital class among the near-Earth asteroids), this body shares the fate of several other NEAs by being lost for many decades and rediscovered only in 1973. A fortuitous pair of very close encounters with Earth in November 1980 and April 1982 offered the opportunity to take a sequence of multiwavelength observations: (i) visible photometry (Harris et al. 1987), (ii) infrared (Lebofsky et al. 1981), and (iii) even radar observations (Goldstein et al. 1981; Ostro et al. 2002). This wealth of data allowed an early determination of many physical parameters (although some needed slight corrections later on), making Apollo exclusive even in its category. The parameters included the rotation state and pole orientation (e.g., Harris et al. 1987), a size between 1.2 and 1.6 km, a geometric albedo of 0.2 to 0.26 (e.g., Harris 1998; Ostro et al. 2002), and a Q-type spectral classification (Tholen 1984; Bus & Binzel 2002). Apollo was recognized to be one of the first targets for a detection of the Yarkovsky effect by Vokrouhlický et al. (2005b; see even previous evidence of nongravitational perturbations of its heliocentric motion by Ziolkowski 1983 and Yeomans 1991, who suspected unseen cometary activity). The actual detection of the Yarkovsky effect was reported in 2008 (see Chesley et al. 2008), even preceded by the detection of the YORP effect (see Kaasalainen et al. 2007; Ďurech et al. 2008). Apollo has been found to have an unusually small satellite in its class (Ostro et al. 2005), which was uniquely detected in the radar observations taken in October and November 2005. Because it is so very small, not much is known about it, but it cannot influence our model determination of either the Yarkovsky or the YORP effects. Apollo was also taken as an exemplary case to show the possible role of planetary close encounters in explaining asteroidal Q-type spectra, which are similar to ordinary chondrite analogs (e.g., Nesvorný et al. 2010). In summary, not much is known about (1862) Apollo, but this body continues to be an inspiration for interesting concepts of asteroidal science.

One such step has been undertaken by Rozitis et al. (2013), who combined information about the detected Yarkovsky and YORP effects and interpreted them using a single thermophysical model of this asteroid (a goal that has been imagined by Čapek & Vokrouhlický 2004, see their Sect. 5, but not achieved before). We adopted their solution for an Apollo size D = 1.55 ± 0.07 km, a geometric albedo p_V = 0.20 ± 0.02, a surface thermal inertia $\Gamma ={140}_{-100}^{+140}\text{\hspace{0.17em}J\hspace{0.17em}}{\text{m}}^{-2}{\text{s}}^{-0.5}{\text{K}}^{-1}$, and a bulk density ${\rho }_{\text{b}}={2.85}_{-0.68}^{+0.48}\text{\hspace{0.17em}g\hspace{0.17em}}{\text{cm}}^{-3}$. In a way, our work is a mere repetition of the Rozitis et al. (2013) calculation, except that we now have updated and more accurate values of both YORP and Yarkovsky effects. As for the latter, the current astrometric data set provides da/dt = −(1.94 ± 0.16) × 10⁻⁴ au My⁻¹, while the value of Rozitis et al. (2013) was higher by about 10%. Similarly, our new value for the YORP effect (Sect. 2.1) is about 10% smaller and much more accurate than the value (5.5 ± 1.2) × 10⁻⁸ rad d⁻² from Ďurech et al. (2008) considered by Rozitis et al. (2013).

Figure 16 shows the comparison between the predicted Yarkovsky drift da/dt from our model and the observed value (the shape variants produce scatter of roughly ±10% about the median value used in Fig. 16 in addition to to the bulk density dependence). Assuming a bulk density of 2.75 g cm⁻³, which we nominally used in the simulation, the predicted semimajor axis drift matches the observed value very well. Small adjustments of size, surface roughness, and/or thermal inertia may further help to match any value in the uncertainty interval of da/dt. The same model predicts a median υ_model ≃ 5.0 × 10⁻⁸rad d⁻² for the 2.75 g cm⁻³ bulk density, and a range between 1.80 × 10⁻⁸ rad d⁻² and 7.25 × 10⁻⁸ rad d⁻². This is fairly consistent with the observed value υ ≃ 4.94 × 10⁻⁸rad d⁻² (Sect. 2.1). In addition to the global shape dependence, further differences may readily be explained by the thermal beaming of a rough surface (Rozitis & Green 2012), and/or simply a slightly larger body size.

3.2 (2100) Ra-Shalom

Ra-Shalom has received a wealth of infrared observations, both from the ground and from space (e.g., Harris et al. 1998; Delbó et al. 2003; Trilling et al. 2010; Usui et al. 2011). In spite of the persisting uncertainty, the results converge to a size of 2–2.8 km, a geometric albedo 0.1–0.18, and an unusually high value of the surface thermal inertia (possibly Γ ≃ 1000 J m⁻² s^−0.5K⁻¹). We adopted the results of Shepard et al. (2008), who combined various data sets, including radar observations, to obtain the most complete picture of this asteroid. Their preferred size was D = 2.3 ± 0.2 km and a geometric albedo p_V = 0.13 ± 0.03. Shepard et al. (2008) also analyzed a composite visible-to-infrared spectrum of Ra-Shalom and found that it was similar to spectra of CV3 meteorites. Based on this similarity, they argued that a K-type classification would be the best match for this asteroid (refining the previous C- or Xc-type classification⁴). Combining the measured grain density of CV3 meteorites and the estimated ~30% porosity of C-class asteroids, Shepard et al. (2008) suggested a bulk density of ρ_b = 2.4 ± 0.6 g cm⁻³. Finally, these authors also interpreted their radar measurements and thermal observations as providing evidence for a coarse or rocky surface with only a thin or unimportant regolith layer.

The large size and favorable orbit of Ra-Shalom, belonging to what Milani et al. (1989) classified as the Toro orbital group residing near or inside the mean motion resonances with the Earth, offered numerous opportunities to observe this asteroid with radar. It belongs to the record-holders in the number of radar observations at different apparitions (in this case, seven observations between 1981 and 2022). Their accuracy and long time-base resulted in a high-quality orbit determination that allowed a firm detection of the Yarkovsky effect (see already Vokrouhlický et al. 2005b). Interestingly, the related drift of the semimajor axis da/dt = −(2.01 ± 0.43) × 10⁻⁴ au My⁻¹ is well matched by our model, adopting a bulk density in between 2 and 3 g cm⁻³ and the suggested high surface thermal inertia (Fig. 17). These two parameters are correlated: If the thermal inertia were lower, the required bulk density would also be lower. Additionally, the unaccounted-for effects of the surface roughness in our model may slightly decrease the bulk density solution, as indicated by the arrow in the top panel of Fig. 17.

Ďurech et al. (2012b, 2018) repeatedly sought the YORP effect in the available photometric data of that date, but were only able to place an upper limit on the value. In this paper, we finally determined the value from a series of light-curve data extended to the last year. In contrast to the conclusion in Ďurech et al. (2018), who speculated that Ra-Shalom might be the first case for which YORP will be found to decelerate the rotation rate, we found a weakly accelerating signal of υ ≃ 2.9 × 10⁻⁹ rad d⁻² (Sect. 2.2). This is even lower than the value found for asteroid (1685) Toro, even though Toro (i) is larger and (ii) has a larger semimajor axis. This comparison confirms the unusual weakness of the YORP effect for Ra-Shalom.

Assuming a bulk density of 2.5 g cm⁻³, fine-tuned to match the Yarkovsky effect, our simulation predicts a median υ_model ≃ −2.1 × 10⁻⁸ rad d⁻². Not only is the absolute value of υ_model higher than the detected υ, but its sign is opposite (similarly to what has been found in Ďurech et al. 2018). Clearly, the nominal convex shape model, which just has a low resolution at large scales, cannot capture important features of the YORP effect. The full range of predicted υ_model values from the sample of 11 equivalent models is −0.15 × 10⁻⁸ rad d⁻² to −6.1 × 10⁻⁸ rad d⁻². This is still negative, but some of the extreme values are now closer to the detected υ.

A possibly important clue to the solution is the evidence of a significant surface roughness (Shepard et al. 2008). This could promote a significant contribution of the fine-scale surface effects, such as mutual facet irradiation and/or transverse heat conduction, of which at least the latter preferentially contributes by an acceleration component in YORP (e.g., Golubov & Krugly 2012). The situation of (2100) Ra-Shalom may therefore be reminiscent of (25143) Itokawa, for which the early analyses, making use of the large-scale resolution shape models of this asteroid, consistently predicted a negative υ value (e.g., Vokrouhlický et al. 2004; Scheeres et al. 2007; Breiter et al. 2009). However, after enough observations were finally available, the detected υ value proved to be positive (Lowry et al. 2014), which implied a spin-up state of Itokawa. Lowry et al. (2014), following an earlier theoretical concept formulated by Scheeres & Gaskell (2008; see also Breiter et al. 2009), suggested that the perturbing effect was due to a difference in the bulk density of the two main shape features, notably the quasi-ellipsoidal lobes in contact. However, the difference required to explain the overall effect appeared rather large. Ševeček et al. (2015) proposed an alternative, and perhaps more natural, explanation based on the novel concept of lateral heat conduction on small-scale surface features discovered by Golubov & Krugly (2012). We hypothesize that the same model may explain the disparity in the observed and modeled υ values for Ra-Shalom. However, we leave a detailed study for the future.

We now briefly return to the observationally revealed precession of the Ra-Shalom rotation pole (Sect. 2.2). We found that the best-fit shape model corresponds to the dynamical flattening value ∆ ≃ 0.26, which provides a theoretical value for the precession constant α ≃ 2720″ yr⁻¹ (which is fully compatible with the best-fit value of this parameter; Fig. 4). From the simplest possible point of view, in which the heliocentric orbit would be fixed in space, the pole would precess with an angular speed of −α cos ϵ ≃ 2522″ yr⁻¹, where ϵ ≃ 158° is the rotation pole obliquity, about the direction of the orbital angular momentum (therefore normal to the orbital plane; e.g., Bertotti et al. 2003). The situation becomes more complicated when the orbital plane evolves in time due to planetary perturbations. In the case of Ra-Shalom, the essence of the generalization may be described using a model that is still very simple, in which the orbital inclination remains approximately constant (I ≃ 15.7°), but the longitude of node precesses in inertial space about the sum of the planetary orbital angular momenta with a frequency s ≃ −22.98″ yr⁻¹. The dynamics of the rotation pole within this setup is described by the Colombo top model (e.g., Colombo 1966; Haponiak et al. 2020), whose implementation for asteroids can be found, for instance, in Vokrouhlický et al. (2006). We refer to these references for more details and note here that the Colombo top model has two fundamental parameters: (i) the orbit inclination I, and (ii) the ratio of frequencies κ = α/(2σ) ≃ −59. The complexity of the pole evolution is determined by the relation of κ to a critical parameter ${\kappa }_{\star }=-\frac{1}{2}{\left({\sin }^{2/3}I+{\cos }^{2/3}I\right)}^{3/2}\simeq -0.823$. The most complicated situation occurs when κ is just slightly smaller than κ_★. In our situation, κ is much smaller than κ_★, and this confines the complexity to a narrow obliquity slab near 90^°, which restores the expected regularity at high and low values of ϵ. The latter corresponds to the case of Ra-Shalom.

To verify the conclusion, we numerically integrated the evolution of the Ra-Shalom rotation pole over the 44 yr time span between the first and last available light-curve observations (September 1978 to September 2022)⁵. The results are summarized in Fig. 18. The 44 yr track of the pole evolution is shown by the solid black segment, and the gray lines show solutions of the Colombo top model (the proximity of the true evolution of the Ra-Shalom pole from our numerical simulation to the Colombo top solution starting at nearby initial conditions shows that the latter is a fairly good but not complete approximation).

The existence of the resonant zone at about 90° obliquity perturbs the Ra-Shalom obliquity only very slightly, producing a negligible variation during the time span of interest. The principal dynamical effect consists of a ⋍ 31° drift in longitudinal angle ϕ associated with the obliquity є (right panel in Fig. 18). The true angular distance δ between the Ra-Shalom pole direction in 1978 and in 2022 in space is smaller, however. We easily find that cosδ = 1 − 2sin² єsin²ϕ/2, and therefore, δ ⋍ 11.5°. Because the inclination of the Ra-Shalom orbit is rather small, most of this effect is projected into a change in ecliptic longitude of the pole, while a small part, about 4^°, is the change in ecliptic latitude.

Fig. 16 Predicted values of the Yarkovsky effect and density for asteroid (1862) Apollo. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (2.5, 2.75, and 3 g cm−3; see the labels and the middle curve for 2.75 g cm−3) are shown by solid curves. The abscissa is the surface thermal inertia in SI units (J m−2 s−0.5 K−1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.1, and an effective size of 1.55 km. The grayscale horizontal region shows the value −(1.94 ± 0.16) × 10−4 au My−1 from the orbit determination. The grayscale vertical region shows the range of the best-fit surface thermal inertia value ${140}_{-100}^{+140}\text{\hspace{0.17em}\hspace{0.17em}J\hspace{0.17em}}{\text{m}}^{-2}{\text{s}}^{-0.5}{\text{K}}^{-1}$ (see Rozitis & Green 2012). Top panel: Model-predicted bulk density to match the observed value of the semimajor axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. A nominal effective size of 1.55 km is used. If this value were higher or lower, the density solution would shift in the direction indicated by the arrows (preserving the ρb D value). The analysis of ten different shape models of (1862) Apollo reveals a variation of ±12% about the median value used in the figure. Effects of small-scale surface roughness, if important, would shift the solution toward a lower value by typically 10–30% (see Rozitis & Green 2012).

Fig. 17 Predicted values of the Yarkovsky effect and density for asteroid (2100) Ra-Shalom. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (2, 2.5, and 3 g cm−3; see the labels) are shown by solid curves. The abscissa is the surface thermal inertia in SI units (Jm−2 s−0.5 K−1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.2 and an effective size of 2.3 km. The grayscale horizontal region shows the value −(2.01 ± 0.43) × 10−4 auMy−1 from the orbit determination. The vertical dashed line indicates the suggested surface high-thermal inertia of ≃ 1000 Jm−2 s−0.5 K−1. Top panel: Model-predicted bulk density to match the observed value of the semimajor axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. A nominal effective size of 2.3 km is used. If this value were higher or lower, the density solution would shift in the direction indicated by the arrows (preserving the ρb D value). The analysis of 10 different shape models of (2100) Ra-Shalom reveals a variation of ±10% about the median value used in the figure. Effects of small-scale surface roughness, if important, would shift the solution toward a lower value by typically 10–30% (see Rozitis & Green 2012).

Fig. 18 Rotation pole evolution of (2100) Ra-Shalom due to the solar torque and its comparison with the simple Colombo top model. The coordinates on the axes are defined in the orbital plane, namely (i) the obliquity є (or cos є) on the right and left ordinates, and (ii) the longitude ϕ reckoned from the direction 90° away from the orbital node (located at ϕ = −90°; e.g., Breiter et al. 2005). The numerically propagated evolution of the Ra-Shalom pole over the 44 yr time span is shown by the solid black line. The gray lines are solutions of the simple Colombo top model. For most of the obliquity values, they are basically simple straight lines of approximately constant є. The exception is the Cassini resonant zone at є ⋍ 90° (with the separatrix shown by the red curves, and the stationary Cassini state 2 shown by the red dot). Unlike the case of asteroid (433) Eros, whose spin axis librates in the resonant zone (Vokrouhlický et al. 2005a), the evolution of the Ra-Shalom pole is simpler and only consists of nearly regular precession. The zoom in the right panel shows that the change in obliquity is minimum, but conforms to the Colombo top solution.

3.3 (85989) 1999 JD6

The most complex study of this asteroid was presented by Marshall (2017), who combined several sets of multiwavelength data from visible photometry and infrared observations to very detailed radar sensing during its 2015 apparition. The results revealed the highly elongated shape of a contact binary with a volume equivalent size of D = 1.45 ± 0.14 km and a surface thermal inertia Γ ⋍ 280 J m⁻² s^−0.5 K⁻¹ (see also Campins et al. 2009). The radar data were complemented with light-curve observations between 1999 and 2015, which allowed determining the sidereal rotation period to 7.6643464 ± 0.0000056 h. The abundant radar data in 2015 and their fortuitous viewing geometry discussed in Marshall (2017) resulted in an unusually accurate determination of the asteroid rotation pole, namely (λ, β) = (220.3°, −73.43°), with only a fraction of a degree of uncertainty. This implies an obliquity of nearly 180^°. Marshall (2017) also searched for the YORP signal in the data available until 2015, but he set only an upper limit |υ| < 2.6 × 10⁻⁸ rad d⁻². The recent analysis of the Yarkovsky effect with Gaia DR3 precise astrometric measurements included by Dziadura et al. (2023) provided da/dt = −(4.26 ± 0.98) × 10⁻⁴ auMy⁻¹ for this asteroid. The taxonomic classification of (85989) 1999 JD6 is not clearly constrained, with suggestions of K class (Thomas et al. 2014), X class (Carry et al. 2016), and L class (Binzel et al. 2019), which means that only some of the primitive classes are excluded. From BVR photometry, Polishook & Brosch (2008) concluded that (85989) 1999 JD6 belongs to the K or S class.

Recently, Tian et al. (2022) analyzed a larger set of photometric observations taken between 2000 and 2020 and claimed to have detected the YORP effect with υ = (2.4 ± 0.3) × 10⁻⁸radd⁻². However, the results we presented in Sect. 2.3 do not confirm their finding. Instead, we obtained υ = −(5 ± 4) × 10⁻⁹ radd⁻², with a zero value still statistically acceptable. This is a tighter constraint than was found by Marshall (2017), but it does not yet prove the YORP signal. The reason for the difference between our results and that in Tian et al. (2022) eludes us. The best we could have done is to double-check our solution, which is compatible for the YORP nondetection with that of Marshall (2017).

We also used our rotation state and convex shape model to estimate the expected value of the semimajor axis secular drift (the Yarkovsky effect) and the sidereal rotation rate secular change (the YORP effect). We assumed a 2 gcm⁻³ bulk density and ran simulations for a wide range of surface thermal inertia. The results are shown in Fig. 19. The predicted semimajor axis drift is generally consistent with the value determined by Dziadura et al. (2023), but if the low thermal inertia of the surface holds, the bulk density also needs to be low (≃1.3–1.7gcm⁻³). The intrinsic bulk density might be higher because by using the convex-hull model, we overestimated the volume that is occupied by matter for this contact binary object.

As expected, the formal median value υ_model ≃ 2.7 × 10⁻⁷ rad d⁻² obtained with our nominal shape model and a bulk density 2 g cm⁻³ is significantly different from the upper limit |υ| ≤ 10⁻⁸radd⁻² (with even a preference for the negative values). However, the full range of υ_model from our 11 shape/spin variant solutions extends from −3.5 × 10⁻⁷ radd⁻² to 3.7 × 10⁻⁷radd⁻², indicating that the YORP effect cannot be accurately predicted from the simple convex models and available photometric observations.

Fig. 19 Predicted values of the Yarkovsky effect and the density for asteroid (85989) 1999 JD6. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (1.5, 2, and 2.5 g cm−3; see the labels) shown by the solid curves; the abscissa is the surface thermal inertia in SI units (Jm−2 s−0.5 K−1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.3 and an effective size of 1.45 km. The grayscale horizontal region shows the value −(4.26 ± 0.98) × 10−4 au My−1 from the orbit determination. Top panel: Model-predicted bulk density to match the observed value of the semi-major axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. The nominal effective size of 1 km is used. If this value were lower or higher, the density solution would shift in the direction indicated by the arrows (preserving the pb D value). The analysis of ten different shape models of (85989) 1999 JD6 reveals a variation of ±10% about the median value used in the figure.

3.4 (138852) 2000 WN10

The population of the Earth coorbitals is an interesting transient subclass of NEAs with a characteristic dynamical lifetime of some tens of thousands of years (e.g., Morais & Morbidelli 2002). Coorbitals, or objects that are orbitally very close to them, may offer a unique possibility among NEAs to be observable every year for a decade or more. Because these annual Earth encounters may be close, even very small NEAs may be observed. In spite of the slight drawback that they typically repeat a similar observing geometry, coorbitals are particularly suitable in situations when a series of observations seeks to determine the cumulative effect of a weak perturbation. This is the case of the Yarkovsky and YORP effects. For this reason, Vokrouhlický et al. (2005b) considered the detection of the Yarkovsky effect for a handful of small coorbital asteroids, including (54509) YORP (for which the Yarkovsky effect was swiftly reported by Chesley et al. 2008) and (138852) 2000 WN10 (for which the Yarkovsky effect was reported by Vokrouhlický et al. 2015, due its larger size and unfortunate absence of the radar astrometry). Even more importantly, (54509) YORP was one of the first two objects for which the YORP effect was detected (Taylor et al. 2007; Lowry et al. 2007). The YORP detection for (138852) 2000 WN10 proved to be harder and had to wait until the present paper. This is due to its larger size and intrinsically weaker YORP strength⁶.

To compare the observed da/dt and υ (from Sect. 2.4) with the theoretical values, we need to know or assume several physical parameters of the asteroid. Because 2000 WN10 is small and lacks radar observations, the information about this asteroid is unfortunately sparse. The light-curve inversion provides complete information about the spin state and a rough (convex) shape model. Based on broadband photometry, Ieva et al. (2018) found (138852) 2000 WN10 to be compatible with the S-type group of asteroids. Perna et al. (2018) conducted a spectroscopic survey of small NEAs and included (138852) 2000 WN10 in their program. The visible spectrum corresponds to the Sq taxonomic type (similarly Binzel et al. 2019, found it to be a Q-type object). Adopting a mean geometric albedo of 0.24 for this group, they argued for a rather small size of ≃250 m. However, without infrared observations, the albedo and size values are simply an educated guess. If we were to assume the mean albedo value of the NEA population 0.15, the size would recalibrate to ≃320 m. At this moment, we adopt D ≃ 300 m for our thermal model. In practice, we rescale unconstrained dimensions of the shape model from Sect. 2.4 such that its volume is equivalent to a sphere with a diameter of 300 m. Finally, we need to adopt some value of (i) the surface thermal inertia Γ and (ii) the bulk density ρ_b. In this case, we ran a series of simulations for different Γ, sufficiently sampling the interval of values from thermophys-ical modeling of small NEAs (see, e.g., Delbó et al. 2015). For the bulk density, we used 2 g cm⁻³. The results for different ρ_b values were readily obtained with the da/dt ∝ ρ⁻¹ relation.

Figure 20 shows our theoretical predictions for the Yarkovsky semimajor axis drift as compared to the observed value da/dt = (17.8 ± 2.9) × 10^-4 auMy^-1 (the variation in the predicted values due to different shape models is small in this case, only ±3%). The plausible range 1.7 to 2.4gcm^-3 for the bulk density is fully acceptable, but is shifted to slightly lower values than those of the small Q-type asteroid (6489) Golevka, as also determined by the Yarkovsky effect (see Chesley et al. 2003). At the same time, the required thermal inertia in between 150 and 600 J m^-2 s^-0.5 K^-1 favorably compares to what is expected from limited cases of small S-type NEAs with good visible and infrared observations (Delbó et al. 2015). There is clearly room for further adjustments due to small changes in the effective size D and the effects of small-scale roughness (the arrows in the upper panel schematically illustrate their influence).

The same simulations that provided satisfactory results for the Yarkovsky effect were also used to compute the YORP effect. The light-curve data analysis resulted in a low and positive value υ ≃ 5.5 × 10^-8radd^-2, indicating that the rotation rate of 2000 WN10 increases slowly. However, using the obtained shape model, the bulk density of 2gcm^-3, and the effective size of 300 m, our simplified nominal model provides υ_model ≃ −8.4 × 10^-8 radd^-2. This is a similar mismatch as for Ra-Shalom, presented in Sect. 3.2. Here, however, we find a much stronger variation of the predicted umodel values computed for 11 possible solutions from light-curve fitting: −15 × 10^-7radd^-2 to 26 × 10^-7 radd^-2. Clearly, the sample of resulting convex shape models cannot predict the YORP effect accurately. This likely originates from two effects: (i) the very low obliquity of × 10^° means that the polar flattening of the shape model from the light-curve inversion is only poorly resolved, and (ii) the coorbital nature of the 2000 WN10 orbit inherently limits the viewing geometry for an Earth-bound observed for subsequent seasons (see discussion in Taylor et al. 2007). Because of the low obliquity, the solar aspect angle of our observations has been always limited to between 82 and 88^°. The shape-modeling procedure clearly interpolates the surface even in the zone that was poorly sampled by the observations. This produces a huge missing piece of information, however, that would be needed to determine υ_model. At this moment, we must satisfy ourselves by observing that υ and υ_model are generally on the same order of magnitude (or at least within the expected limits).

Fig. 20 Predicted values of the Yarkovsky effect and density for asteroid (138852) 2000 WN10. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (1.5, 2, and 2.5gcm−3; see the labels) are shown by the solid curves; the abscissa is the surface thermal inertia in SI units (Jm−2 s−05 K-1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.4 and an effective size of 300 m. The grayscale horizontal region shows the value (17.8 ± 2.9) × 10-4auMy-1 from the orbit determination. Top panel: Model-predicted bulk density to match the observed value of the semi-major axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. The nominal effective size of 300 m is used. If this value were higher or lower, the density solution would shift in the direction indicated by the arrows (preserving the ρb D value). The analysis of ten different shape models of (138852) 2000 WN10 reveals a variation of ±3% about the median value used in the figure. The effects of small-scale surface roughness, if important, would shift the solution toward a lower value by typically 10–30% (see Rozitis & Green 2012).

3.5 (161989) Cacus

This asteroid was discovered during its close encounter with the Earth in February 1978, but it was subsequently lost until its rediscovery in January 2003 (the close encounter in September 1997 was missed because the position of Cacus was at a low declination on the southern sky). Nevertheless, the 1978 apparition is important by providing an astrometric tie point, but mainly by providing the early photometric observations (Degewij et al. 1978; Schuster et al. 1979). Cacus has been regularly observed since 2003, and both astrometry and photometry were accumulated. The radar observations taken from the Goldstone DSS-14 dish in August and September 2022 have significantly increased the orbit accuracy. A low-accuracy detection of the YORP effect was reported by Ďurech et al. (2018) and was improved in this paper by adding photometric observations taken in 2022. The post-radar orbit also allows us to constrain the Yarkovsky effect fairly well, that is, the semimajor axis secular drift of da/dt = −(4.45 ± 0.78) × 10^-4 auMy^-1. Cacus was also fortuitously observed by the WISE spacecraft in its cryogenic phase (Mainzer et al. 2011), which allowed Ďurech et al. (2018) to determine its effective size D = 1.0 ± 0.2 km and the most likely range of the surface thermal inertia Γ = 500–800 J m^-2 s^-0.5 K^-1. Finally, Thomas et al. (2014) and Binzel et al. (2019) reported spectroscopic observations, classifying Cacus consistently as a Q- or S-class object. Overall, in a quite short interval of time, Cacus has been elevated to the rank of NEAs with fairly complete information.

We used the best-fit spin state and shape model found in Sect. 2.5, together with the assumed effective size of 1 km, to compute the Yarkovsky effect for a plausible range of surface thermal inertia and bulk density values. The results are shown in Fig. 21. The observed semimajor axis drift da/dt = −(4.45 ± 0.78) × 10^-4auMy^-1 is well compatible with a bulk density between 1.7 and 2.4gcm^-3 (Fig. 21). Because the taxonomic class is similar to (1862) Apollo and (138852) 2000 WN10, the inferred bulk density for Cacus appears to be acceptable (a slightly smaller size and/or lower thermal inertia may push the value higher and close to the value inferred for Apollo or Golevka).

As above, we were less successful when we tried to match the observed acceleration of the rotation rate υ ≃ 1.94 × 10^-8 rad d^-2 (Sect. 2.5). Using the nominal model as above for the Yarkovsky effect (and specifically, 2gcm^-3 bulk density), we obtained a median υ_model ≃ −2 × 10^-8 radd^-2 and a full range between −20 × 10^-8 rad d^-2 and 18 × 10^-8 rad d^-2 for a sample of ten variant shape or spin solutions. Similarly to the case of 2000 WN10, we suspect that the overall uncertainty in the shape model plays a fundamental role. The obliquity of Cacus is basically 180^°, implying that the spin axis stretch of the model, as well as other details of the shape model along this direction, are very loosely constrained. This means that the result strongly depends on model details, even at medium and large scales. At this moment, the data simply do not provide us with an accurate enough model to firmly predict the υ_model value.

Fig. 21 Predicted values of the Yarkovsky effect and density for asteroid (161989) Cacus. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (1.5, 2 and 2.5gcm-3; see the labels) are shown by the solid curves; the abscissa is the surface thermal inertia in SI units (J m-2 s-0.5 K-1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.5 and an effective size of 1 km. The grayscale horizontal region shows the value −(4.45 ± 0.78) × 10-4 auMy-1 from the orbit determination. The grayscale vertical region shows the range of the acceptable surface thermal inertia value 200 to 2000 Jm-2 s-0.5 K-1 and the best-fitting core 500 to 800 Jm-2 s-0.5 K-1 (see Ďurech et al. 2018). Top panel: Model-predicted bulk density to match the observed value of the semimajor axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. The nominal effective size of 1 km is used. If this value were higher or lower, the density solution would shift in the direction indicated by the arrows (preserving the ρb D value). The analysis of ten different shape models of (161989) Cacus reveals a variation of ±3% about the median value used in the figure. The effects of small-scale surface roughness, if important, would shift the solution toward a lower value by typically 10–30% (see Rozitis & Green 2012).

4 Discussion and conclusions

It has become tradition over the past years to conclude a paper that reported a new detection of the YORP effect by recalling the previous detections and noting that in all cases, the rotation rate increases (i.e., υ is positive for all asteroids). We continue this tradition. For brevity, we recall Table 3 of Tian et al. (2022), who listed the relevant last YORP detections. Our work represents two modifications: (i) an extension by two new detections for (2100) Ra-Shalom and (138852) 2000 WN10, but (ii) also one retraction, at least as far as the situation stands now, namely a nondetection of YORP for (85989) 1999 JD6. Taken together, the YORP effect has been detected for 12 asteroids, and the rotation rate has indeed been found to increase in all cases⁷. Based on our new results, we provide a list of YORP detections in Table 1. In addition to the spin parameters and the size, the table also lists the parameter C_Y, which is a nondimensional coefficient introduced by Rozitis & Green (2013b; see also Rossi et al. 2009) to describe the YORP strength for a particular asteroid. The C_Y parameter is the υ parameter normalized for the semimajor axis a, the eccentricity e, the size D, and the density ρ (we used the same density of 2500 kg cm^-3) according to the formula $$C_Y=v\frac{a^2\sqrt{1-e^2}\rho D^2}{G_1},$$(2)

where G₁ ≈ 6.4 × 10¹⁶ kgms^-2 is a modified constant of solar radiation (2/π times the solar radiation pressure per unit area at a unit distance 1 × 10¹⁷ kg m s^-2). The values of C_Y allowed us to directly compare the YORP strength after removing the effects of the size and heliocentric distance. Apollo has the highest value of C_Y = 0.025, which is also likely the reason for the good agreement of its theoretical value of υ_model with the observed value and why it is stable with respect to bootstrap shape variants. For the other three asteroids we modeled in this work, Ra-Shalom, 2000 WN10, and Cacus, C_Y is smaller by at least one order of magnitude. This means that the detected YORP is weak and thus is more sensitive to the shape details. The υ_model values are therefore spread much more widely for these three cases.

While not an exception from the rule, we may mention for the sake of interest that (138852) 2000 WN10 is the first case in which the YORP effect has been determined for a prograde-rotating asteroid. This is expected, however, because the theory does not expect any difference between the prograde-and retrograde-rotating cases. Retrograde rotators predominate in the NEA population and represent a fraction of about two-thirds of the total. This is confirmed by direct observations of their rotation pole (e.g., La Spina et al. 2004; Kryszczyńska et al. 2007), but it is also fairly well reflected in the statistics of the Yarkovsky detections (e.g., Farnocchia et al. 2013; Tardioli et al. 2017; Greenberg et al. 2020). In this respect, it seems a small portion that the YORP effect is detected for only one prograde-rotator in a sample of 12 objects (if YORP acts equally on prograde- and retrograde-rotating bodies). Using a binomial distribution, we find that the probability that this occurs by chance is ≃ 5%. As far as we know, no apriori selection bias exists to preferentially detect the YORP effect for retrograde-rotating NEAs. As a result, the situation may still indicate that the sample of NEAs for which the YORP effect was detected is still too small to draw bold conclusions.

Following the same line of argument, the binomial distribution may indicate the likelihood that none of the YORP detections revealed negative u if the model were not to give a preference to its sign. Simple algebra provides a probability of only 0.05% that 12 cases have the same sign of υ by chance alone. This would all become more realistic if the model were to give a preference to spinning up, rather than spinning down, by YORP. Golubov & Krugly (2012) discovered this hidden element in the theory by accounting for the lateral heat conduction in the small-scale surface features of the asteroid (an interesting attempt for a combined model of the YORP effect may be found in Golubov & Scheeres 2019). If the still uncertain fraction to accelerate the rotation rate is about 70–80%, the chance of missing a deceleration case out of 12 trials would increase to ≃ (1.5–7)%. Here again, we conclude that the sample of asteroids in which the YORP effect has been detected is small. This motivates further efforts to search for asteroids for which the YORP effect may be detected. Of particular interest may be a future follow-up of (85989) 1999 JD6. If the tendency from our solution in Sect. 2.3 is confirmed, it may become the first asteroid for which the detected YORP effect decelerates the rotation rate.

Our new detections add new insight into the properties of the sample of objects for which the YORP effect has been determined in several respects. For instance, υ = (2.9 ± 2.0) × 10^-9 rad d^-2 for (2100) Ra-Shalom is the lowest value detected so far. This confirms that it might be possible to extend YORP detections to large asteroids, eventually even beyond the near-Earth asteroid group if enough data are available in the future. Moreover, Ra-Shalom has the longest rotation period for which the YORP effect has been detected so far. At the same time, the Ra-Shalom YORP detection came together with the detection of the spin-axis precession. In this respect, we note that the detection of the YORP effect for an asteroid in a tumbling state is yet another novelty. An indication of such detection can be found in the analysis of the rotation state of the very small asteroid 2012 TC4 (see Lee et al. 2021). The accurately monitored rotation state of (99942) Apophis (e.g. Pravec et al. 2014; Lee et al. 2022; Ďurech et al. 2022a) might indicate that this is a very suitable candidate, especially if high-quality data were obtained from both space-borne and ground-based observations before its very close approach to the Earth in April 2029.

The accuracy of the observationally determined YORP strength for (1862) Apollo, expressed by a value υ/δυ ≃ 55 for the signal-to-noise ratio, is superior to all currently known cases (followed with not much worse solution for (1620) Geographos; see Ďurech et al. 2022b). Future observations, maybe within the next decade or two, will certainly continue to improve the solution. Additionally, for asteroids as large as (1862) Apollo, these observations do not require large instruments: 0.5–1 m scale telescopes may be suitable. As an example, if a good-quality light-curve observation of (1862) Apollo is obtained during its close approach in October 2030, the signal-to-noise ratio of the YORP detection would increase to ≃ 70–75. The question arises, however, whether this is even needed and what the scientific justification would be.

A straightforward answer is provided by continuing a high-accuracy theoretical modeling work along the line of Rozitis et al. (2013). If the theoretical prediction of the Yarkovsky and YORP strength could still be improved, spectacular constraints might be derived on parameters such as the bulk density (see Farnocchia et al. 2021, where only the Yarkovsky side was considered). In cases like (2100) Ra-Shalom (Sect. 2.2), (1620) Geographos, and/or (1685) Toro (Ďurech et al. 2022b), further observations are clearly motivated to better decorrelate the YORP effect (thus also to improve its accuracy) from the effect of the regular precession of the spin axis. The observationally determined value of the precession constant may also contribute significantly to the consistency of the whole approach and might even help to answer interesting scientific questions (e.g., improve the shape model or suggest a possible inhomo-geneity of the density distribution). One additional aspect of the YORP detection methods that have been used so far may be clarified in view of the prospects for further significant improvements in the detection accuracy.

The baseline model with the constant sidereal rotation frequency ω(t) = ω₀ has always been generalized in the simplest way, namely using the linear model ω(t) = ω₀ + υt so far (implying that ω₀ only holds at some arbitrary epoch t = 0). However, when the asteroid orbit is highly eccentric, all this is slightly more general because the YORP torque depends at any moment along the orbits on the instantaneous heliocentric distance, on the orientation of the body frame with respect to the heliocentric position vector, and on other factors (e.g., Rubincam 2000; Čapek & Vokrouhlický 2004). The empirical factor υ in the formula above is only the result of averaging the YORP effect over one revolution about the Sun. As a result, the correct model for the rotation frequency evolution in time should read ω(t) = ω₀ + υt + ∆ω(t), where ∆ω(t) is a periodic function with the baseline periodicity T of the asteroid revolution about the Sun. We may rearrange the formula to provide the variation in the rotation rate, so that it reads ω(t) − ω₀ = υ [t + T ζ (t)], where we introduced a nondimensional factor ζ(t) = ∆ω(t)/(T υ). Figure 22 shows ζ as a function of heliocentric distance rather than time, which is possible because the heliocentric motion on a fixed ellipse is periodic (here we clearly may neglect planetary perturbations). We note that the variation has a non-negligible amplitude of ⋍ 0.3. However, the principal point is that it does not accumulate in time. At maximum, its positive value lasts half of the orbital period T/2, and the variation in the rotation frequency is − υ ζ_max T/2, producing an accumulated phase shift of − υ ζ_max T²/4. For (1862) Apollo, this amounts to ⋍ 0.1°. By repeating the same calculation for (101955) Bennu, we found that the periodic part of the YORP signal only has an amplitude of ⋍ 0.006^° (because Bennu has a lower eccentricity and shorter orbital period). This is likely too small even for the very accurate observations of the OSIRIS-REx spacecraft mission.

We may conclude that in the case of large asteroids, such as Apollo, with currently reachable YORP detection, we may only rely on the secular change in the rotation rate υ, neglecting variations in the instantaneous radiation torques along the orbit. The situation would be reversed in a hypothetical experiment, when an extremely strong torque would affect the rotation of a very small body and the available observations would cover a timescale comparable to or even shorter than the revolution period T about the Sun. The future may perhaps offer this interesting possibility.

Fig. 22 Nondimensional factor ζ (on the ordinate) related to the periodic part of the rotation frequency change ∆ω (see the text), expressed as a function of the heliocentric distance (at the abscissa). The sense of motion is clarified by the arrows, which imply that the upper branch of the curve corresponds to the motion from perihelion to aphelion, and vice versa for the lower branch of the curve. This prediction corresponds to the specific case of (1862) Apollo; the secular part with the υ value determined in Sect. 2.1 has been subtracted from the rotation-rate ω(t) evolution and used for the normalization of the ζ function. The high eccentricity (⋍ 0.56) of the Apollo orbit implies that the amplitude of the effect is not too small. The red symbols indicate the epochs along the orbits for which light-curve observations are available. Some of the observations were taken during close encounters on extreme values of phase angle exceeding 90° at heliocentric distances smaller than 1 au (indicated by the vertical dashed line). However, Apollo is large enough to also be easily observable near the aphelion of its orbit at a heliocentric distance of 2.29 au (see also Table A.2).

Acknowledgements

This work has been supported by the Grant Agency of the Czech Republic, grants 20–04431S and 23–04946S. The work at Modra was supported by the Slovak Grant Agency for Science VEGA, Grant 1/0530/22. Y.K. thanks for supporting the European Federation of Academies of Sciences and Humanities (grants ALLEA EFDS-FL1-18). A.R. and C.S. acknowledge support from the UK Science and Technology Facilities Council. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 870403 (NEOROCKS). D.P. is thankful for the WISE Observatory staff.

Appendix A New photometric observations

The aspect tables below list the asteroid distance from the Sun r and from the Earth ∆, the solar phase angle α, the geocentric ecliptic coordinates of the asteroid (λ,β), and the observatory at which the data were taken (see Table A.1) or a reference to the original publication.

References

All Tables

All Figures

	Fig. 1 Bootstrap distribution of the pole direction in the ecliptic longitude λ and latitude β (top), the YORP parameter υ, and the rotation period P0 (bottom) for asteroid (1862) Apollo (the latter corresponding to epoch JD0 2452276.0). The color-coding corresponds to the density of points. The black cross marks the best-fit solution based on the full original (not bootstrapped) data set.
In the text

	Fig. 2 Convex shape model of asteroid (1862) Apollo shown from the equatorial level (left and center, 90° apart) and pole-on (right).
In the text

	Fig. 3 Example light curves of (1862) Apollo. The epoch of observation is shown by the label at the top of each panel. The blue points are the observed data, and the red curve is our best model including YORP. The geometry of the observation is described by the aspect angle θ, the solar aspect angle θ0, and the solar phase angle α.
In the text

	Fig. 4 Color map of the reduced χ2 values of the light-curve fit for asteroid Ra-Shalom for different values of the precession rate α and the YORP strength υ. The red point marks the lowest χ2 value, the inner yellow contour is the estimated 1σ uncertainty limit, and the outer white contour is a 3σ boundary.
In the text

	Fig. 5 Convex shape model of asteroid (2100) Ra-Shalom shown from the equatorial level (left and center, 90° apart) and pole-on (right).
In the text

	Fig. 6 Example light curves of (2100) Ra-Shalom. The blue points show observed data, and the red curve shows our best model with YORP and precession. The last two light curves are examples of Unistellar eVscope data that were not used for the modeling.
In the text

	Fig. 7 Convex shape model of asteroid (85989) 1999 JD6 shown from the equatorial level (left and center, 90° apart) and pole-on (right).
In the text

	Fig. 8 Bootstrap distribution of the pole direction in the ecliptic longitude λ and latitude β (top) and the YORP parameter υ and the rotation period P0 (bottom) for asteroid (85989) 1999 JD6. The color-coding corresponds to the density of points. The cross marks the solution based on the original (not bootstrapped) data set.
In the text

	Fig. 9 Color map of the reduced χ2 values of the light-curve fit for asteroid (85989) 1999 JD6 for different values of the precession rate α and the YORP strength υ. The inner yellow contour is the estimated lσ uncertainty limit, and the outer white contour is a 3σ boundary.
In the text

	Fig. 10 Example light curves of (85989) 1999 JD6. The blue points show observed data, the red curve shows our best model without YORP, and the dashed curve shows the best model with the spin and YORP parameters fixed at values given by Tian et al. (2022; clearly not matching the available observations).
In the text

	Fig. 11 Convex shape model of asteroid (138852) 2000 WN10 shown from the equatorial level (left and center, 90° apart) and pole-on (right).
In the text

Fig. 12 Bootstrap distribution of the pole direction in the ecliptic longitude λ and latitude β (top), the YORP parameter υ, and the rotation period P0 (bottom) for asteroid (138852) 2000 WN10 (the latter at epoch JD0 2457022.5). The color-coding corresponds to the density of points, and the cross marks the solution based on the original (not bootstrapped) data set. Pole ecliptic latitudes β higher than 90° correspond to pole directions β − 90°, λ − 180°.

In the text

	Fig. 13 Example light curves of (138852) 2000 WN10. The blue points show the observed data, the red curve shows our best model with YORP, and the dashed curve shows the best model without YORP.
In the text

	Fig. 14 Convex shape model of asteroid (161989) Cacus shown from the equatorial level (left and center, 90° apart) and pole-on (right).
In the text

	Fig. 15 Bootstrap distribution of the pole direction in the ecliptic longitude λ and latitude β (top), the YORP parameter υ, and the rotation period P0 (bottom) for asteroid (161989) Cacus (the P0 value is given at epoch JD0 2451544.5). The color-coding corresponds to the density of points, and the cross marks the solution based on the original (not bootstrapped) data set.
In the text

Fig. 16 Predicted values of the Yarkovsky effect and density for asteroid (1862) Apollo. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (2.5, 2.75, and 3 g cm−3; see the labels and the middle curve for 2.75 g cm−3) are shown by solid curves. The abscissa is the surface thermal inertia in SI units (J m−2 s−0.5 K−1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.1, and an effective size of 1.55 km. The grayscale horizontal region shows the value −(1.94 ± 0.16) × 10−4 au My−1 from the orbit determination. The grayscale vertical region shows the range of the best-fit surface thermal inertia value ${140}_{-100}^{+140}\text{\hspace{0.17em}\hspace{0.17em}J\hspace{0.17em}}{\text{m}}^{-2}{\text{s}}^{-0.5}{\text{K}}^{-1}$ (see Rozitis & Green 2012). Top panel: Model-predicted bulk density to match the observed value of the semimajor axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. A nominal effective size of 1.55 km is used. If this value were higher or lower, the density solution would shift in the direction indicated by the arrows (preserving the ρb D value). The analysis of ten different shape models of (1862) Apollo reveals a variation of ±12% about the median value used in the figure. Effects of small-scale surface roughness, if important, would shift the solution toward a lower value by typically 10–30% (see Rozitis & Green 2012).

In the text

Fig. 17 Predicted values of the Yarkovsky effect and density for asteroid (2100) Ra-Shalom. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (2, 2.5, and 3 g cm−3; see the labels) are shown by solid curves. The abscissa is the surface thermal inertia in SI units (Jm−2 s−0.5 K−1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.2 and an effective size of 2.3 km. The grayscale horizontal region shows the value −(2.01 ± 0.43) × 10−4 auMy−1 from the orbit determination. The vertical dashed line indicates the suggested surface high-thermal inertia of ≃ 1000 Jm−2 s−0.5 K−1. Top panel: Model-predicted bulk density to match the observed value of the semimajor axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. A nominal effective size of 2.3 km is used. If this value were higher or lower, the density solution would shift in the direction indicated by the arrows (preserving the ρb D value). The analysis of 10 different shape models of (2100) Ra-Shalom reveals a variation of ±10% about the median value used in the figure. Effects of small-scale surface roughness, if important, would shift the solution toward a lower value by typically 10–30% (see Rozitis & Green 2012).

In the text

Fig. 18 Rotation pole evolution of (2100) Ra-Shalom due to the solar torque and its comparison with the simple Colombo top model. The coordinates on the axes are defined in the orbital plane, namely (i) the obliquity є (or cos є) on the right and left ordinates, and (ii) the longitude ϕ reckoned from the direction 90° away from the orbital node (located at ϕ = −90°; e.g., Breiter et al. 2005). The numerically propagated evolution of the Ra-Shalom pole over the 44 yr time span is shown by the solid black line. The gray lines are solutions of the simple Colombo top model. For most of the obliquity values, they are basically simple straight lines of approximately constant є. The exception is the Cassini resonant zone at є ⋍ 90° (with the separatrix shown by the red curves, and the stationary Cassini state 2 shown by the red dot). Unlike the case of asteroid (433) Eros, whose spin axis librates in the resonant zone (Vokrouhlický et al. 2005a), the evolution of the Ra-Shalom pole is simpler and only consists of nearly regular precession. The zoom in the right panel shows that the change in obliquity is minimum, but conforms to the Colombo top solution.

In the text

Fig. 19 Predicted values of the Yarkovsky effect and the density for asteroid (85989) 1999 JD6. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (1.5, 2, and 2.5 g cm−3; see the labels) shown by the solid curves; the abscissa is the surface thermal inertia in SI units (Jm−2 s−0.5 K−1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.3 and an effective size of 1.45 km. The grayscale horizontal region shows the value −(4.26 ± 0.98) × 10−4 au My−1 from the orbit determination. Top panel: Model-predicted bulk density to match the observed value of the semi-major axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. The nominal effective size of 1 km is used. If this value were lower or higher, the density solution would shift in the direction indicated by the arrows (preserving the pb D value). The analysis of ten different shape models of (85989) 1999 JD6 reveals a variation of ±10% about the median value used in the figure.

In the text

Fig. 20 Predicted values of the Yarkovsky effect and density for asteroid (138852) 2000 WN10. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (1.5, 2, and 2.5gcm−3; see the labels) are shown by the solid curves; the abscissa is the surface thermal inertia in SI units (Jm−2 s−05 K-1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.4 and an effective size of 300 m. The grayscale horizontal region shows the value (17.8 ± 2.9) × 10-4auMy-1 from the orbit determination. Top panel: Model-predicted bulk density to match the observed value of the semi-major axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. The nominal effective size of 300 m is used. If this value were higher or lower, the density solution would shift in the direction indicated by the arrows (preserving the ρb D value). The analysis of ten different shape models of (138852) 2000 WN10 reveals a variation of ±3% about the median value used in the figure. The effects of small-scale surface roughness, if important, would shift the solution toward a lower value by typically 10–30% (see Rozitis & Green 2012).

In the text

Fig. 21 Predicted values of the Yarkovsky effect and density for asteroid (161989) Cacus. Bottom panel: Predicted semimajor axis drift da/dt (ordinate) due to the Yarkovsky effect from our model for three different values of the bulk density (1.5, 2 and 2.5gcm-3; see the labels) are shown by the solid curves; the abscissa is the surface thermal inertia in SI units (J m-2 s-0.5 K-1). We assumed the rotation state and shape model from the light-curve inversion in Sect. 2.5 and an effective size of 1 km. The grayscale horizontal region shows the value −(4.45 ± 0.78) × 10-4 auMy-1 from the orbit determination. The grayscale vertical region shows the range of the acceptable surface thermal inertia value 200 to 2000 Jm-2 s-0.5 K-1 and the best-fitting core 500 to 800 Jm-2 s-0.5 K-1 (see Ďurech et al. 2018). Top panel: Model-predicted bulk density to match the observed value of the semimajor axis drift shown by the grayscale region. The solid line in the middle shows the exact correspondence surrounded by a map of the sigma interval of the Yarkovsky drift. The nominal effective size of 1 km is used. If this value were higher or lower, the density solution would shift in the direction indicated by the arrows (preserving the ρb D value). The analysis of ten different shape models of (161989) Cacus reveals a variation of ±3% about the median value used in the figure. The effects of small-scale surface roughness, if important, would shift the solution toward a lower value by typically 10–30% (see Rozitis & Green 2012).

In the text

Fig. 22 Nondimensional factor ζ (on the ordinate) related to the periodic part of the rotation frequency change ∆ω (see the text), expressed as a function of the heliocentric distance (at the abscissa). The sense of motion is clarified by the arrows, which imply that the upper branch of the curve corresponds to the motion from perihelion to aphelion, and vice versa for the lower branch of the curve. This prediction corresponds to the specific case of (1862) Apollo; the secular part with the υ value determined in Sect. 2.1 has been subtracted from the rotation-rate ω(t) evolution and used for the normalization of the ζ function. The high eccentricity (⋍ 0.56) of the Apollo orbit implies that the amplitude of the effect is not too small. The red symbols indicate the epochs along the orbits for which light-curve observations are available. Some of the observations were taken during close encounters on extreme values of phase angle exceeding 90° at heliocentric distances smaller than 1 au (indicated by the vertical dashed line). However, Apollo is large enough to also be easily observable near the aphelion of its orbit at a heliocentric distance of 2.29 au (see also Table A.2).

In the text