© The Authors 2024

1. Introduction

Binary asteroids are estimated to account for approximately 15% of the asteroid population (Margot et al. 2015; Virkki et al. 2022; Minker & Carry 2023; Liberato et al. 2024). The dynamical lifetime of near-Earth asteroids and the collisional lifetime of main-belt asteroids are typically comparable to or longer than 10 million years (Farinella et al. 1998; Gladman et al. 2000). Understanding the long-term dynamics of binary asteroids is crucial for unravelling their evolutionary paths to their current configurations and deciphering their history.

In classical theory, the tidal effect and the binary Yarkovsky-O’Keefe-Radzievskii-Paddack (BYORP) effect governs the long-term evolution of binary asteroid systems (Ćuk & Burns 2005; Ćuk 2007; Ćuk & Nesvorný 2010; Jacobson & Scheeres 2011; Jacobson et al. 2014). The tidal effect results from energy dissipation due to deformation, which drives the spin of the object towards a synchronous state and moves the orbit out of a synchronous state (Murray & Dermott 1999). A synchronous state occurs when the spin period equals the mutual orbital period. The BYORP effect is a net radiative force with a random direction and magnitude averaged over a period and is caused by anisotropic radiation from the asteroid’s irregular shape (Ćuk & Burns 2005).

However, some aspects are not satisfactorily explained by classical theory. The tidal synchronization timescale for small binaries is considerably long (e.g. ∼10 Myr) unless a small tidal parameter, Q/k, is assumed, where k is the Love number and Q is the quality factor (Murray & Dermott 1999; Ćuk & Burns 2005; Quillen et al. 2022; Zhou et al. 2024). The Yarkovsky-O’Keefe-Radzievskii-Paddack (YORP) effect predicts that half of secondaries are synchronous and the other half asynchronous; however, in reality the vast majority of known binaries (i.e.,∼90%) are singly synchronous where only the secondary is synchronous, like the Earth-Moon system (Zhou et al. 2024). Furthermore, four binary asteroids (namely the 1996 FG₃, 1999 KW₄, 2001 SL₉, and Didymos systems) have been detected with secular mutual orbital drifts. BYORP and tidal effects could potentially explain these drift rates, but refinements to shape models and tidal parameters are needed (Scheirich et al. 2015, 2021, 2024; Scheirich & Pravec 2022; Naidu et al. 2024; Cueva et al. 2024; Richardson et al. 2024). This suggests that additional factors may be influencing the long-term dynamics of binary asteroids.

Recently, Zhou et al. (2024) revisited the Yarkovsky effect on binary asteroid systems, building on pioneering work related to Earth’s satellites and planetary rings (Rubincam 1982, 2006, 2014; Milani et al. 1987; Farinella & Vokrouhlický 1996; Métris et al. 1997; Vokrouhlický et al. 2005, 2007). The binary Yarkovsky effect arises from eclipse-induced thermal perturbations and thermal radiation from the other component in the system. By reorganizing the solution provided by Vokrouhlický et al. (2007), Zhou et al. (2024) found that the binary Yarkovsky effect on a secondary asteroid with a low inclination can drive the mutual orbit towards a synchronous state. The typical timescale for this synchronization process is approximately 0.1 Myr, which can be shorter than the timescales associated with tidal effects and the YORP effect. Thus, it has been proposed that the binary Yarkovsky effect may account for the synchronous state observed in the majority of known binary asteroids (Zhou et al. 2024). The binary Yarkovsky effect ceases to operate on the secondary once it enters a synchronous state.

However, there should also be a mirror binary Yarkovsky effect on the primary asteroid (hereafter the BYP) similar to the binary Yarkovsky effect on the secondary asteroid (BYS). The BYP is expected to continue modifying the mutual orbit after the BYS is inactive for the synchronized secondary, suggesting that most observed binary asteroids (∼90%) should be influenced by the BYP in addition to tidal effects and the BYORP effect. Preliminary estimates indicate that the orbital drift rate caused by the BYP is weaker than that of the BYS, by roughly a factor of (r_s/r_p)², where r_s is the radius of the secondary and r_p is the radius of the primary. However, because the primary experiences only a partial eclipse due to the smaller size of the secondary, unknown non-linear effects may be at play.

This study explores the behaviours of the BYP and its influence on singly synchronous binary asteroids and evaluates its typical strength using numerical methods. This paper is organized as follows: Section 2 overviews the mechanism and main equations of the binary Yarkovsky effect; Sect. 3 describes the numerical method used in this work to investigate the BYP; Sect. 4 discusses the numerical result, develops the empirical approximate formula for the BYP, and presents the estimated orbit drift rates of known synchronous binaries due to the BYP.

2. Analytical consideration

In this paper, subscripts ‘s’ and ‘p’ denote parameters related to the secondary and primary, respectively. According to Zhou et al. (2024), for binary asteroids with aligned spin vectors, mutual orbital vectors, and heliocentric orbits, the BYS-induced orbital drift rate is

$$\begin{array}{c}\hfill {\dot{a}}_{\mathrm{BY},\mathrm{s}}=\frac{2f_{\mathrm{BY},\mathrm{s}}{\mathcal{F}}_{\mathrm{s}}}{n},\end{array}$$(1)

with

$$\begin{array}{c}\hfill f_{\mathrm{BY},\mathrm{s}}=f_{\mathrm{YS},\mathrm{s}}(1-\frac{\pi r_{\mathrm{p}}}{4a}).\end{array}$$(2)

Here f_BY and f_YS are the binary Yarkovsky coefficient and the Yarkovsky-Schach (YS) coefficient, respectively. The n is the mean motion, and the a is the semi-major axis of the mutual orbit. The nominal radiation pressure per unit mass, ℱ, is defined as

$$\begin{array}{c}\hfill \mathcal{F}=\frac{\mathrm{\Phi }(1-A)\pi r^2}{\mathit{mc}}·\end{array}$$(3)

The coefficient f_YS is a complicated function of the physical properties of the components. However, in the limit that the size of the binary components is larger than a few centimetres, which is true in the context of binary asteroids, f_YS can be approximated as

$$\begin{array}{c}\hfill f_{\mathrm{YS},\mathrm{s}}=\frac{4r_{\mathrm{p}}}{9\pi a}\frac{{\mathrm{\Theta }}_{\mathrm{s}}}{2+2{\mathrm{\Theta }}_{\mathrm{s}}+{\mathrm{\Theta }}_{\mathrm{s}}^2}\mathrm{Sign}(n-{\omega }_{\mathrm{s}}).\end{array}$$(4)

Here ω is the spin rate, and r_p/πa accounts for the time fraction of the eclipse. The thermal parameter, Θ_s, is defined as

$$\begin{array}{c}\hfill {\mathrm{\Theta }}_{\mathrm{s}}=\frac{\mathrm{\Gamma }\sqrt{|\mathrm{\Delta }|}}{\epsilon \sigma T_{\mathrm{sub}}^3}·\end{array}$$(5)

Here Γ is the thermal inertia, Δ_s = ω_s − n is the relative frequency, ϵ is the emissivity, σ = 5.67 × 10⁻⁸ W m⁻² k⁻⁴ is the Stefan-Boltzmann constant, and T_sub = ((1 − A)Φ/ϵσ)^1/4 is the subsolar temperature.

The total binary Yarkovsky effect is the combined contribution of both the primary and secondary asteroids. For the primary component, in principle, there is a mirror effect to the secondary. The formula for the primary can be obtained by simply replacing the subscript ‘s’ with ‘p’ in Eqs. (1) and (4), which gives

$$\begin{array}{c}\hfill f_{\mathrm{YS},\mathrm{p}}=\frac{4r_{\mathrm{s}}}{9\pi a}\frac{{\mathrm{\Theta }}_{\mathrm{p}}}{2+2{\mathrm{\Theta }}_{\mathrm{p}}+{\mathrm{\Theta }}_{\mathrm{p}}^2}\mathrm{Sign}(n-{\omega }_{\mathrm{p}}).\end{array}$$(6)

However, due to the complexity of the geometry of the shadow cast on the primary component by the secondary component, this may not be accurate. In the derivation of Eq. (4), we assumed that the entire body stops receiving sunlight once its mass centre enters the shadow region, for the sake of simplifying the calculation. This assumption is reasonable for the secondary but not for the primary component since only a portion of the primary enters the eclipse at any given time (similar to a solar eclipse on Earth caused by the Moon). This partial eclipse reduces the strength of the YS effect, which heavily depends on the shadow; because of its induced non-linearity, a numerical approach is required to solve it.

This study ignores the planetary Yarkovsky (pY) effect on the primary, which is caused by the radiation from the secondary asteroid, as it is expected to be minimal. It would only induce a correction factor of (1 − πr_s/4a)∼1 for the overall effect for the primary.

3. Numerical model

The radiation force of a surface element is ∝σT⁴S, where T denotes the temperature and S denotes the surface area. Thus, to obtain the radiation force, we need to obtain the surface temperature of the shape model.

An irregular shape could produce the BYORP torque (Ćuk & Burns 2005). To eliminate the influence of the BYORP effect, we chose a spherical shape model approximated by a polyhedron with 1280 triangulated facets as the shape model. This symmetric shape model introduces a BYORP-induced error of only 10⁻⁶ to f_YS, which is negligible compared to f_YS ∼ 10⁻³.

The temperature, T, of the surface and the layer beneath is governed by

$$\begin{array}{c}\hfill \frac{\partial T}{\partial t}=\frac{\kappa }{C\rho }\frac{{\partial }^{2}T}{\partial z^2},\end{array}$$(7)

with two boundary conditions,

$$\begin{array}{cc}\hfill \kappa \frac{\partial T}{\partial z}|{}_{z=0}& =E(t)-e\sigma T^4{|}_{z=0},\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill \kappa \frac{\partial T}{\partial z}|{}_{z\to \infty }& =0,\hfill \end{array}$$(9)

where t is the time, κ is the thermal conductivity, C is the specific heat capacity, ρ is the bulk density of the asteroid, and e is the emissivity. In this study we set κ = 0.1 W m⁻¹ K⁻¹, C = 550 J K⁻¹ kg⁻¹, and ρ = 2000 kg m⁻³ for both the primary and secondary components by default. This gives a thermal inertia $\mathrm{\Gamma }=\sqrt{C\rho \kappa }\sim 330\mathrm{tiu}$ [J m⁻² K⁻¹ s^−1/2], which is close to the measured mean value of ∼200 tiu for kilometre-sized near-Earth objects (NEOs; Delbo’ et al. 2007). According to in situ measurements, asteroids Ryugu and Bennu have thermal inertia values of 225 ± 45 tiu (Shimaki et al. 2020) and 350 ± 20 tiu (DellaGiustina et al. 2019), respectively. The one-dimensional heat conduction equation (Eq. (7)) is sufficient for our problem since our interests lie in binary asteroid systems whose components are much larger than the thermal penetration depth; the latter is approximately a few centimetres. A full three-dimensional heat conduction equation is needed when investigating small (sub-centimetre-scale) dust.

We solved the temperature T_i, j, k at a depth of (j − 1)δz below the i-th facet at the k-th time step using the numerical scheme described in Zhou & Michel (2024). The thermal equilibrium state, where the intake energy equals the released energy of the system, is usually established after ∼50 mutual orbits unless the thermal inertia is considerably high (e.g. > 1000 tiu). A snapshot of our simulation is shown in Fig. 1. With the surface temperature, one can obtain the radiation force in the k-th timestep for the whole body via

$$\begin{array}{c}\hfill {\mathit{F}}_k=-\frac{2ϵ\sigma }{3c}\sum _i{T}_{i,0,k}^4{S}_i{\mathit{n}}_i,\end{array}$$(10)

Fig. 1. Left: Snapshots of the temperature field of the binary asteroids. Diagrams (a)–(d) illustrate the anti-clockwise orbit of the secondary asteroid around the primary. Both the primary and secondary have spin rates of 3 hours. The other properties of the binary system are detailed in Sect. 4.1. In phase (a) the primary is partially eclipsed by the secondary, while in phase (c) the primary fully eclipses the secondary. Right: Tangential accelerations due to thermal forces for the primary (blue) and secondary (red). The eclipse periods are represented by shaded areas.

where S_i and n_i are the area and the normal vector of the i-th surface element, respectively. The tangential component of the force, which is the source of the Yarkovsky effect, in the k-th timestep is then easily obtained as F_k · v_k/v_k, assuming a circular orbit. The YS coefficient is calculated as

$$\begin{array}{c}\hfill f_{\mathrm{YS}}=\frac{1}{k_{\max }}\sum _k\frac{{\mathit{F}}_k·{\mathit{v}}_k·c}{v_k\mathrm{\Phi }(1-A){\sum }_i{S}_i},\end{array}$$(11)

with A, the albedo of the asteroid, being ∼0.1. The mutual heating between the primary and the secondary was ignored due to its minimal contribution to the overall effect on the primary.

4. Results and implications

We investigate the YS coefficient as a function of the frequency ratio, ω_p/n, in Sect. 4.1 and propose an empirical modified formula for BYP in Sect. 4.2. Finally, we discuss the implications of the BYP on the long-term evolution of binary asteroids and estimate its magnitude for known synchronous binary asteroids in Sect. 4.3.

For simplicity, we assumed the spin vector of the primary asteroid aligns with the mutual orbit vector and considered the mutual orbit to be co-planar with the heliocentric orbit, given the low inclinations (i.e. i < 30°) observed in real binary asteroids (Pravec et al. 2012). Low inclinations are likely common as these systems are believed to form via YORP-driven spin-up and mass shedding followed by re-accumulation in the equatorial plane of the primary (Walsh et al. 2008; Pravec et al. 2012; Zhang et al. 2021; Agrusa et al. 2024). We note that this assumption leads to an upper limit for our estimated Yarkovsky effect. The dependence of the BYS and BYP on non-zero inclinations is discussed in Appendix A.

4.1. Synchronization of the mutual orbit

We considered a typical small binary asteroid system on a heliocentric circular orbit with a heliocentric semi-major axis a_h = 1 au. We set the semi-major axis of the mutual a ≃ 3.1r_p to keep the orbital period P = 12 h. By varying the spin periods of the primary and the secondary, we obtain the YS coefficient as a function of the frequency ratio, ω/n, which can be easily translated to Δ_s via Eq. (4) and Δ_p via Eq. (6).

Figure 2 shows the YS coefficient, f_YS, p, for the primary, obtained from the numerical model and the simple analytical model, as a function of the frequency ratio, ω_p/n. The Yarkovsky force has the same sign as the YS coefficient. The YS coefficient is positive when ω_p/n is smaller than 1, indicating the mutual orbit will expand under the Yarkovsky force. On the other hand, when ω_p/n is larger than 1, the YS coefficient is negative, implying a shrinking orbit. The Yarkovsky force vanishes when the primary is in a synchronous state, ω_p = n. Therefore, we conclude that the BYP pushed the mutual orbit towards the synchronous orbit of the primary component.

Fig. 2. YS coefficient for the primary as a function of the ratio of spin velocity to mean motion. The primary asteroid has a radius of 1 km, and the secondary has a radius of 0.3 km. Gold circles represent the numerical results. The solid gold curve shows the simple analytical solution from Eq. (6), while the dashed gold curve shows the modified solution from Eq. (12) with α = 1.7. The ωp/n distribution of confirmed binary asteroids is depicted in the blue histogram, which is truncated at 10.

However, we observe that the magnitude of the numerical results deviates slightly from the analytical results, which may be due to the partial eclipse on the primary, as discussed in the following section.

4.2. Empirical formula for the primary

Ideally, the ratio of f_YS, p to f_YS, s should be proportional to r_s/r_p. However, due to the non-linear effect produced by the partial eclipse stated above, the YS coefficient for the primary, assuming the same physical properties as the secondary except the size, is affected by a non-linear dependence on r_s. As a result, the ratio of f_YS, p to f_YS, s is better described as

$$\begin{array}{c}\hfill f_{\mathrm{YS},\mathrm{p}}={\left(\frac{r_{\mathrm{s}}}{r_{\mathrm{p}}}\right)}^{\alpha }{f}_{\mathrm{YS},\mathrm{s}}={\left(\frac{r_{\mathrm{s}}}{r_{\mathrm{p}}}\right)}^{\alpha -1}{f}_{\mathrm{YS},\mathrm{p}},\end{array}$$(12)

where f_YS, s and f_YS, p can be obtained from Eqs. (4) and (6), respectively.

We searched the best-fit α for the numerical results by minimizing the sum of the squares of the difference between simulation results and the fitted values, as shown in the left panel in Fig. 3. The right panel of Fig. 3 shows the values of α under different frequency ratios, ω/n, and thermal inertia. We find that α is insensitive to ω/n but shows a clear decreasing trend with increasing thermal inertia. The values of α are 1.81 ± 0.05, 1.77 ± 0.05, 1.75 ± 0.05, 1.70 ± 0.04, and 1.65 ± 0.03 for Γ = 100,  200,  300,  500, and 1000 tiu [J m⁻² s^−1/2 K⁻¹], respectively. One-kilometre-sized objects typically have a thermal inertia of around 300 tiu (Delbo’ et al. 2007), as seen in Didymos (Rivkin et al. 2023; Rozitis et al. 2024). Once α was determined, we obtained the modified f_YS, p solution. An example of the modified f_YS, p solution with α = 1.7 is shown in Fig. 2 (dashed gold curve). It fits the simulation much better than the simple solution with α = 1.

Fig. 3. Left: YS coefficient as a function of the secondary-to-primary size ratio. The primary asteroid has a radius of 1 km. The orbital period is 12 hours, and the spin periods of both components are 8 hours. The ratio of the primary-to-secondary YS coefficient is indicated by red triangles. The solid red line represents the best-fit curve: fYS, p/fYS, s = (rs/rp)1.80. The dashed red line shows the theoretical formula without accounting for the non-linear effect of partial eclipses: fYS, p/fYS, s = (rs/rp). Right: Value of α as a function of the frequency ratio and thermal inertia. The dashed black line represents the theoretical value, excluding the non-linear effects of partial eclipses. For a typical small binary system with Γ ∼ 200 tiu and a primary spin period of approximately three hours, α is around 1.8.

4.3. Implications

We have demonstrated that the BYP moves the primary asteroid towards a synchronous orbit, albeit on a longer timescale than the BYS. Together with the BYS (Zhou et al. 2024) and the BYORP effect (Ćuk & Burns 2005), this forms a unified description of radiative forces acting on binary asteroid systems, as summarized in Table 1: For an asynchronous component, the binary Yarkovsky effect is at play; for a synchronous component with an asymmetric shape, the binary YORP effect is dominant; and for a synchronous component with a symmetric shape, there is no significant secular radiative effect.

Based on this framework, we summarize the primary dynamics for different types of binary asteroids as follows. For synchronous binaries where both components are asynchronous, the BYS, the BYP, and tidal effects are all active. The BYS dominates over the BYP and, with the assistance of tidal effects, drives the secondary component towards a synchronous state. The synchronization timescale for a binary asteroid system with a primary radius of approximately 1 km is roughly 0.1 million years (Zhou et al. 2024). For larger objects, the tidal effect predominates over the BYS. The YORP effect may also contribute to synchronization, though its stochastic nature is not yet well constrained (Statler 2009; Golubov & Krugly 2012; Bottke et al. 2015; Zhou et al. 2022).

For singly synchronous binaries where the secondary is synchronous, which is the most common type of known binary asteroids, BYP, BYORP (acting on the secondary), and tidal effects are all active. The BYP is expected to shrink the orbit of most binary asteroids given that the majority of primary asteroids in the main belt or in near-Earth orbit rotate rapidly, close to the spin limit of 2.2 hours (Walsh & Jacobson 2015). The timescale of the BYP is typically an order of magnitude longer than that of the BYS and depends on the size ratio of the two components. The BYORP effect is predicted to act on a timescale ranging from 0.01 to 10 million years with a random direction and depends on factors such as shape (Ćuk & Burns 2005; McMahon & Scheeres 2010; Steinberg 2011; Jacobson & Scheeres 2011), surface roughness (Ćuk et al. 2024), and rotational state (Quillen et al. 2022).

For doubly synchronous binaries (i.e. both components are synchronous), the BYORP effect acts on both components. Tidal effects also become significant in eccentric orbits (Wisdom 2008; Goldberg & Batygin 2024) or for librating objects (Scheeres 2009; Jacobson et al. 2014).

While the magnitude of the BYORP effect remains largely uncertain due to the unconstrained shape and fine surface structures of most asteroids, the BYP can be estimated quickly using the formula in this work (see Eq. (12)). We have estimated the YS coefficients and orbital drift rates for known synchronous binary asteroids with detected spin periods and orbital periods (Pravec & Harris 2007; Warner et al. 2009; Pravec et al. 2012, 2016; Monteiro et al. 2023). As an order-of-magnitude approximation, we assumed spherical binaries that have relatively low inclinations, which ensures that these binaries experience eclipses during each mutual orbit. This assumption is appropriate for most binary asteroids (Pravec et al. 2012). The binary Yarkovsky effect is expected to be insensitive to shape, except in cases of extreme irregularity, which is unrealistic for binary asteroids. These bodies are expected to have relatively regular shapes following rotational disruption. All of these binary asteroids shrink the mutual orbit under the BYP due to their faster spin relative to the mean motion (see also the histogram in Fig. 2).

Figure 4 shows that for the known small binary asteroids, the absolute value of f_YS, p ranges from 10⁻⁴ to a few ×10⁻³ and the BYP-induced orbital drifting rate from −0.001 cm yr⁻¹ to −1 cm yr⁻¹. For comparison, the typical BYORP coefficient is ∼10⁻³ (McMahon & Scheeres 2010; Jacobson & Scheeres 2011) and decreases with reduced relative surface roughness (Ćuk et al. 2024). There is a clear trend that larger and colder binary asteroids have smaller values of f_YS, p and ȧ (located towards the left and bottom of the plot). Thermal inertia influences these results in a complex manner. Generally, within the range of 100–500 tiu, as thermal inertia increases, the Yarkovsky effect decreases for more distant objects (e.g. middle and outer main belt objects), while it strengthens for closer objects (e.g. NEOs and inner main belt objects).

Fig. 4. YS coefficient and YS-induced orbital drift rate for confirmed small binary asteroids, assuming thermal inertia of 100 tiu (left) and 500 tiu (right). The colours indicate the heliocentric semi-major axis, with bluer colours representing greater distances from the Sun. The size of each circle corresponds to the size of the primary asteroid. As thermal inertia increases, blue dots tend to move downwards and red dots upwards, although the outcome is a complex function of thermal inertia.

We compared our results with four binary asteroids that have observed orbital migration. The estimated binary Yarkovsky effect versus observed values are as follows: –0.05 versus −0.08 ± 0.02 cm yr⁻¹ (Scheirich & Pravec 2022; Scheirich et al. 2024; Naidu et al. 2024) for pre-impact Didymos, –0.079 versus −0.07 ± 0.34 cm yr⁻¹ (Scheirich et al. 2015) for 1996 FG₃, –0.13 versus −2.8 ± 0.2 cm yr⁻¹ (Scheirich et al. 2021) for 2001 SL₉, and –0.19 versus 1.2 cm yr⁻¹ (Scheirich et al. 2021) for 1999 KW₄. For pre-impact Didymos and 1996 FG₃, the BYP predicts orbital drift rates close to the observed values. However, for asteroid 2001 SL₉, our estimated BYP-induced orbital drift rate is an order of magnitude lower than the observed value. Furthermore, the outward drift of 1999 KW₄ cannot be explained by the BYP alone, suggesting the presence of additional mechanisms, such as a strong tidal effect, a BYORP effect on a synchronous secondary, or the BYS on an asynchronous secondary. Further observational data are needed to better understand the complex long-term dynamics of binary asteroids.

In Vokrouhlický et al. (2005)’s pioneering work on the Yarkovsky effect in binary asteroids, simulations of asteroid 2000 DP₁₀₇ showed that the mean transverse acceleration of the mutual orbital motion is ∼ − 6 × 10⁻¹⁵ m s⁻² assuming K = 0.1 W m⁻¹ K⁻¹, C = 800 J kg⁻¹ K⁻¹, and ρ = 1.7 g cm⁻³ for both components. This translates to an orbital drift rate of about ȧ ∼ −0.85 cm yr⁻¹. In their simulation, the primary and secondary are both spherical polyhedral shape models with the assumption of a synchronous secondary. Therefore, their results reflect the BYP. Our theoretical estimate of −0.27 cm yr⁻¹ is consistent in both sign and order of magnitude with the previous simulation result despite some deviation that may arise from the system’s eccentricity.

5. Conclusion

This work, together with our previous work on the BYS (Zhou et al. 2024), completes the basic theoretical framework of the binary Yarkovsky effect. The BYP is shown to modify the mutual orbit after the secondary asteroid reaches synchronization. The timescale of the BYP is generally an order of magnitude longer than that of the BYS and depends on the size ratio between the secondary and primary.

We propose an empirical modified formula to estimate the BYP: applying the traditional binary Yarkovsky formula (Eq. (6)) and then multiplying it by (r_s/r_p)^{(α − 1)} (see Eq. (12)). Our numerical results indicate that α is relatively insensitive to the frequency ratio, ω/n, but decreases with increasing thermal inertia. For a typical small binary asteroid system with Γ = 200 tiu, the best-fit value for α is approximately 1.7. We summarize the primary mechanisms for binary asteroid systems as follows: the BYS and tidal effects are active for asynchronous binaries, while BYS, BYORP, and tidal effects operate in singly synchronous binaries, which are the most commonly observed binary asteroids.

We estimated BYP-induced drift rates for known small binary asteroids with primary radii r_p < 1 km. These drift rates range from –0.001 to –1 cm yr⁻¹. For middle and outer main belt objects, the BYP tends to decrease with increasing thermal inertia, whereas for NEOs and inner main belt objects, the BYP increases with increasing thermal inertia.

We compared our results with the observed orbital drift rates of four binary systems. Our findings are consistent with the pre-impact Didymos and 1996 FG₃ systems but show discrepancies for the 2001 SL₉ and 1999 KW₄ systems. We suggest that a complicated model involving BYP, BYORP, and tidal effects is required to fully understand the long-term dynamics of binary systems.

Acknowledgments

I thank the referee for the constructive suggestions. I thank David Vokrouhlický for his valuable suggestions. I am also grateful to Patrick Michel and Seiji Sugita for their support of this research. I would like to acknowledge the mobility aid from the Université Côte d’Azur and funding support from the Chinese Scholarship Council (No. 202110320014).

References

Appendix A: Discussion on non-zero inclination cases

The inclination i is defined as the angle between the vectors of the mutual orbit and the heliocentric orbit. A non-zero inclination complicates the thermal perturbation caused by the shadow. We offer an estimate of the non-zero inclination and justify the importance of BYP on real binary asteroids. A comprehensive numerical investigation on how a non-zero inclination affects both BYS and BYP is out of the scope of this paper and is left for future study due to its complexity. Our following simple estimation is based on the principle that the strength of the YS effect is generally proportional to the solar flux loss due to the shadow. We will use an inclination-dependent factor f_i to describe the ratio of the YS effect for a non-zero inclination relative to the case with zero inclination.

To begin, we considered how inclination influences the YS effect on the secondary. When the inclination is non-zero, the shadow cast by the secondary on the primary no longer remains fixed at the equator but oscillates around it, reaching a maximum altitude, β. Because the primary is approximately spherical, the duration it experiences shadowing is roughly proportional to cos β, leading to a weakening factor of

$$\begin{array}{c}\hfill f_{i,\mathrm{s}}=cos\beta .\end{array}$$(A.1)

The parameter β can be estimated via simple geometry,

$$\begin{array}{c}\hfill sin\beta =\frac{atani}{r_{\mathrm{p}}}.\end{array}$$(A.2)

A sin β < 1 naturally necessitates a low inclination, i < arctan(r_p/a)≃r_p/a, for the occurrence of the shadow.

The case for the primary would be a bit different. The shadow cast by the secondary on the primary is above the equator rather than at the equator when the inclination is not zero. These high-altitude areas affected by the shadow have lower temperatures than the equator, as the total incident solar energy decreases by a factor of cos β. Consequently, the overall YS effect would be reduced by a factor of

$$\begin{array}{c}\hfill f_{i,\mathrm{p}}={cos}^2\beta .\end{array}$$(A.3)

Recall that the binary Yarkovsky effect also includes the pY effect in addition to the YS effect. The strength of pY effect is approximated by the YS coefficient multiplied by −πr_p/a for the secondary and −πr_s/a for the primary. Therefore, compared to Eq. 2, the total binary Yarkovsky coefficient becomes more complicated for non-zero inclination cases:

$$\begin{array}{c}\hfill f_{\mathrm{BY},\mathrm{s}}=(cos\beta -\frac{\pi r_{\mathrm{p}}}{4a})f_{\mathrm{YS},\mathrm{s}}\end{array}$$(A.4)

for the secondary and

$$\begin{array}{c}\hfill f_{\mathrm{BY},\mathrm{p}}=({cos}^2\beta -\frac{\pi r_{\mathrm{s}}}{4a})f_{\mathrm{YS},\mathrm{p}}\end{array}$$(A.5)

for the primary. We note that the typical values for f_YS, s and f_YS, p are 0.01 (Zhou et al. 2024) and 0.001 (Fig. 2), respectively.

Figure A.1 shows the value of the binary Yarkovsky coefficient as a function of inclination, described by Eqs. A.4 and A.5. It is seen that YS decreases with the inclination and vanishes with a critical inclination, within which the eclipse always occurs during a mutual orbit. In contrast, the pY effect remains independent of inclination, as it is driven by radiation from the companion object within the binary system. There is a specific inclination at which pY cancels out YS, resulting in the disappearance of the binary Yarkovsky effect.

Fig. A.1. Binary Yarkovsky coefficient as a function of the inclination for the secondary (red) and the primary (blue). The dashed black line is the critical inclination (i.e. arctan(rp/a)), within which the eclipse always occurs during a mutual orbit. The values of parameters to calculate fBY (Eqs. A.4 and A.5) are set as follows: a/rp = 2.5, fBY, s = 0.01 and fBY, p = 0.001. In this figure, a positive value of fBY represents the migration direction towards the synchronous orbit and a negative value denotes the opposite direction to the synchronous orbit.

When the inclination is low and YS dominates, the binary Yarkovsky effect tends to synchronize the orbit. Conversely, at high inclinations where pY is dominant, the effect causes de-synchronization. The assumption of zero inclination in the main text maximizes the binary Yarkovsky effect and ensures mutual orbital synchronization. Although the strength of the binary Yarkovsky effect diminishes with increasing inclination, it does not change by more than an order of magnitude except near the critical inclination (i.e. arctan(r_p/a)), justifying the importance of the binary Yarkovsky effect for most of the binary asteroids.

We note that a more detailed numerical investigation of the role of the inclination is required in the future, especially considering the non-zero obliquities that complicate the problem significantly.

All Tables

All Figures

Fig. 1. Left: Snapshots of the temperature field of the binary asteroids. Diagrams (a)–(d) illustrate the anti-clockwise orbit of the secondary asteroid around the primary. Both the primary and secondary have spin rates of 3 hours. The other properties of the binary system are detailed in Sect. 4.1. In phase (a) the primary is partially eclipsed by the secondary, while in phase (c) the primary fully eclipses the secondary. Right: Tangential accelerations due to thermal forces for the primary (blue) and secondary (red). The eclipse periods are represented by shaded areas.

In the text

Fig. 2. YS coefficient for the primary as a function of the ratio of spin velocity to mean motion. The primary asteroid has a radius of 1 km, and the secondary has a radius of 0.3 km. Gold circles represent the numerical results. The solid gold curve shows the simple analytical solution from Eq. (6), while the dashed gold curve shows the modified solution from Eq. (12) with α = 1.7. The ωp/n distribution of confirmed binary asteroids is depicted in the blue histogram, which is truncated at 10.

In the text

Fig. 3. Left: YS coefficient as a function of the secondary-to-primary size ratio. The primary asteroid has a radius of 1 km. The orbital period is 12 hours, and the spin periods of both components are 8 hours. The ratio of the primary-to-secondary YS coefficient is indicated by red triangles. The solid red line represents the best-fit curve: fYS, p/fYS, s = (rs/rp)1.80. The dashed red line shows the theoretical formula without accounting for the non-linear effect of partial eclipses: fYS, p/fYS, s = (rs/rp). Right: Value of α as a function of the frequency ratio and thermal inertia. The dashed black line represents the theoretical value, excluding the non-linear effects of partial eclipses. For a typical small binary system with Γ ∼ 200 tiu and a primary spin period of approximately three hours, α is around 1.8.

In the text

Fig. 4. YS coefficient and YS-induced orbital drift rate for confirmed small binary asteroids, assuming thermal inertia of 100 tiu (left) and 500 tiu (right). The colours indicate the heliocentric semi-major axis, with bluer colours representing greater distances from the Sun. The size of each circle corresponds to the size of the primary asteroid. As thermal inertia increases, blue dots tend to move downwards and red dots upwards, although the outcome is a complex function of thermal inertia.

In the text

Fig. A.1. Binary Yarkovsky coefficient as a function of the inclination for the secondary (red) and the primary (blue). The dashed black line is the critical inclination (i.e. arctan(rp/a)), within which the eclipse always occurs during a mutual orbit. The values of parameters to calculate fBY (Eqs. A.4 and A.5) are set as follows: a/rp = 2.5, fBY, s = 0.01 and fBY, p = 0.001. In this figure, a positive value of fBY represents the migration direction towards the synchronous orbit and a negative value denotes the opposite direction to the synchronous orbit.

In the text