© The Authors 2025

1 Introduction

A small rotating body orbiting around the Sun is heated by the Solar radiation, and a specifically anisotropic temperature distribution appears on its surface. The absorbed energy is then released into space via thermal re-radiation, and the net recoil force of the thermal radiation may influence the asteroid’s orbital evolution.

Such a thermodynamic effect was first proposed by and later named after Polish Russian engineer Yarkovsky based on the erroneous ether theory in 1901, but it did not attract widespread attention (see e.g. Beekman 2005, for a brief review). Öpik (1951) reintroduced this concept and noted that under the Yarkovsky effect, the semi-major axis of a small-sized object increases if it rotates progradely and decreases in the case of retrograde rotation. Then, after Peterson (1976) proposed a two-step mechanism based on this effect to explain the origin of meteorites, the Yarkovsky effect came into focus. Rubincam (1987) discovered the seasonal Yarkovsky effect while explaining the observed orbital decay of the satellite LAGEOS. Finally, the diurnal component (arising from the body’s rotation) and the seasonal component (from orbital revolution) together constitute the complete Yarkovsky effect.

The Yarkovsky effect can be calculated either by analytical method or by numerical simulation. The former, often referred to as the ‘linear theory’, is based on the assumption that the temperature at any location on the surface of a small object does not differ much from its average temperature. Vokrouhlicky (1998a) linearised the emission term (∝T⁴) in the equation of boundary condition, and obtained a complete linear solution for the diurnal Yarkovsky effect of a spherical body with an arbitrary spin axis obliquity. Later, this linear theory was extended to the biaxial ellipsoidal body (Vokrouhlicky 1998b). Vokrouhlický & Brož (1999) reduced a small celestial body to a high-thermal-conductivity spheroid covered with a thin lowconductivity surface layer and obtained an analytical solution for the seasonal Yarkovsky effect. Finally, a complete solution that includes both the diurnal and seasonal components was proposed by Vokrouhlický (1999).

Numerical computations of the Yarkovsky effect were developed with different thermal models for small bodies. Constrained by limitations in computing capability and observational accuracy, early thermal models for small celestial bodies often introduced certain simplified assumptions and approximations. The standard thermal model (STM) considers asteroids as smooth, non-rotating spheres, while the fast rotation model (FRM) is applicable to rapidly rotating asteroids with large thermal inertia (Lebofsky et al. 1986; Lebofsky & Spencer 1989). Both STM and FRM are primarily applicable to mainbelt asteroids, while Harris (1998) established the near-Earth asteroid thermal model (NEATM). All these three models assume zero thermal emission on the night side of an asteroid. Spitale & Greenberg (2001) and Spitale & Greenberg (2002) used a finite-difference method to solve the three-dimensional (3D) heat conduction equation in an isotropic spherical body. They simultaneously calculated the diurnal and seasonal components, and their effects on the semi-major axis, eccentricity, and longitude of periapse. Lagerros (1996) proposed the thermal physical model (TPM), approximating asteroids as triaxial ellipsoids and establishing a one-dimensional (1D) heat conduction model. Subsequently, Lagerros (1997) further extended the model to irregularly shaped bodies, considering the effect of shadowing and self-heating. Rozitis & Green (2011) proposed a well-developed advanced thermophysical model (ATPM), which incorporates the effects of shadowing, scattering of sunlight, self-heating and previously neglected thermal-infrared beaming caused by surface roughness. In ATPM, the lateral heat conduction is still neglected and only 1D heat conduction perpendicular to the surface is taken into account. Then Rozitis & Green (2012) and Rozitis & Green (2013) adapted ATPM to both Gaussian random spheres and shapes of real asteroids to predict the semimajor axis drift caused by the Yarkovsky effect in the presence of surface roughness and global self-heating. Basart et al. (2012) modelled the 3D heat transfer incorporating the heat transfer and surface-to-ambient radiation modules in COMSOL¹. Using a sphere as an example, the study computed the dependence of the diurnal Yarkovsky effect on radius, rotation period and obliquity of incident radiation. Xu et al. (2022) calculate the Yarkovsky effect of real asteroids and propose an index called ‘effective area’ to measure the influence of asteroids’ shape, which has been shown to be a fairly good index of the Yarkovsky effect on irregularly shaped asteroids.

The linear theory gives an analytical solution for asteroids of regular shapes such as spheres and biaxial ellipsoids. Whereas for irregularly shaped asteroids, we can either approximate them as regular shapes, which introduces certain errors, or numerically compute the surface temperature and then the recoil force of thermal irradiation, which gives more accurate results but requires a heavy calculation. The shape index can measure the effect of shape on the diurnal Yarkovsky effect and the application of such an index allows for the rapid and accurate estimation of Yarkovsky effect on arbitrarily shaped asteroids. We notice that its applicability has not been clearly clarified, and its accuracy can be further improved. In this paper, we improve the definition of the shape index by taking into account the effects of projected shadows and self-heating, and we also investigate its applicability under different thermal parameters. The rest of this paper is organised as follows. In Section 2, we briefly review the heat conduction equation, and set up the numerical models. In Section 3, we define the shape index, improve it by taking into account the shadowing effect, and test the range of parameters for which the shape index works well. We also analyse the influence of lateral heat transfer, scattering and self-heating effects on the performance of the shape index. Finally we conclude our investigation in Section 4.

2 Model and methods

2.1 Yarkovsky effect

The Yarkovsky effect of an asteroid arises from the recoil force of the surface thermal radiation, which depends on the temperature distribution. The temperature in turn is determined by the heat conduction as well as the radiation absorbed and emitted by the surface. Since generally the diurnal Yarkovsky effect is much stronger than the seasonal effect (see e.g. Vokrouhlický 1999; Xu et al. 2020), we focus on the former in this paper. In the following, the Yarkovsky effect simply refers to the diurnal effect, unless otherwise stated.

In a rotating body-fixed co-ordinate system, the 3D heat conduction equation can be written as: $$\rho C\frac{\mathrm{\partial }T}{\mathrm{\partial }t}=K{\mathrm{\nabla }}^{2}T,$$(1)

where T represents the temperature, and ρ, C, and K are the density, specific heat capacity, and thermal conductivity, respectively. The boundary conditions are the conservation of energy: $$ϵ\sigma T^4+K(\mathbf{n}\cdot \mathrm{\nabla }T)=\alpha \mathcal{E},$$(2)

where ϵ, σ, α, and $\mathcal{E}$ are the emissivity, Stefan–Boltzmann constant, absorption coefficient, and external radiation flux, respectively, while n denotes the unit vector normal to the surface. Numerically, the heat conduction equation can be solved by the finite element method and the surface temperature is then obtained.

For a Lambert’s isotropic thermal emission geometry, the recoil force generated by the small object’s radiation is: $$\mathbf{F}=-\int \frac{2}{3c}ϵ\sigma T^4\mathbf{n}\mathrm{d}S.$$(3)

Substituting the force into the Gauss perturbation equation, we obtain the drift rate of the semi-major axis due to the Yarkovsky effect: $$\frac{\mathrm{d}a}{\mathrm{d}t}=\frac{2a^{3/2}}{\sqrt{\mu (1-e^2)}}[F_{r}e\sin v+F_t(1+e\cos \nu )],$$(4)

where μ = G (M_⊙ + m) is the reduced mass of the Sun and the small body, e is the orbital eccentricity, and v the true anomaly. In addition, F_r and F_t are the radial and tangential components of the Yarkovsky force (F) respectively.

Equation (4) gives the instantaneous acceleration of the semi-major axis (da/dt), which might change with time as the asteroid rotates and revolves along its orbit. Generally, the rotation period, P, is much shorter than the revolution period, P_orb, so it is reasonable to assume that the temperature distribution on the asteroid’s surface reaches a state of dynamic equilibrium at any phase of its orbit. Therefore, the surface temperature can be calculated numerically and the corresponding Yarkovsky force can be averaged over one rotation period (P) at each position along the orbit. And the semi-major axis drift rate in long term should be averaged again over the revolution period (P_orb).

2.2 Finite element method

The linear theory of Yarkovsky effect works quite well for objects of regular shape, but it might bring considerable error if we treat an irregular asteroid as an ellipsoid that the linear theory can handle (Xu et al. 2022). In this paper, we numerically calculate the temperature of an asteroid of any given shape using the software COMSOL, in which the heat transfer equation is computed in 3D model; that is, the heat conduction in both lateral and perpendicular directions to the asteroid’s surface is taken into account. The model also incorporates shadowing, scattering of sunlight and self-heating effect. But the roughness-induced thermal-infrared beaming effect is not included because the surface is always set as Lambertian surface.

Our main interest in this paper is to quantify the influence of an asteroid’s shape on the Yarkovsky effect. Currently, the number of asteroids whose shapes are well known is still limited. We obtain from the Planetary Data System² 34 real asteroids as our samples of irregularly shaped objects. For the sake of direct comparison, these shape models are then scaled to have the same volume as a sphere of radius R = 10 m.

The triangular facet shape models of these samples are then constructed and fed the software. Requiring a reasonable accuracy of temperature, the resolution of these shape models is automatically set by COMSOL. Thousands of meshes for the surface of each model were generated. An asteroid is assumed to be a core composed of free tetrahedra covered by a shell of certain thickness. The heat exchange occurs barely not in the core, but only in the shell (see e.g. Xu et al. 2022). Composed of triangular prisms, the shell is set to be 5 times thick as the thermal penetration depth, l_d, and is then divided into 25 layers to make the mesh for numerical calculations. We note that the penetration depth, $l_{\mathrm{d}}=\sqrt{KP/(2\pi \rho C)}$, is the e-folding depth at which the temperature variation amplitude is decreased by a factor of 1 / e.

We set the asteroid at the initial moment as an isothermal object, with the temperature T_ini, determined by the balance between the outward thermal radiation and the absorption of solar radiation; that is, $4\pi R^2\cdot ϵ\sigma T_{\text{ini }}^4=\pi R^2\cdot \alpha \mathcal{E}$. The time step used to calculate the variation of temperature is taken to be one hundredth of the rotation period (P/100) in the simulation. Generally, the surface temperature reaches a dynamic equilibrium after several cycles of rotation, then the temperature at any point on the surface varies periodically. When the difference of temperature distribution at a certain moment from the one at the same phase in the next rotational cycle is less than 0.1%, the dynamic equilibrium state is regarded as reached. The number of rotation periods required for an asteroid to reach the equilibrium depends on thermal parameters. Particularly, an asteroid attains the equilibrium quickly when the value of thermal conductivity is small, and vice versa.

The actual thermal processes could be very complex. The thermal parameters are not constants but functions of temperature. The albedo may be a function of wavelength, which means that the surface reflects the self-heating radiation in infrared with a different efficiency from the Solar radiation that is described by the Bond albedo. However, to simplify our analyses, we set in this paper the thermal parameters as constants and take the albedo in the infrared band as the Bond albedo. The symbols and values of the parameters we use in this article are summarised in Table 1.

The rotation period P = 1800 s, in Table 1 is relatively short compared with the observed periods in common asteroids. We note that it does not affect our analyses, but it brings great convenience in our computation using COMSOL. When the thermal conductivity is low, at a fixed time step (P/100), a shorter rotation period serves to reduce errors in the calculation of temperature distribution. Additionally, for asteroids of large size (R ≫ l_d), the Yarkovsky effect is determined by the ‘thermal parameter’ Θ (Vokrouhlicky 1998a) defined as: $$\mathrm{\Theta }=\frac{\mathrm{\Gamma }\sqrt{\omega }}{ϵ\sigma T_{\star }^3},$$(5)

where $\mathrm{\Gamma }=\sqrt{\rho CK}$ is the thermal inertia, ω = 2π/P is the angular velocity of the body’s rotation, and T_⋆ is the subsolar temperature defined by $ϵ\sigma T_{\star }^4=\alpha \mathcal{E}$. Therefore, the value of period, P, in fact can be traded with other parameters like ρ and C, and the impact of a too-short rotation period can be mitigated by adjusting other parameters. In this paper, we varied the thermal conductivity across four orders of magnitude, ensuring that the choice of rotation period does not compromise the reliability of our conclusions.

3 Shape index

We have found in our previous paper (Xu et al. 2022) an approximately linear dependence of the diurnal Yarkovsky effect on the ‘effective area’ of an irregularly shaped asteroid, in which the Yarkovsky effect was measured by the drift rate of semi-major axis, da/dt. And the effective area was used to assess the characteristics of an asteroid’s shape that affects the Yarkovsky effect. Below we briefly revisit the concept of the effective area and improve it by taking into account some geometrical and physical effects that were ignored previously.

3.1 Definition of shape index

Since the rate of an asteroid’s migration caused by the Yarkovsky effect depends explicitly on its spin obliquity, γ, as da/dt ∝ cos γ (Vokrouhlicky 1998a; Xu et al. 2022), we describe the strength of the Yarkovsky effect in the simplest case with the spin axis being perpendicular to the orbital plane; that is, γ = 0°. In this case, the equatorial plane of the body coincides with the orbital plane.

For a polyhedral shape model consisting of N facets, if the i-th surface element with an area, a_i, is illuminated by the Sun at moment, t, then the Solar radiation energy absorbed by the surface element at this moment is: $$E_{\mathrm{i}\mathrm{n}}=-\alpha \mathcal{E}{a}_i({\mathbf{n}}_i\cdot {\mathbf{n}}_{\mathrm{i}\mathrm{n}}),$$(6)

where n_i and n_in denote the directions of this surface element and the Solar incidence. We note that these directions are defined in an asteroid-centred frame, in which n_i rotates and n_in is fixed.

When the temperature distribution on the surface of an asteroid has reached the dynamic equilibrium, the temperature of any surface element should satisfy T (t + P) = T(t), and the amounts of heat absorbed and released by a certain surface element should be equal in every rotation period. For an asteroid of size much larger than the penetration depth (R ≫ l_d), the heat budget of any surface element is approximately completed locally but not globally; that is, only in the same surface element. Ignoring the heat exchange between neighbouring surface elements, a surface element absorbs the Solar radiation, stores the energy in it, and then re-radiates the energy from the same surface element after a delay of time, δt. Given the rotation frequency, ω, the delay in time, δt, means that the re-radiation will be omitted after the surface element rotates by a delay angle θ = ωδt. Therefore the recoil force produced by the thermal radiating of energy absorbed by this surface element at this moment is in the direction n_rf = −n_ie^iθ, where $\mathrm{i}=\sqrt{-1}$ and e^iθ denotes a rotation by angle θ. Since the delay angle, θ, is determined by the thermal parameter, Θ, for a homogeneous asteroid, we expect that all surface elements have the same δt and therefore the same delay angle.

In our model, the Solar incidence always lies on the equatorial plane of the body, thus the tangential direction of the circular orbit is n_t = n_ine^iπ/2. Finally, the magnitude of the recoil force projected on the tangential direction produced by this surface element at this moment is proportional to the energy released in the same direction: $$E_{\text{out }}\propto -({\mathbf{n}}_{\mathrm{r}\mathrm{f}}\cdot {\mathbf{n}}_{\mathrm{t}})E_{\text{in }}\propto \alpha \mathcal{E}{a}_i({\mathbf{n}}_i\cdot {\mathbf{n}}_{\text{in }})({\mathbf{n}}_i{e}^{\mathrm{i}\theta }\cdot {\mathbf{n}}_{\text{in }}{e}^{\mathrm{i}\pi /2}).$$(7)

The net recoil force of the thermal radiation produced by the i-th surface element in a rotation period can be obtained by averaging over a rotation.

If the shape of an asteroid is fully convex, there will be no projected shadows and self-heating effect. Each surface element is exposed to the Solar radiation for exactly half rotation, and the average contribution of each facet to the Yarkovsky effect over a rotation period is proportional to: $$F_{t,i}\propto \frac{1}{P}{\int }_{t_0}^{t_0+P/2}\alpha \mathcal{E}{a}_i({\mathbf{n}}_i\cdot {\mathbf{n}}_{\mathrm{i}\mathrm{n}})({\mathbf{n}}_i{e}^{\mathrm{i}\theta }\cdot {\mathbf{n}}_{\mathrm{i}\mathrm{n}}{e}^{\mathrm{i}\pi /2})\mathrm{d}t,$$(8)

where t₀ is the sunrise moment of the i-th surface element. Denote all unit vectors by direction cosines. Particularly, the Solar incident is in n_in = (cosβ_x, cosβ_y, cosβ_z). In the body-centred frame rotating with the asteroid, the direction cosine can be reduced to n_in = (cosβ, sinβ, 0), where β is the angle between the x axis and n_in. The average over time in a rotation period in Eq. (8) is equivalent to an average over the angle β, and further algebraic calculations show: $$\begin{array}{rl}{F}_{t,i}\propto & \frac{1}{\pi }{\int }_{-\pi /2}^{+\pi /2}\alpha \mathcal{E}{a}_i({\mathbf{n}}_i\cdot {\mathbf{n}}_{\text{in }})({\mathbf{n}}_i{e}^{\mathrm{i}\theta }\cdot {\mathbf{n}}_{\text{in }}{e}^{\mathrm{i}\pi /2})\mathrm{d}\beta \\ & =\frac{1}{\pi }\alpha \mathcal{E}{a}_i\sin \theta {\int }_{-\pi /2}^{+\pi /2}{({\mathbf{n}}_i\cdot {\mathbf{n}}_{\text{in }})}^2\mathrm{d}\beta \\ & =\frac{1}{4}\alpha \mathcal{E}{a}_i\sin \theta {({\mathbf{n}}_i\cdot {\mathbf{n}}_{\text{noon }})}^2.\end{array}$$(9)

Obviously, in this integration, we have set β = 0 when this surface element is at its local noon, and n_noon = n_in(β = 0).

The contributions from all surface elements can be summarised to represent the overall Yarkovsky force, and the summation can then be normalised by dividing the corresponding value obtained from a spherical asteroid of the same thermal parameters and with the same bulk volume. The constant ‘parameters’ including α, $\mathcal{E}$, and the delay angle sin θ, are cancelled out in the normalization, and finally a shape index S₁ is defined as such a normalised summation. In fact, it’s the ratio of the ability of an object with an arbitrary shape to produce the Yarkovsky effect to that of a sphere with the same volume and thermal parameters, and it reads: $$S_1=\frac{A}{4}\sum _i^N{a}_i{({\mathbf{n}}_i\cdot {\mathbf{n}}_{\mathrm{n}\mathrm{o}\mathrm{o}\mathrm{n}})}^2,$$(10)

where A = [6/(πV²)^1/3 and V is the volume of the asteroid. Surely, according to the definition, S₁ = 1 for a spherical asteroid. The index S₁ describes qualitatively the influence of shape on the Yarkovsky effect. We note that S₁ has the same physical essence as the ‘effective area’ defined in Xu et al. (2022), only except the former has been normalised and thus is dimensionless.

So far, we have assumed that the delay (δt in time or θ in angle) between the absorption and emission of energy on all surface elements are the same on an asteroid. To verify this assumption, using COMSOL we numerically calculate the surface temperature of two model asteroids, namely, a spherical asteroid and the asteroid YORP (scaled to have the same volume as a sphere of radius 10 m). The delay occurs all the time, but in practice, it is not easy to define it at each moment. Instead, we calculate the radiation force of each surface element and average it over a rotation period. As comparison, we set the heat conductivity K = 0 (therefore no heat conduction and no delay at all) and calculate the same averaged force from the same surface element. Then the delay angle is defined as the angular difference between directions of the averaged forces in these two cases. We plot the statistics of delay angles of two asteroid models in Fig. 1.

Since the force component perpendicular to the orbital (equatorial) plane is not of interest in this paper, the delay angles and the magnitudes of the force shown in Fig. 1 are in fact their projections on the orbital plane. We note that considering different sizes of facets used in calculations, we use the pressure (force per area) instead of the force to represent the magnitude of radiation force in Fig. 1. With the same thermal parameters for both sphere and YORP model, regardless of their shapes, the delay angles of surface elements from both models display a similar distribution. Most of the delay angles are concentrated at 13° ~ 14° as shown in Fig. 1a where thermal parameter Θ = 0.736 is adopted. In addition, these facets with nearly the same delay angle have the largest radiation pressures and contribute the most to the Yarkovsky effect. Although the delay angles of some surface elements that are located in the polar region or greatly affected by the shadowing effect (in the YORP model) deviate significantly from 13° ~ 14°, their radiation pressures are very low, and they contribute very little to the Yarkovsky effect.

More intuitively, we show in Fig. 1b the surface temperature distribution on the sphere model and the radiation forces at different latitudes. Because of the geometric symmetry of the sphere model, all facets at a certain latitude on a sphere are identical. Thus the direction of the radiation force at different latitudes tells directly the delay angle (from the subsolar point). From Fig. 1b, we see that the delay angles in the polar region are larger than those in low latitudes. However, the magnitudes of the radiation forces at high latitude region are much smaller than the ones at low latitude region. Therefore, the total radiation force are mainly determined by the forces produced in the equatorial region, and their delay angles are nearly identical to each other. So, the assumption that all surface elements of a body have the same delay angle, θ, determined by the thermal parameter, Θ, is valid.

In the definition of the shape index S₁, in Eq. (10), the illumination of a surface element is determined only by its own direction, which is true only if the body is fully convex. But, small bodies of metres to kilometres in radii, which are most affected by the Yarkovsky effect, are likely to have irregular shapes due to their weak self-gravity, and they might not be simply approximated as fully convex. Particularly, many of them tend to have non-negligible projected shadows, under which the surface area apparently has lower temperature and thus weaker thermal radiation. Taking into account the shadowing effect, we may improve the shape index as follows.

For surface elements that are in the projected shadows of other elements at a given moment t, they are unable to absorb energy from the Solar radiation, and thus they are assumed to have barely thermal radiation and contribute negligibly to the Yarkovsky force before the Sun rises (locally for them). Therefore, a correction to Eq. (6) can be made as: $$E_{\text{in }}\propto -\alpha \mathcal{E}{a}_i[1-s_i(t)]({\mathbf{n}}_i\cdot {\mathbf{n}}_{\text{in }}),$$(11)

where s_i(t) indicates whether the i-th surface element is in shadow at time t, with s (t) = 1 if it is and s (t) = 0 if it is not. Following the same derivation as shape index S₁, a corrected shape index, S₂, can be defined as: $$S_2=\frac{A}{P}{\int }_0^P\left[\sum _i^N{a}_i(1-s_i(t)){({\mathbf{n}}_i\cdot {\mathbf{n}}_{\mathrm{i}\mathrm{n}})}^2\right]\mathrm{d}t.$$(12)

To calculate S₂, we need to run the shadow tests and compute the integration over a rotation period.

As we have shown above, the power of thermal radiation in the direction of orbital velocity is proportional to the shape index. We note that the drift rate of semi-major axis (da/dt) of an asteroid is actually proportional to the changing rate of the mechanical energy (the power) of an orbit, which in turn is driven by the recoil force of the thermal radiation. Therefore, a linear relationship between the Yarkovsky effect (measured by da/dt) and the shape index is expected.

Fig. 1 Directions and magnitudes of the averaged radiation forces of surface elements on sphere model and YORP model. A thermal conductivity K1 = 0.0015 W m−1 K−1 (correspondingly Θ = 0.736) is adopted for both models. The direction and magnitude have been projected in the orbital (equatorial) plane here (see text). (a) The magnitudes of radiation forces of surface elements versus the delay angle (upper panel), and the frequency of delay angles of all surface elements (lower panel). The delay angle of a surface element is the angle between the radiation force averaged over a rotation period and the corresponding averaged force of the same surface element but with K = 0. Blue and orange indicate the results for sphere and YORP model, respectively. (b) The direction and magnitude of the radiation force at different latitudes on the surface of a sphere. The longitude 0° is at noon. The longitude of the black dot at certain latitude indicates that the radiation force at this latitude points to the direction given by the longitude. The size of black dots represents the magnitude of the radiation force (in N/m2, see text) projected in the orbital plane. The yellow circle indicates the direction and magnitude (F/4πR2) of the total radiation force. The background colour gives the surface temperature.

3.2 Linear relation

The relationship between the index S₁ and the strength of the diurnal Yarkovksy effect has been investigated by Xu et al. (2022), in which a linear dependence of the semi-major axis drift rate on S₁ has been verified using 34 shape models of real asteroids.

When the shape of an asteroid is known, the indices S₁ and S₂ can be easily computed according to their definitions in Eqs. (10) and (12). While the surface temperature of an asteroid with given thermal parameters can be calculated using COMSOL and thus the Yarkovsky effect (represented by the drift rate of semi-major axis) can be obtained. All 34 models adopted in Xu et al. (2022) have been assumed to spin around the z axis. It is true that most of asteroids do rotate around the principal axis (that is the z axis) except for those tumblers. However, to extend the sample size of shape models, in this paper we assume that these asteroids also rotate around the x axis and y axis (while the Sun is always on the equatorial plane). Since rotating around different axes produces different Yarkovsky effects as well as different shape indices, by this way we finally have 102 (= 34 × 3) shape models.

We adopt the model of asteroids as described in Section 2 with parameters given in Table 1. A thermal conductivity K = 0.0015 W m⁻¹ K⁻¹, which is the typical value for the regolithcovered asteroids (Farinella et al. 1998), is adopted here. With these parameters, the thermal parameter Θ = 0.736 and the scaled radius of the body R′ = R/l_d ~ 10⁴, and an asteroid experiences approximately 50 rotation periods before the thermal equilibrium is reached across the body.

We present in Fig. 2a the semi-major axis drift rates (da/dt) versus the two shape indices (S₁, S₂) for these 102 shape models. The linear fits to the relationship between da/dt and S_1,2 are also plotted in Fig. 2a, and the residuals of the fits are shown in Fig. 2b.

Just as in Xu et al. (2022), the linear dependence of da/dt on the shape index S₁ (cyan diamonds in Fig. 2a) can be easily recognised. The linear fit has a fairly impressive coefficient of determination of R² = 0.984 and the mean standard deviation is σ_ϵ1 = 3.8 × 10⁻⁴ AU/Myr. However, the distribution of residuals (cyan diamonds in Fig. 2b) exhibits anomalies, and some significant outliers, all with negative values, can be seen in the lower part of the panel.

We suspect that the obvious deviation from the linear relation is due to the shadowing effect that has been neglected in S₁. In fact, when a facet on the asteroid surface is shadowed by other facets during rotation, it is unable to absorb energy from the solar radiation, and accordingly, its contribution to the diurnal Yarkovsky effect is weakened. For an asteroid that deviates significantly from fully convex shape, ignoring the effect of projected shadow may lead to a significant overestimation of its shape index S₁. An overestimated S₁ corresponds to a biased high semi-major axis drift in the fit, thus results in a negative residual. We examine the shape models that are clearly outliers in Fig. 2b (blue diamonds for S₁), and find that they are really affected by the shadowing effect. As examples, we present four shape models with the most anomalous residuals in Fig. 2c, all of which, as expected, are extremely irregular shapes with significant projected shadows.

Taking the shadowing effect into account, the shape index S₂ behaves better in indicating the strength of Yarkovsky effect than S₁, demonstrating a more distinct linear relationship with da/dt (orange dots in Fig. 2a). All the points, including those outliers in the fitting results of S₁, now fit closely the line. The coefficient of determination of the linear fit now is R² = 0.995, and the corresponding residuals in Fig. 2b form a normal distribution, with the mean standard deviation decreasing to σ_ϵ2 = 2.2 × 10⁻⁴ AU/Myr.

Fig. 2 Relationship between the semi-major axis drift rate (da/dt) and the shape index (S1,2). (a) da/dt versus S1 (cyan diamonds) and S2 (orange dots). (b) Fitting residuals. (c) Four shape models with large residuals: YORP, Kleopatra, Ida, and Golevka.

3.3 Dependence on K

When defining the shape index, the absorption, conduction and radiation of energy via a certain surface element are assumed to occur only in the same surface element; that is, the lateral heat transfer has been neglected. When establishing the linear relationship (Fig. 2) between the diurnal Yarkovsky effect and the shape indices, a thermal conductivity K₁ = 0.0015 W m⁻¹ K⁻¹ was chosen, under which the heat diffusion should be so slow that the assumption of null lateral heat transfer between surface elements is fairly reasonable. However, for basalt or iron-rich asteroids with much higher thermal conductivity, the lateral heat transfer cannot be ignored and the temperature distribution on the surface may be significantly influenced.

For large (R ≫ l_d) asteroids, the Yarkovsky effect is basically determined by the thermal parameter Θ. In the above calculations, the K₁ = 0.0015 W m⁻¹ K⁻¹ with other parameters given in Table 1 corresponds to a relatively small thermal parameter Θ₁ = 0.736, although Θ is typically of the order of unity (Vokrouhlický et al. 2015) for objects found near Earth or in the main belt. The bulk density and heat capacity of asteroids may vary over a range of multiples, while the value of conductivity, K, has a much greater range of more than four orders of magnitude (Delbo et al. 2015). To investigate how the heat conduction may influence the relationship between da/dt and S₂, using the same shape models and the same ρ and C values as in Fig. 2, we calculated the Yarkovsky effect of asteroids with two additional conductivity values K₂ = 0.15 and K₃ = 15 W m⁻¹ K⁻¹, corresponding to thermal parameter Θ₂ = 7.36 and Θ₃ = 73.6, and check whether the linear relationship between the strength of Yarkovsky effect and the shape index can be verified. We note that these parameters adopted here including the Θ may be not typical values for real asteroids. But our goal is to verify the applicability of the shape index, even under very extreme conditions. In this sense, these parameters are acceptable.

To reach the dynamic equilibrium of temperature distribution on the surface, 100 and 150 rotation periods are required in the numerical simulations for the cases of K₂ and K₃, respectively. Limited by our computational resource, this time we calculated only 34 asteroid models spinning all around the z axis in these simulations. The results, including those for K₁, are summarised in Fig. 3.

As K increases by orders of magnitude, Θ deviates away from 1, and for a given shape (thus a given S₂), the semimajor axis drift rate decreases. For example, at S₂ = 0.95, we obtain from Fig. 3 that the da/dt is 0.01827, 0.01176, and 0.001513 AU/Myr for Θ₁, Θ₂, and Θ₃, respectively. In fact, a relative larger Θ (> 1) indicates that the energy is being efficiently stored in the surface layer and the temperature is more likely to be evenly distributed on the surface, resulting in a weaker Yarkovsky effect. However, as shown in Fig. 3, the linear relationship between da/dt and S₂ holds very well for all three Θ values. The least-squares fitting gives the relations as: $$\frac{\mathrm{d}a}{\mathrm{d}t}=\{\begin{array}{cc}0.01804S_2+0.001130& \text{ for }{\mathrm{\Theta }}_1,\\ 0.01371S_2-0.001264& \text{ for }{\mathrm{\Theta }}_2,\\ 0.001773S_2-0.0001709& \text{ for }{\mathrm{\Theta }}_3.\end{array}$$(13)

Note that all units (AU/Myr) have been omitted for simplicity in the equations. The coefficient of determination, R², and mean standard deviation, σ_ϵ, are 0.992 and 1.5 × 10⁻⁴ AU/Myr for Θ₂. For the case of Θ₃, they are R² = 0.993 and σ_ϵ = 1.8 × 10⁻⁵ AU/Myr, respectively. In addition, comparing the da/dt calculated using Eq. (13) with the real values obtained from COMSOL simulations, we find that the maximum deviation is always less than 3% for all these 34 models.

Encouraged by these very good R² values and small deviations σ_ϵ for cases with very different K (thus different Θ) values, we calculate the Yarkovsky effect for more cases with K in a wide range, and fit the results as the linear relation: $$\mathrm{d}a/\mathrm{d}t=k{S}_2+b.$$(14)

To save computations, for each K value we selected from the 34 samples arbitrarily 5 models (with shape indices S₂ ≈ 0.805, 0.833, 0.962, 0.993, 1.05), calculated their surface temperature using COMSOL, derived the da/dt, and finally fitted the results as Eq. (14). The slope, k, and the intercept, b, are summarized in Fig. 4. We note that the k and b for cases of K₁, K₂, and K₃ (thus Θ₁, Θ₂, and Θ₃) in Fig. 4 were re-computed from the same 5 shape models.

Ideally, if the Yarkovsky effect depends on the shape of an asteroid characterised by the shape index, we would expect an intercept b = 0 in Eq. (14). But some approximations have been unavoidably introduced in the definition of the shape index, for instance, as shown in Fig. 1, the delay angles of different surface elements on a body are not strictly identical. In addition, a 3D model is used in calculating the Yarkovsky effect, but the lateral heat conduction is ignored in the definition of the shape index. However, as shown in Fig. 3, a linear relationship is very well obeyed for shape models with the shape index S₂ from 0.8 to 1.5. And, as we can see in Fig. 4 and Eq. (13), the intercept is at least one order of magnitude smaller than the slope, which is the exact da/dt at S₂ = 1. Therefore, for practical use, the linear relationship is still an excellent approximation.

In Fig. 4, as the conductivity (thus Θ) increases, the slope k of the linear fit first increases and then decreases, reaching its maximum when Θ ~ 1.74. In fact, Θ is the ratio of the thermal relaxation time to the rotational period (Farinella et al. 1998). Theoretically, for Θ < 1, energy is rapidly lost by radiation, and the delay angle is small, while for Θ > 1, the temperature gradient on the surface cannot be efficiently established. Both lead to a smaller Yarkovsky effect.

When the slope, k, reaches its maximum, we notice the intercept, b, approaches to zero (Fig. 4), implying the optimal condition for the definition of the shape index is attained.

Fig. 3 Semi-major axis drift rate da/dt for different thermal parameter Θ versus shape index S2. The orange dots, green squares, and cyan diamonds represent the da/dt for asteroids with different K values, thus of Θ1 = 0.736, Θ2 = 7.36, and Θ3 = 73.6, respectively. The lines are the corresponding least-squares fits.

Fig. 4 Parameters (k and b in Eq. (14)) of linear fitting for different conductivities K (and thus thermal parameter values Θ). Note the errors for large K values are so small that the error bars are narrower than the symbols standing for data points.

3.4 Lateral heat conduction

The rate of heat conduction inside a body is determined not only by the value of K, but also by the temperature gradient. A higher Θ leads to a lower temperature difference between adjacent surface elements. Thus, as we have shown in Figs. 3 and 4, a larger Θ value does not necessarily fail the mechanism of the linear relationship between da/dt and S₂. We take the asteroid 1998 KY26 as a shape example to show the variation of temperature distribution on the surface at different Θ values. Again, three K values (thus three thermal parameters Θ₁, Θ₂, Θ₃) are adopted and the results are presented in Fig. 5.

The small K value (corresponding to Θ₁, the top panels of Fig. 5) prohibits the diffusion of energy in the body, and the surface elements around the equator absorb and keep much more Solar radiation than the ones in the polar region. As a result, the global temperature difference across the surface reaches nearly 200 K; that is, the temperature gradient on the surface is large. But the thermal conductivity is so small that the lateral heat transfer is nearly blocked. Thus the assumption of null lateral heat transfer in the definition of the shape index is fulfilled, and the linear dependence of Yarkovsky effect on the shape index is soundly obeyed.

When K increases by 4 orders of magnitude to K₃ (thus Θ increases by 2 orders of magnitude to Θ₃) the global temperature difference across the surface is only ~40 K (bottom panels of Fig. 5). The small temperature gradient limits the lateral heat transfer, and therefore, the vast majority of heat received by a surface element is in fact stored in and then emitted after a delay from the same surface element. The assumption of negligible lateral diffusion of heat is still available, and the linear relationship between the semi-major axis drift and S₂ still holds.

Comparing the surface temperatures in the three cases in Fig. 5, we find that the most obvious change due to different Θ values is the difference between the temperatures of day side and night side. For any given Θ value, the greatest temperature difference on the surface is always between the equatorial and polar regions. As the thermal conductivity increases, the magnitude of temperature variation along the latitude circles (as an example, the temperature along the equator is plotted in Fig. 6) decreases much more than it does along the longitude circles. This again indicates that the lateral heat transfer contributes little to the temperature distribution, otherwise the polar region should obtain more energy from lateral diffusion of heat to reach a higher temperature and a lower temperature variation along the longitude circle. In addition, as Θ increases, the temperature distribution on the surface, particularly along the east-west direction, maintains nearly the same pattern in the three cases (Fig. 5), implying also that the heat diffusion is enhanced mainly in the direction into the deeper part inside the asteroid, which leads to a longer delay of thermal emission and a lower temperature fluctuation on the surface.

Finally, the influence of the lateral heat transfer can also be estimated by comparing the Yarkovsky effects of ideal spherical asteroid calculated using both 1D and 3D models. In the 1D model, the lateral heat transfer is completely neglected, while in the 3D model it is fully taken into account. Regarding the Yarkovsky effect (da/dt) calculated from the 3D model as accurate, we find the relative errors of the 1D model are 1.6 × 10⁻⁴, 3.2 × 10⁻⁴, and 2.3 × 10⁻³ when the conductivity is set as K₁, K₂, and K₃, respectively. Although the relative error (due to neglecting the lateral heat transfer) increases somewhat as the thermal conductivity greatly increases (by four orders of magnitude), it remains at a low level of 10⁻³. This suggests that the lateral heat transfer does not influence the Yarkovsky effect much, and a high K will not invalidate the linear relationship between the semimajor axis drift rate (da/dt) and the shape index (S₂), as long as the size of asteroid is not too small (R ≫ l_d).

Fig. 5 Temperature distribution on the surface of asteroids with different thermal parameter Θ, taking the shape of 1998 KY 26 as an example. The small body has a retrograde rotation around the z axis, with the Solar radiation from the direction nin = (1, 0, 0). From top to bottom are the cases for Θ1, Θ2, and Θ3, and from left to right are the isometric views, day hemispheres and night hemispheres. Note that the colour bars have different temperature scales.

Fig. 6 Surface temperature on equator of the asteroid model in Fig. 5. The longitude of 0° (midnight) and 180° (noon) are on the meridian plane. Three K values (corresponding to Θ1, Θ2, Θ3) are adopted and the results are given in different colours.

3.5 Scattering effect and self-heating effect

In addition to absorbing the Solar radiation, a surface element may also obtain energy from other elements by their scattering (the Solar radiation) and thermal irradiation (this effect is called self-heating). This part of energy is ignored in the definition of the shape index S₂.

The contribution from the scattering and self-heating effects can be numerically estimated. We select the four shape models in Fig. 2c as examples, since they deviate significantly from fully convex shapes and the influences of scattering and self-heating effects are expected to be most evident in these models. The ‘real’ Yarkovsky effects of these four models have been calculated and presented in Fig. 3, where the scattering and self-heating effects were taken into account. We regard the tangential component of the recoil force (F_t in Eq. (4)) as the accurate measurement of the Yarkovsky effect, and numerically re-calculate the force component, $F_t^{\mathrm{\prime }}$, for the same four models but without the scattering and self-heating effects (by turning off these effects in the COMSOL simulations). The relative differences between them $\left(\right|1-F_t^{\mathrm{\prime }}/F_t\left|\right)$ are summarized in Table 2.

As shown in Table 2, the relative error due to neglecting the scattering and self-heating effects increases as K increases. However, even for the highest possible conductivity K₃ and for the most irregular shapes, the error induced by neglecting scattering and self-heating is still very small, only of the order of 10⁻³.

The scattering Solar radiation absorbed by the i-th surface element at moment t can be written as: $$E_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}}=-\alpha (1-\alpha )\mathcal{E}\sum _{j\ne i}{f}_{i,j}(1-s_j(t))({\mathbf{n}}_j\cdot {\mathbf{n}}_{\mathrm{i}\mathrm{n}}).$$(15)

The f_i,j in Eq. (15) is the view factor from the i-th surface element to the j-th one, and it reads: $$f_{i,j}=v_{i,j}\frac{\cos {\theta }_i\cos {\theta }_j}{\pi d_{i,j}^2}{a}_j,$$(16)

where v_i,j denotes whether there is line-of-sight visibility between the i-th and j-th surface elements, d_i,j is the distance between them, θ_i and θ_j are the incidence and emission angles of the i-th and j-th surface element, respectively, and a_j is the area. For more details on the definition and calculation of view factor, see Lagerros (1998) and Rozitis & Green (2011).

The self-heating effect depends on the surface temperature distribution in addition to the surface geometry of the asteroid. For the i-th surface element, the radiation that it absorbs from other elements is: $$E_{\mathrm{r}\mathrm{a}\mathrm{d}}=\alpha ϵ\sigma \sum _{j\ne i}{f}_{i,j}{T}_j^4(t).$$(17)

The temperature T_j at the j-th surface element can be estimated by assuming an instantaneous heat absorption and emission; that is, an equilibrium: $$-\alpha \mathcal{E}(1-s_j(t))({\mathbf{n}}_j\cdot {\mathbf{n}}_{\mathrm{i}\mathrm{n}})=ϵ\sigma T_j^4(t).$$(18)

All the energy absorbed by the i-th surface element at time t, including contributions from the direct Solar radiation as in Eq. (6) (taking into account the shadowing effect), scattering as in Eq. (15), and self-heating as in Eq. (17), now is: $$\begin{array}{rl}{E}_{\text{in }}+E_{\text{scat }}+E_{\text{rad }}=& -\alpha \mathcal{E}[a_i(1-s_i(t))({\mathbf{n}}_i\cdot {\mathbf{n}}_{\text{in }})\\ & +\sum _{j\ne i}{f}_{i,j}(1-s_j(t))({\mathbf{n}}_j\cdot {\mathbf{n}}_{\text{in }})].\end{array}$$(19)

Finally, an improved shape index S₃ that takes into account the scattering effect and self-heating effect can be written as: $$\begin{array}{rl}{S}_3=\frac{A}{P}{\int }_0^P\{& \sum _i^N({\mathbf{n}}_i\cdot {\mathbf{n}}_{\mathrm{i}\mathrm{n}})[a_i(1-s_i(t))({\mathbf{n}}_i\cdot {\mathbf{n}}_{\mathrm{i}\mathrm{n}})\\ & -\sum _{j\ne i}{f}_{i,j}(1-s_j(t))({\mathbf{n}}_j\cdot {\mathbf{n}}_{\mathrm{i}\mathrm{n}})\left]\right\}\mathrm{d}t\end{array}$$(20)

To check whether this new shape index S₃ improves the estimation of Yarkovsky effect, we compute the S₃ values for the 34 shape models and calculate the linear fitting between da/dt and S₃ as we did for S₁ and S₂. The coefficients of determination (R²) of the fitting for both S₂ and S₃ are listed in Table 3.

The smaller R² values for S₃ than the ones for S₂ (Table 3) reveal that this new shape index S₃ does not make a better linear relationship. In fact, when we calculate the self-heating effect, E_rad, the surface temperature has been obtained by simply assuming an instantaneous equilibrium (Eq. (18)). This approximation is equivalent to the assumption K = 0, and it might introduce extra error. This explains that the index S₃ behaves even worse when a higher Θ value is adopted, as shown in Table 3. Without knowing the specific temperature distribution, it is difficult to evaluate the influence of self-heating on the Yarkovsky effect.

Considering the fact that the errors caused by the scattering and self-heating effects are ignorable in practice (Table 2), we therefore abandon the idea of incorporating these effects into the definition of the shape index.

3.6 Shape index and shapes

For a given asteroid, the Yarkovsky effect is determined by its thermal parameter, Θ. When focusing on the shape of an asteroid with given thermal parameter, two factors may influence the Yarkovsky effect. One is the orientation of surface elements, and the other is the effect of shadowing. For the former, if the normal of a facet on the body surface is parallel to the equatorial plane, as the body rotates this facet will have the largest cross sectional area with respect both to the Solar radiation and to the tangential direction of the orbit, and it will produce the most effective Yarkovsky force. Therefore, the more a body has surface elements with normal being parallel to the equatorial plane, the stronger the Yarkovsky effect is, and according to the definition of the shape index, the bigger S₁ and S₂ are. For example, a cylinder with radius R and height h rotating around its axis, according to Eq. (10), has a shape index S₁ = (3h/4R)^1/3. The higher the cylinder is, the bigger S₁ is, and vice versa. On the contrary, if a body is oblate, and most of its surface elements have normal nearly perpendicular to the equatorial plane, its shape index must be small. A very low (flat) cylinder is an extreme example.

To show more clearly the shape index value of different shapes, we consider triaxial ellipsoids with three semi-axes (a, b, c). Two parameters e₁, e₂ describing the shape are defined as: $$e_1=\frac{c}{\sqrt{ab}},e_2=\frac{b}{a}.$$(21)

Apparently, a triaxial ellipsoid reduces to a biaxial ellipsoid when e₂ = 1, and ellipsoids characterised by e₂ and 1/e₂ are identical to each other. Assuming the ellipsoid rotates around the axis with semi-major axis c, we calculate the S₁ of them and present the results in Fig. 7. We note that ellipsoids are fully convex, so S₁ = S₂ here.

For biaxial ellipsoids, the shape index increases as the oblateness e₁ increases. S₁ < 1 for oblate spheroid with e₁ < 1, and on the contrary S₁ > 1 for e₁ > 1. For triaxial ellipsoids, the shape index S₁ > 1 for all the ellipsoids with e₁ > 1, indicating that all ‘tall’ spheroids have stronger Yarkovsky effect than a sphere. While for oblate triaxial ellipsoids with e₁ < 1, the shape index can also be larger than 1, especially when e₂ is small. This phenomenon of oblate triaxial ellipsoids having stronger Yarkovsky effect has already been observed in the numerical simulations and explained in Xu et al. (2022). It should be noted that the rotation of an ellipsoid with e₁ > 1 around the axis c is generally unstable.

Asteroids 1994 KW4 Beta and Geographos (Fig. 8) are good examples to show the intuitive connection between the shape index and their shape. Approximately regard them as ellipsoids, the 1994 KW4 Beta with dimensions (0.571 × 0.463 × 0.349) km has e₁ = 0.68, e₂ = 0.81, and Geographos of (5.0 × 2.0 × 2.1) km has e₁ = 0.66, e₂ = 0.4 (Ostro et al. 2006; Hudson & Ostro 1999). They have similar e₁ but quite different e₂ values. The shape indices of these two ellipsoids are 0.811 and 0.954. The values calculated from the real shape models are S₁ = 0.82 (S₂ = 0.81) for 1994 KW 4 Beta and S₁ = 1.10 (S₂ = 1.06) for Geographos, respectively. Because 1994 KW4 Beta has a relatively smooth surface, the ‘real’ shape index S₁ of it roughly agrees with the ellipsoid approximation (0.82 ≈ 0.811). But, more irregularities can be seen on Geographos’ surface, and they introduce a bigger ‘real’ shape index than the ellipsoid approximation (1.10 > 0.954).

Top-shaped asteroids are also quite common, with Ryugu and Bennu being the ready examples. A top-shaped asteroid can be reduced to a combination of two cones with a common base and oppositely directed tops. Denote the ratio between the equatorial radius R and height h of the northern and southern cones by r, that is, r = R/h, we can write the shape index S₁ as: $$S_1(r)=\frac{3}{4}{\left(\frac{2}{r}\right)}^{1/3}{\left(\frac{1}{1+r^2}\right)}^{1/2}.$$(22)

Clearly, S₁ decreases monotonically as r increases. We find that S₁(r) = 1 when r = 0.56, implying that this top-shaped asteroid has equivalent Yarkovsky effect as a sphere asteroid when the cone angle is ~30°. Otherwise, when the cone angle is smaller than 30° (r < 0.56), the object will have a stronger Yarkovsky effect and the shape index S₁ > 1. For Ryugu and Bennu, their equatorial to polar axis ratios of 1.147 and 1.104 (Watanabe et al. 2019; Lauretta et al. 2019) give S₁ = 0.593, 0.614, respectively. These values are much smaller than their actual shape indices S₁ = 0.955, 0.960 calculated from their real shapes, implying that the jointed-cones is not a good model for these two asteroids. In fact, some ‘plains’ can be found in the equatorial region of both asteroids, and they might contribute considerably to the real shape index.

So far, the cylinders, ellipsoids, and spinning top shape discussed above are all fully convex, which means S₁ = S₂ for them. In the real world, the boulders and craters on the surface of asteroids may increase the specific surface area, alter the orientation of surface elements, and produce shadows. Although the influences of these irregularities on shape index are highly complex, a general estimation can still be made by approximating the real shape as a deformed version of its inertia ellipsoid. For the 34 shape models, we calculate the shape index both for real shapes (S₂) and for their inertia ellipsoids ($S_2^{\mathrm{\prime }}$), and find that for almost every model the S₂ of the real shape is larger than that of its inertia ellipsoid. Especially for those of very irregular shapes, their S₂ are significantly larger than the corresponding $S_2^{\mathrm{\prime }}$. Again, taking the four asteroids in Fig. 3 as examples, we listed in Table 4 the two shape indices. We note the difference between S₂ and $S_2^{\mathrm{\prime }}$ for Kleopatra is over 20%.

Roughly speaking, data in Table 4 implies that the irregularities in shape are likely to increase the shape index, and a stronger Yarkovsky effect might be expected if the asteroid is very irregularly shaped.

Fig. 7 Shape index S1 (given by colour) for triaxial ellipsoids characterised by oblateness e1 and e2. The dashed line is for S1 = 1.

Fig. 8 Shapes of asteroids (a) 1994 KW4 Beta and (b) Geographos.

4 Conclusion

The Yarkovsky effect plays an important role in the orbital evolution of asteroids. In addition to the thermal, physical, and dynamical parameters, the Yarkovsky effect of an asteroid strongly depends on its shape, which brings great complexity to both the analytical and numerical calculations of the effect. To obtain a quick estimation of the diurnal Yarkovsky effect, Xu et al. (2022) proposed the idea of ‘effective area’, which can be easily calculated if the shape of an asteroid is known. In this paper, we improve the idea and refine an index that can quantitatively measure the influence of asteroids’ shape on the Yarkovsky effect.

We adopt 34 real asteroids with shape data available as the sample of shape models, rescale the sizes of these models so that they all have the same volume as a sphere with a radius of 10 meters, assign them with specific thermal parameters, and put them all in a circular orbit with semi-major axis of 1 AU. We established the numerical model in the COMSOL software package to calculate the temperature on their surface, from which the Yarkovsky force and then the semi-major axis drift rate (da/dt) can be derived. The migration rate da/dt was used to measure the strength of the Yarkovsky effect.

After revisiting the definition of effective area in Xu et al. (2022), which is roughly a double projection of the surface area onto a plane parallel to the spin axis, we normalised it to a shape index S₁ (Eq. (10)). When defining S₁, we have assumed that all energy absorbed by a surface element will be released by thermal radiation after a delay angle θ from the same surface element. For a given asteroid, the Yarkovsky effect is determined by the dimensionless thermal parameter Θ. Although the delay θ can be calculated if Θ is known, we proved that the shape index S₁ does not depend on the delay θ specifically. We also showed by numerical simulations that the assumption of all surface elements having the same delay angle is valid (Fig. 1). Therefore, this shape index S₁ is determined solely by the shape of an asteroid, and a linear relationship between da/dt and S₁ as proposed in Xu et al. (2022) was reaffirmed using a larger sample of asteroids (Fig. 2).

We found that the outliers in the fitted linear relationship are those that might be most influenced by the shadowing effect. We therefore improved the shape index by removing the contributions from surface elements when they are under the projected shadows. The linear relationship between da/dt and the improved shape index S₂ (Eq. (12)) is much better than the former one with respect to S₁, and the outliers are tamed (Fig. 2). As for the linear fit, the coefficient of determination (R²) increases from 0.984 for S₁ to 0.995 for S₂.

When defining the shape index, a 1D heat transfer model is assumed and we totally ignore the lateral heat transfer. This assumption might be violated if the asteroid has a large thermal conductivity. We adopted three conductivity values in a wide range over four orders of magnitude (K₁ = 0.0015, K₂ = 0.15, K₃ = 15 W m⁻¹ K⁻¹), and tested the performance of shape index S₂. We found that the strength of Yarkovsky effect changes remarkably as the thermal parameter Θ varies for different K, but the linear relationship between da/dt and S₂ is robust (Fig. 3).

We extended the calculations and obtained the fitting parameters (the slope k and intercept b) of the linear relationship for other Θ values. As Θ increases, the slope k increases first and then decreases, reaching its maximum when Θ ~ 1.74. The linear relationship is very well kept for all Θ values. After these calculations, if the thermal parameter Θ is known, the Yarkovsky effect of an asteroid of any shape can be directly estimated by its shape index S₂ using the linear relationship with the parameters given in Fig. 4.

The rate of heat transfer inside an asteroid depends on both the thermal conductivity and temperature gradient. An increase in K (thus in Θ) results in quicker heat diffusion inside the body and more evenly distributed temperature over the surface (Fig. 5) that often leads to a weaker Yarkovsky effect. Furthermore, a higher Θ reduces the temperature gradient along the east-west direction much more than it does along the north-south direction (Figs. 5, 6). We argue that the lateral heat conduction inside an asteroid is not significant, even when it has a high thermal conductivity. This was also confirmed by the negligible difference of Yarkovsky effect between a 1D model (without the lateral transfer) and a 3D model (with the lateral transfer) for the same spherical body. Therefore, the estimation of Yarkovsky effect by its linear dependence on the shape index S₂ is valid for all bodies satisfying R ≫ l_d.

The radiation scattering and the self-heating effects may affect the Yarkovsky effect and their influences may be enhanced when the thermal conductivity increases. However, through numerical simulations, we found that the contributions from these effects to the Yarkovsky effect are ignorable. Even for very irregularly shaped asteroids with very large K₃ = 15 W m⁻¹ K⁻¹ (thus large Θ = 73.6) value, their contributions makes only a relative difference in the Yarkovsky force of ~10⁻³ (Table 2). After incorporating these effects into the definition of the shape index (now S₃), we found no improvements to the linear relationship between da/dt and S₃ (Table 3). In fact, without knowing the details of surface temperature, the influence of self-heating cannot be included precisely in the estimation. Therefore, we conclude that it is not necessary to take into account the scattering and self-heating effects in our definition of shape index.

To demonstrate the dependence of shape index on specific shapes, we computed the S₁ and S₂ for several sets of simple shape models, including cylinders, ellipsoids (Fig. 7) and top-shaped asteroids (Eq. (22)), and found that prolate shapes, such as ellipsoids with oblateness e₁ > 1 and spinning tops with equatorial to polar axis ratio r < 0.56 would have large shape index S₁ > 1, while the oblate shapes have smaller S₁. In addition, for irregularly shaped asteroids, we calculated the S₂ both for their real shapes and for their inertia ellipsoids (Table 4), and found that irregularities in shape generally increase the shape index S₂.

As we have shown, the shape index (S₂) is a good indicator of the strength of diurnal Yarkovsky effect. It is solely determined by the shape of an asteroid, and can be easily computed, bringing great convenience to the estimations of Yarkovsky effect for a large number of asteroids. Knowing the shape of an asteroid by observations, the Yarkovsky effect can be easily obtained by comparing its shape index with the one of an ‘ideal standard asteroid’ (a spherical asteroid) with the same thermal parameters. Meanwhile, the shape of an asteroid is generally inverted from its light curves. With the shape index as a bridge, we may be able to find a direct link between the Yarkovsky effect of an asteroid and its light curves, on which we are working through machine learning now.

Acknowledgements

Our great appreciation goes to Prof. Vokrouhlicky for his insightful comments that helped us in improving our manuscript. This work has been supported by the National Natural Science Foundation of China (NSFC, Grants No. 12373081 and No. 12150009) and the China Manned Space Program with grant No.CMS-CSST-2025-A16. We also acknowledge the support from National Key R&D Program of China (2019YFA0706601).

References

All Tables

All Figures

Fig. 1 Directions and magnitudes of the averaged radiation forces of surface elements on sphere model and YORP model. A thermal conductivity K1 = 0.0015 W m−1 K−1 (correspondingly Θ = 0.736) is adopted for both models. The direction and magnitude have been projected in the orbital (equatorial) plane here (see text). (a) The magnitudes of radiation forces of surface elements versus the delay angle (upper panel), and the frequency of delay angles of all surface elements (lower panel). The delay angle of a surface element is the angle between the radiation force averaged over a rotation period and the corresponding averaged force of the same surface element but with K = 0. Blue and orange indicate the results for sphere and YORP model, respectively. (b) The direction and magnitude of the radiation force at different latitudes on the surface of a sphere. The longitude 0° is at noon. The longitude of the black dot at certain latitude indicates that the radiation force at this latitude points to the direction given by the longitude. The size of black dots represents the magnitude of the radiation force (in N/m2, see text) projected in the orbital plane. The yellow circle indicates the direction and magnitude (F/4πR2) of the total radiation force. The background colour gives the surface temperature.

In the text

	Fig. 2 Relationship between the semi-major axis drift rate (da/dt) and the shape index (S1,2). (a) da/dt versus S1 (cyan diamonds) and S2 (orange dots). (b) Fitting residuals. (c) Four shape models with large residuals: YORP, Kleopatra, Ida, and Golevka.
In the text

	Fig. 3 Semi-major axis drift rate da/dt for different thermal parameter Θ versus shape index S2. The orange dots, green squares, and cyan diamonds represent the da/dt for asteroids with different K values, thus of Θ1 = 0.736, Θ2 = 7.36, and Θ3 = 73.6, respectively. The lines are the corresponding least-squares fits.
In the text

	Fig. 4 Parameters (k and b in Eq. (14)) of linear fitting for different conductivities K (and thus thermal parameter values Θ). Note the errors for large K values are so small that the error bars are narrower than the symbols standing for data points.
In the text

Fig. 5 Temperature distribution on the surface of asteroids with different thermal parameter Θ, taking the shape of 1998 KY 26 as an example. The small body has a retrograde rotation around the z axis, with the Solar radiation from the direction nin = (1, 0, 0). From top to bottom are the cases for Θ1, Θ2, and Θ3, and from left to right are the isometric views, day hemispheres and night hemispheres. Note that the colour bars have different temperature scales.

In the text

	Fig. 6 Surface temperature on equator of the asteroid model in Fig. 5. The longitude of 0° (midnight) and 180° (noon) are on the meridian plane. Three K values (corresponding to Θ1, Θ2, Θ3) are adopted and the results are given in different colours.
In the text

	Fig. 7 Shape index S1 (given by colour) for triaxial ellipsoids characterised by oblateness e1 and e2. The dashed line is for S1 = 1.
In the text

	Fig. 8 Shapes of asteroids (a) 1994 KW4 Beta and (b) Geographos.
In the text